fortran - Function with more arguments and integration

Question

I have I simple problem but I cannot find a solution anywhere. I have to integrate a function (for example using a Simpson's rule subroutine) but I am obliged to pass to my function more than one argument: one is the variable that I want to integrate later and another one is just a value coming from a different calculation which I cannot perform inside the function.

The problem is that the Simpson subroutine only accept f(x) to perform the integral and not f(x,y).

After Vladimir suggestions I modified the code.

Below the example:

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

I think I'm pretty close to the solution but I cannot figure it out... If I call simpson2(F_int, ..) without putting the argument in F_int I receive this message:

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

Any help? Thanks in advance!

Answer 1

Now you have a code we can work with, good job!

You need to tell the compiler, that F_int is a function. That can be done by

external F_int

but it is much better to learn Fortran 90 and use modules or at least interface blocks.

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

Now you can easily use the module and integrate the function

use my_functions

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

But you will find that F_int still does not know what y is! It has it's own y with undefined value! You should put y into the module instead so that everyone can see it.

module my_functions

  implicit none

  double precision :: y

contains

Don't forget to remove all other declarations of y! Both in function F_int and in the main program. Probably it is also better to call it differently.

Don't forget to set the value of y somewhere inside your main loop!