Using Triple Integration in IDL

Calling the QROMB, QROMO, or QSIMP routine within a user-supplied function of one of these same routines allows the user to calculate double and triple integrals. This Help Article contains an example demonstrating triple integration.

Each term of an integrand (the equation within the integral) can contain another integration method. The double and triple integrations are performed over each term of the integrand. The following example uses triple integration to determine the volume of a three-dimensional equation over a specified region.

Integrating Over a Solid (Triple Integration): this example calculates the volume of the solid defined by



over the region defined by 0 <= x <= 1, 0 <= y <= 1, and 0 <= z <= 1. The solution can be found by evaluating the triple integral



Note that the correct solution to this integral is 3.

This example program is made up of six routines: the main routine, the integration in the z-direction, the second integration of the xy term, the second integration of the second x^2*y^2 term, the third integration in the x term, and the third integration in the x^2 term. The main routine is the last routine in the program. The file containing this program should be named the same as the main routine (integratingOverAVolume.pro).

Code example:

FUNCTION xSquaredTerm, x
 
   ; Integration of the x squared term.
   thirdIntegration = 9.*x^2
   RETURN, thirdIntegration
 
END
 
FUNCTION xTerm, x
 
   ; Integration of the linear x term.
   thirdIntegration = x
   RETURN, thirdIntegration
 
END
 
FUNCTION xSquaredySquaredTerm, y
 
   ; Integration of the y squared term.
   secondIntegration = QROMB('xSquaredTerm', 0., 1.)*y^2
   RETURN, secondIntegration
 
END
 
FUNCTION xyTerm, y
 
   ; Integration of the linear y term.
   secondIntegration = QROMB('xTerm', 0., 1.)*y
   RETURN, secondIntegration
 
END
 
FUNCTION zDirection, z
 
   ; Re-write equation to consider both x terms.
   firstIntegration = QROMB('xSquaredySquaredTerm', 0., 1.) + $
   8.*(QROMB('xyTerm', 0., 1.))*z + 1.
   RETURN, firstIntegration
 
END
 
PRO integratingOverAVolume
 
   ; Determine the area of the surface represented
   ; by 9x^2y^2 + 8xyz + 1.
   volume = QROMB('zDirection', 0., 1. )
 
   ; Output results.
   PRINT, 'Resulting volume:  ', volume
 
END

_________________________________________
Reviewed by BC on 09/05/2014