时间戳

首页 > 解决方案 > How to fix solution with double integration

fortran - How to fix solution with double integration

问题描述

The solution for this double integration is -0.083 but in the final compliation it appears -Infinity. It seems that the error is very simple, but I really can't find it.

I have been searching specially in the module section but I don't see why it appears like -Infinity. For example, if you change the two functions between them (x in f2 and x^2 in f1) the solution for the integration is 0.083 and the code gives it correct. Can annyone find the error? Thanks a lot.

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

The main program is simple, just calling the subroutine and using the module. It also prints the final result.

标签： fortrangfortranfortran90fortran95

解决方案

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

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

  • real*8是非标准的 Fortran。改为使用real(real64)。real64 必须由use :: iso_fotran_env, only: real64.
  • do-loops 中的非整数表达式 ( do i=1d0,n-1d0) 是现代 Fortran 中已删除的功能。请改用整数。
  • 代码应使用空格和缩进进行格式化
  • print*,应该替换为write(*,*)
  • 代码应始终使用英文名称
  • 写implicit none在模块的开头，而不是每个函数。
  • 使用语句private、public和使模块/程序接口清晰only
  • 如果要转换为 type 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

推荐阅读