2

I have a function to compute Gaussian quadrature of a function $f(x)$ over the region $x \in [a,b]$. Here, $f(x)$ takes only one argument. I would want to know what a good practice would be to use gaussquad with a function which might take more arguments, for example getlaser below.

Laser is a derived type, and calling gaussquad(mylaser%getlaser, a, b) obviously does not work.

double precision function gaussquad(f, a, b) result(I)
  implicit none
  double precision :: f
  double precision, intent(in) :: a, b
  
  I = 2.d0*f(b-a)
end function

double precision function getlaser(this, t)
  implicit none
  class(Laser), intent(in) :: this
  double precision, intent(in) :: t

  getlaser = dsin(this%omega*t)
end function getlaser
mbertolino
  • 53
  • 3

1 Answers1

1

The getlaser procedure has a user-defined dummy argument this which makes it impossible to define a general integration module. In the following I will explain how to implement such a general integration module assuming standard data types.

One option would be to define an optional parameter array in gaussquad which can be passed through to the procedure f.

Following is a possible implementation for the integration module

! integ.f90
module integ_m
  implicit none

  private
  public gaussquad

  abstract interface
    real function finter(x, p)
      real,           intent(in) :: x
      real, optional, intent(in) :: p(:)
    end function
  end interface

contains

  function gaussquad(f, a, b, p) result(int)
    !! compute integral: int_a^b f(x; p) dx
    procedure(finter)          :: f
      !! function to integrate
    real,           intent(in) :: a, b
      !! integration bounds
    real, optional, intent(in) :: p(:)
      !! parameter array
    real                       :: int
      !! integral value

    int = (b-a)*f(0.5*(a+b), p=p)
  end function

end module

One would use it like in this program

! main.f90
program main
  use integ_m, only: gaussquad

  implicit none

  print *, 'integrate x^2',         gaussquad(parabola, 0.0, 1.0        )
  print *, 'integrate laser (sin)', gaussquad(getlaser, 0.0, 1.0, [10.0])


contains

  real function parabola(x, p)
    real,           intent(in) :: x
    real, optional, intent(in) :: p(:)

    if (present(p)) error stop "function doesnt use parameters"
    parabola = x*x
  end function

  real function getlaser(t, p)
    real,           intent(in) :: t
    real, optional, intent(in) :: p(:)

    getlaser = sin(p(1)*t)
  end function

end program

Compilation and running yields

$ gfortran -g -Wall -fcheck=all integ.f90 main.f90 && ./a.out 
 integrate x^2  0.250000000    
 integrate laser (sin) -0.958924294
jack
  • 1,658
  • 1
  • 7
  • 18