具有更多参数和集成的功能

Function with more arguments and integration

我有一个简单的问题,但我无法在任何地方找到解决方案。 我必须集成一个函数(例如使用辛普森规则子例程)但我不得不向我的函数传递不止一个参数:一个是我稍后要集成的变量,另一个只是来自不同的值我无法在函数内部执行的计算。

问题是辛普森子程序只接受 f(x) 来执行积分而不接受 f(x,y)。

根据弗拉基米尔的建议,我修改了代码。

下面的例子:

Program main2
!------------------------------------------------------------------
! Integration of a function using Simpson rule 
! with doubling number of intervals
!------------------------------------------------------------------
! to compile:
! gfortran main2.f90 -o simp2

implicit none
double precision r, rb, rmin, rmax, rstep, integral, eps
double precision F_int
integer nint, i, rbins
double precision t

rbins = 4
rmin = 0.0
rmax = 4.0
rstep = (rmax-rmin)/rbins
rb = rmin

eps = 1.0e-8
func = 0.0

t=2.0


do i=1,rbins
   call func(rb,t,res)
   write(*,*)'r, f(rb) (in main) = ', rb, res
   !test = F_int(rb)
   !write(*,*)'test F_int (in loop) = ', test    
   call simpson2(F_int(rb),rmin,rb,eps,integral,nint)     
   write(*,*)'r, integral = ', rb, integral
   rb = rb+rstep
end do

end program main2

subroutine func(x,y,res)
!----------------------------------------
! Real Function 
!----------------------------------------
implicit none
double precision res
double precision, intent(in) ::  x
double precision y
res = 2.0*x + y 
write(*,*)'f(x,y) (in func) = ',res
return
end subroutine func


function F_int(x)
!Function to integrate

implicit none
double precision F_int, res
double precision, intent(in) ::  x
double precision y
call func(x,y,res)
F_int = res
end function F_int



Subroutine simpson2(f,a,b,eps,integral,nint)
!==========================================================
! Integration of f(x) on [a,b]
! Method: Simpson rule with doubling number of intervals  
!         till  error = coeff*|I_n - I_2n| < eps
! written by: Alex Godunov (October 2009)
!----------------------------------------------------------
! IN:
! f   - Function to integrate (supplied by a user)
! a   - Lower limit of integration
! b   - Upper limit of integration
! eps - tolerance
! OUT:
! integral - Result of integration
! nint     - number of intervals to achieve accuracy
!==========================================================
implicit none
double precision f, a, b, eps, integral
double precision sn, s2n, h, x
integer nint
double precision, parameter :: coeff = 1.0/15.0 ! error estimate coeff
integer, parameter :: nmax=1048576              ! max number of intervals
integer n, i



! evaluate integral for 2 intervals (three points)
h = (b-a)/2.0
sn = (1.0/3.0)*h*(f(a)+4.0*f(a+h)+f(b))

write(*,*)'a, b, h, sn (in simp) = ', a, b, h, sn 

! loop over number of intervals (starting from 4 intervals)
n=4
do while (n <= nmax)
s2n = 0.0   
h = (b-a)/dfloat(n)
do i=2, n-2, 2
  x   = a+dfloat(i)*h
  s2n = s2n + 2.0*f(x) + 4.0*f(x+h)
end do
s2n = (s2n + f(a) + f(b) + 4.0*f(a+h))*h/3.0
if(coeff*abs(s2n-sn) <= eps) then
  integral = s2n + coeff*(s2n-sn)
  nint = n
  exit
end if
sn = s2n
n = n*2
end do
return
end subroutine simpson2

我想我已经很接近解决方案了,但我想不出来... 如果我调用 simpson2(F_int, ..) 而没有将参数放入 F_int 我收到此消息:

call simpson2(F_int,rmin,rb,eps,integral,nint)
              1
Warning: Expected a procedure for argument 'f' at (1)

有什么帮助吗? 提前致谢!

现在您有了我们可以使用的代码,干得好!

你需要告诉编译器,F_int 是一个函数。这可以通过

来完成
external F_int

但是最好学习 Fortran 90 并使用模块或至少使用接口块。

module my_functions

  implicit none

contains


subroutine func(x,y,res)
!----------------------------------------
! Real Function 
!----------------------------------------
implicit none
double precision res
double precision, intent(in) ::  x
double precision y
res = 2.0*x + y 
write(*,*)'f(x,y) (in func) = ',res
return
end subroutine func


function F_int(x)
!Function to integrate

implicit none
double precision F_int, res
double precision, intent(in) ::  x
double precision y
call func(x,y,res)
F_int = res
end function F_int

end module

现在您可以轻松使用模块并集成功能

use my_functions

call simpson2(F_int,rmin,rb,eps,integral,nint)

但是你会发现F_int还不知道y是什么!它有自己的 y 和未定义的值!您应该将 y 放入模块中,以便每个人都能看到它。

module my_functions

  implicit none

  double precision :: y

contains

不要忘记删除 y 的所有其他声明!在函数 F_int 和主程序中。或许换个称呼也比较好。

不要忘记在主循环中的某处设置 y 的值!