Using IDL: Mathematics |
Numerical methods of approximating integrals are important in many areas of pure and applied science. For a function of a single variable, f (x), it is often the case that the antiderivative F = Ú f (x) dx is unavailable using standard techniques such as trigonometric substitutions and integration-by-parts formulas. These standard techniques become increasingly unusable when integrating multivariate functions,
f (x, y) and f (x, y, z). Numerically approximating the integral operator provides the only method of solution when the antiderivative is not explicitly available. IDL offers the following numerical methods for the integration of uni-, bi-, and trivariate functions:
Note IDL's iterated Gaussian Quadrature routines, INT_2D and INT_3D, follow the dy-dx and dz-dy-dx order of evaluation, respectively. Problems not conforming to this standard must be changed as described in the following example. |
Suppose that we wish to evaluate
The order of integration is initially described as a dx-dy region in the x-y plane. Using the diagram below, you can easily change the integration order to dy-dx.
The integral is now of the form
The new expression can be evaluated using the INT_2D function.
To use INT_2D, we must specify the function to be integrated and expressions for the upper and lower limits of integration. First, we write an IDL function for the integrand, the function f (x, y):
FUNCTION fxy, X, Y RETURN, Y * COS(X^5) END
Next, we write a function for the limits of integration of the inner integral. Note that the limits of the outer integral are specified numerically, in vector form, while the limits of the inner integral must be specified as an IDL function even if they are constants. In this case, the function is:
FUNCTION pq_limits, X RETURN, [0.0, X^2] END
Now we can use the following IDL commands to print the value of the integral expressed above. First, we define a variable AB_LIMITS containing the vector of lower and upper limits of the outer integral. Next, we call INT_2D. The first argument is the name of the IDL function that represents the integrand (FXY, in this case). The second argument is the name of the variable containing the vector of limits for the outer integral (AB_LIMITS, in this case). The third argument is the name of the IDL function defining the lower and upper limits of the inside integral (PQ_LIMITS, in this case). The fourth argument (48) refers to the number of transformation points used in the computation. As a general rule, the number of transformation points used with iterated Gaussian Quadrature should increase as the integrand becomes more oscillatory or the region of integration becomes more irregular.
ab_limits = [0.0, 2.0] PRINT, INT_2D('fxy', ab_limits, 'pq_limits', 48)
IDL prints:
0.055142668
This is the exact solution to 9 decimal accuracy.
Suppose that we wish to evaluate
This integral can be evaluated using the INT_3D function. As with INT_2D, we must specify the function to be integrated and expressions for the upper and lower limits of integration. Note that in this case IDL functions must be provided for the upper and lower integration limits of both inside integrals.
For the above integral, the required functions are the integrand f (x, y, z):
FUNCTION fxyz, X, Y, Z RETURN, Z * (X^2 + Y^2 + Z^2)^1.5 END
The limits of integration of the first inside integral:
FUNCTION pq_limits, X RETURN, [-SQRT(4.0 - X^2), SQRT(4.0 -X^2)] END
The limits of integration of the second inside integral:
FUNCTION uv_limits, X, Y RETURN, [0.0, SQRT(4.0 - X^2 - Y^2)] END
We can use the following IDL commands to determine the value of the above integral using 6, 10, 20 and 48 transformation points.
For 6 transformation points:
PRINT, INT_3D('fxyz', [-2.0, 2.0], $ 'pq_limits', 'uv_limits', 6)
IDL prints:
57.417720
For 10 transformation points:
PRINT, INT_3D('fxyz', [-2.0, 2.0], $ 'pq_limits', 'uv_limits', 10)
IDL prints:
57.444248
20 transformation points:
PRINT, INT_3D('fxyz', [-2.0, 2.0], $ 'pq_limits', 'uv_limits', 20)
IDL prints:
57.446201
48 transformation points:
PRINT, INT_3D('fxyz', [-2.0, 2.0], $ 'pq_limits', 'uv_limits', 48)
IDL prints:
57.446265
The exact solution to 6-decimal accuracy is 57.446267.
See Differentiation and Integration (in the functional category Mathematics) for a brief description of IDL routines for differentiation and integration. Detailed information is available in the IDL Reference Guide.
IDL Online Help (March 06, 2007)