How do I use nested integral functions to find a double integral in fortran?

Question

I am trying to find the volume under a shape defined by h(x,y)

I have the following function which uses the Simpson Rule to integrate a one-dimensional function

!----------------------------------------------------------------------------
double precision function INTEGRAL2(FUNC,a,b,N)
!----------------------------------------------------------------------------
! ***  numerical integration (Simpson-rule) with equidistant spacing      ***
!----------------------------------------------------------------------------
  implicit none
  double precision,external :: FUNC                    ! the function to be integrated
  double precision,intent(in) :: a,b                   ! boundary values
  integer,intent(in) :: N                    ! number of sub-intervals
  double precision :: dx,x1,x2,xm,f1,f2,fm,int         ! local variables
  integer :: i
  dx  = (b-a)/DBLE(N)                        ! x subinterval
  x1  = a                                    ! left   
  f1  = FUNC(a)
  int = 0.d0
  do i=1,N
    x2  = a+DBLE(i)*dx                       ! right
    xm  = 0.5d0*(x1+x2)                      ! midpoint
    f2  = FUNC(x2)
    fm  = FUNC(xm)
    int = int + (f1+4.d0*fm+f2)/6.d0*(x2-x1) ! Simpson rule
    x1  = x2
    f1  = f2                                 ! save for next subinterval
  enddo
  INTEGRAL2 = int
end

How do I use this to find the double integral of h(x,y)? Note h(x,y) is not reducible to h(x)*h(y), and I want to use this function nested, rather than modifying it/writing a function to do double integration.

I'm fundamentally stuck on the concept of writing the code, but I suspect using a module will be crucial.

Answer 1

Think how Simpson's rule can be derived in 1-d. You divide the domain into pairs of intervals, fit a quadratic curve to each set of 3 points and integrate that. Well, you can do the same in 2-d (or higher) - fit a bi-quadratic by Lagrange interpolation and integrate that; the weights turn out to be the same in each direction - [1,4,1]dx/3 for a 2.dx interval - as for Simpson's rule.

function func( x, y )                            ! function to be integrated
   use iso_fortran_env
   implicit none
   real(real64) func
   real(real64), intent(in) :: x, y

   func = exp( -0.5 * ( x ** 2 + y ** 2 ) )      ! bi-gaussian - should integrate to 2.pi for large domains

end function func

!=======================================================================

function integrate2d( f, x1, x2, y1, y2, nx, ny ) result( res )
   use iso_fortran_env
   implicit none
   real(real64), external :: f                   ! a function of the 2 variables x and y
   real(real64) res
   real(real64), intent(in) :: x1, x2, y1, y2    ! limits of integration
   integer, intent(in) :: nx, ny                 ! numbers of intervals (MUST be EVEN)
   real(real64) dx, dy                           ! interval widths
   real(real64), dimension(-1:1), parameter :: wt = [ 1.0_real64 / 3, 4.0_real64 / 3, 1.0_real64 / 3 ]
                                                 ! weights in each direction (same as Simpson's rule)
   integer i, j, p, q
   real(real64) x, y

   dx = ( x2 - x1 ) / nx
   dy = ( y2 - y1 ) / ny

   res = 0.0
   do j = 1, ny - 1, 2
      y = y1 + j * dy
      do i = 1, nx - 1, 2
         ! Consider a local rectangle 2.dx by 2.dy
         x = x1 + i * dx
         do p = -1, 1
            do q = -1, 1
               res = res + wt(p) * wt(q) * f( x + p * dx, y + q * dy )
            end do
         end do
      end do
   end do

   res = res * dx * dy

end function integrate2d

!=======================================================================

program main
   use iso_fortran_env
   implicit none

   integer nx, ny
   real(real64) :: x1, x2, y1, y2
   real(real64), external :: func
   real(real64), external :: integrate2d
 
   nx = 100;   ny = 100
   x1 = -10.0;   x2 = 10.0
   y1 = -10.0;   y2 = 10.0

   print *, "Integral value is ", integrate2d( func, x1, x2, y1, y2, nx, ny )

end program main

!=======================================================================

Incidentally, there are better forms of quadrature (e.g. gaussian quadrature) that you might like to look up.