3D integration

Average rating
(1 vote)

Igor provides a function "Integrate1D()" for integrating functions of one-dimensional arguments. The associated help file gives a method for using it to integrate functions of two-dimensional arguments, with the statement "This method can be extended to higher dimensions." I found the effort to be not entirely trivial, and thought to save others the work with the following code. I have shown an interestingly complex user function of (x,y,z) that has a simple analytic result for large integration limits after converting to a radial coordinate system

#pragma rtGlobals=1		// Use modern global access method.
 
Function do3dIntegration(xmin,xmax,ymin,ymax,zmin,zmax)
	Variable xmin,xmax,ymin,ymax,zmin,zmax	
	Variable/G globalXmin=xmin
	Variable/G globalXmax=xmax
	Variable/G globalYmin=ymin
	Variable/G globalYmax=ymax
	Variable/G globalY, globalZ
	return  Integrate1D(userFunction3,zmin,zmax,1)		// Romberg integration
End
 
Function userFunction3(inZ)
	Variable inZ	
	NVAR globalZ=globalZ
	globalZ=inZ
	NVAR globalYmin=globalYmin
	NVAR globalYmax=globalYmax	
	return integrate1D(userFunction2,globalYmin,globalYmax,1)	// Romberg integration
End
 
Function userFunction2(inY)
	Variable inY	
	NVAR globalY=globalY
	globalY=inY
	NVAR globalXmin=globalXmin
	NVAR globalXmax=globalXmax	
	return integrate1D(userFunction1,globalXmin,globalXmax,1)	// Romberg integration
End
 
Function userFunction1(inX)
	Variable inX	
	NVAR globalY=globalY
	NVAR globalZ=globalZ	
	return ((inX^2+globalY^2+globalZ^2)^(3/2)) * exp(-(inX^2+globalY^2+globalZ^2))
//  r^3 exp(-r^2)  --> 2*pi * u^2 *exp(-u) --> 4 pi integral
//  r exp(-r^2)     --> 2*pi * u    *exp(-u) --> 2 pi integral
End

printf "%8.6f\r", do3dIntegration(-5,5,-5,5,-5,5)
  12.566372

The exact result should be 4 *pi.

Back to top