如何通过双重集成修复解决方案

Question

这种双重积分的解决方案是 -0.083，但在最终的编译中它显示为 -Infinity。貌似错误很简单，但是实在找不到。

我一直在专门搜索模块部分，但我不明白为什么它看起来像 -Infinity。例如，如果您更改它们之间的两个函数（f2 中的 x 和 f1 中的 x^2），则积分的解为 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

首先，就像 cmets 中提到的那样：您的问题不清楚。您使用了哪些输入参数a、b 和n，您期望得到哪个结果？

除此之外：您发布的代码使用了已弃用的功能和非标准类型以及错误的代码样式。 一些一般提示：

• real*8 是非标准 Fortran。请改用real(real64)。 real64 必须由use :: iso_fotran_env, only: real64. 导入
• do 循环中的非整数表达式 (do i=1d0,n-1d0) 是现代 Fortran 中已删除的功能。请改用整数。
• 代码应使用空格和缩进进行格式化
• print*, 应替换为 write(*,*)
• 代码应始终使用英文名称
• 在模块开头写implicit none，而不是每个函数。
• 使用private、public和only语句使模块/程序界面清晰
• 如果你想转换成real，使用函数REAL而不是DBLE
• 我更喜欢使用result 的更简洁的函数定义
• 使用 intent 关键字：intent(in) 将变量作为 const 引用传递。
• 主程序中的变量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