Как исправить решение с двойной интеграцией

Question

Решением для этой двойной интеграции является -0,083, но в конечном компиляции оно появляется -Infinity.Кажется, что ошибка очень проста, но я действительно не могу ее найти.

Я специально искал в разделе модуля, но не понимаю, почему она выглядит как -Infinity.Например, если вы измените две функции между ними (x в f2 и x ^ 2 в f1), решение для интегрирования будет 0.083, и код даст это правильно.Может ли кто-нибудь найти ошибку?Большое спасибо.

module funciones

contains

function f(x,y)

implicit none

real*8:: x,y,f

f=2d0*x*y

end function

function f1(x)

real*8::x,f1

f1=x

end function


function f2(x)

real*8::x,f2

f2=x**2d0

end function

function g(x,c,d,h)

implicit none

integer::m,j

real*8::x,y,c,d,k,s,h,g

m=nint(((d-c)/h)+1d0)

k=(d-c)/dble(m)

s=0.

do j=1d0,m-1d0

y=c+dble(j)*k

s=s+f(x,y)

end do

g=k*(0.5d0*(f(x,c)+f(x,d))+s)

return

end function


subroutine trapecio(a,b,n,integral)

implicit none

integer::n,i

real*8::a,b,c,d,x,h,s,a1,a2,b1,b2,integral

h=(b-a)/dble(n)

s=0d0

do i=1d0,n-1d0

x=a+dble(i)*h

c=f1(x)

d=f2(x)

s=s+g(x,c,d,h)

end do


a1=f1(a)

a2=f2(a)

b1=f1(b)

b2=f2(b)

integral=h*(0.5d0*g(a,a1,a2,h)+0.5d0*g(b,b1,b2,h)+s)

end subroutine

end module



program main

use funciones

implicit none

integer::n,i

real*8::a,b,c,d,x,s,h,integral

print*, "introduzca los valores de a, b y n"

read(*,*) a, b, n

call trapecio (a,b,n,integral)

print*,integral

end program

Основная программа проста, просто вызывая подпрограмму и используя модуль.Он также печатает окончательный результат.

Answer 1

Прежде всего, как упомянуто в комментариях: ваша проблема не ясна.Какие входные параметры a, b и n вы используете и какого результата ожидаете?

Кроме этого: в опубликованном вами коде используются устаревшие функции, нестандартные типы и стиль плохого кода,Некоторые общие советы:

• real*8 - это нестандартный фортран.Вместо этого используйте real(real64).real64 должен быть импортирован с помощью use :: iso_fotran_env, only: real64.
• нецелых выражений (do i=1d0,n-1d0) в циклах do - удаленная функция в современном Фортране.Вместо этого используйте целые числа.
• код должен быть отформатирован с пробелами и отступами
• print*, должен быть заменен на write(*,*)
• код должен всегда использовать английские имена
• писать implicit none в начале модуля, не для каждой функции.
• сделать интерфейс модуля / программы понятным, используя операторы private, public и only
• Если вы хотите преобразовать в тип real, используйте функцию REAL вместо DBLE
• Я предпочитаю определение функции очистки, используя result
• use intent ключевые слова: intent(in) передает переменную как постоянную ссылку.
• переменные c,d,x,s,h в основной программе не используются.Скомпилируйте с предупреждениями для обнаружения неиспользуемых переменных.

Этот код был изменен в соответствии с предложенными мною предложениями:

module funciones
use :: iso_fortran_env, only: real64
implicit none

private
public :: trapecio, r8

   integer, parameter :: r8 = real64

contains
   function f(x,y) result(value)
      real(r8), intent(in) :: x,y
      real(r8) :: value

      value = 2._r8*x*y
   end function

   function f1(x) result(value)
      real(r8), intent(in) :: x
      real(r8) :: value

      value = x
   end function

   function f2(x) result(value)
      real(r8), intent(in) :: x
      real(r8) :: value

      value = x**2._r8
   end function

   function g(x,c,d,h) result(value)
      real(r8), intent(in) :: x, c, d, h
      real(r8) :: value

      real(r8) :: y, k, s
      integer :: m, j

      m = NINT(((d-c)/h)+1._r8)
      k = (d-c)/REAL(m, r8)
      s = 0._r8
      do j = 1, m-1
         y = c + REAL(j,r8)*k
         s = s + f(x,y)
      end do

      value = k*(0.5_r8*(f(x,c)+f(x,d))+s)
   end function

   subroutine trapecio(a, b, n, integral)
      real(r8), intent(in) :: a, b
      integer, intent(in) :: n
      real(r8), intent(out) :: integral

      integer :: i
      real(r8) :: c, d, x, h, s, a1, a2, b1, b2
      h = (b-a)/REAL(n,r8)
      s = 0._r8

      do i = 1, n-1
         x = a + REAL(i,r8)*h
         c = f1(x)
         d = f2(x)
         s = s + g(x,c,d,h)
      end do

      a1 = f1(a)
      a2 = f2(a)
      b1 = f1(b)
      b2 = f2(b)
      integral = h*(0.5_r8*g(a,a1,a2,h) + 0.5_r8*g(b,b1,b2,h) + s)
   end subroutine
end module

program main
   use funciones, only: trapecio, r8

   implicit none

   integer :: n,i
   real(r8) :: a,b,integral

   write(*,*) "introduzca los valores de a, b y n"
   read(*,*) a, b, n
   call trapecio (a,b,n,integral)
   write(*,*) integral
end program