##
QSIMP

The
QSIMP
function performs
numerical integration of a function over the closed interval [*
A, B*
] using
Simpson's rule. The result will have the same structure as the smaller of *
A*
and *
B*
, and the resulting type will be single- or double-precision floating, depending on the input types.

QSIMP is based on the routine ```
qsimp
```

described in section 4.2 of *
Numerical Recipes in C: The Art of Scientific Computing*
(Second Edition), published by Cambridge University Press, and is used by permission.

###
Calling Sequence

Result = QSIMP(*
Func, A, B*
)

###
Arguments

####
Func

A scalar string specifying the name of a user-supplied IDL function to be integrated. This function must accept a single scalar argument *
X*
and return a scalar result. It must be defined over the closed interval [*
A, B*
].

For example, if we wish to integrate the fourth-order polynomial

*
y*
= (*
x*
^{
4}
- 2*
x*
^{
2}
) sin(*
x*
)

we define a function SIMPSON to express this relationship in the IDL language:

FUNCTION simpson, X

RETURN, (X^4 - 2.0 * X^2) * SIN(X)

END

####
A

The lower limit of the integration. *
A*
can be either a scalar or an array.

####
B

The upper limit of the integration. *
B*
can be either a scalar or an array.

Note: If arrays are specified for *
A*
and *
B*
, then QSIMP integrates the user-supplied function over the interval [*
A*
_{
i}
, *
B*
_{
i}
] for each *
i*
. If either *
A*
or *
B*
is a scalar and the other an array, the scalar is paired with each array element in turn.

###
Keywords

####
DOUBLE

Set this keyword to force the computation to be done in double-precision arithmetic.

####
EPS

The desired fractional accuracy. For single-precision calculations, the default value is 1.0 *
¥*
10^{
-6}
. For double-precision calculations, the default value is 1.0 *
¥*
10^{
-12}
.

####
JMAX

2^{
(JMAX - 1)}
is the maximum allowed number of steps. If not specified, a default of 20 is used.

###
Example

To integrate the SIMPSON function (listed above) over the interval [0, *
p*
/2] and print the result:

A = 0.0

B = !PI/2.0

PRINT, QSIMP('simpson', A, B)

IDL prints:

-0.479158

The exact solution can be found using the integration-by-parts formula:

FB = 4.*B*(B^2-7.)*SIN(B) - (B^4-14.*B^2+28.)*COS(B)

FA = 4.*A*(A^2-7.)*SIN(A) - (A^4-14.*A^2+28.)*COS(A)

exact = FB - FA

PRINT, exact

IDL prints:

-0.479156