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.
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.