The QSIMP function performs numerical integration of a function over the closed interval [*A, B*] using Simpsonâ€™s rule.

## Examples

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

; Define lower limit of integration:

A = 0.0

; Define upper limit of integration:

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`

## Syntax

*Result* = QSIMP( *Func*, *A*, *B* [, /DOUBLE] [, EPS=*value*] [, JMAX=*value*] )

## Return Value

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.

## 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} - 2x^{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`

**Note: **If QSIMP is complex then only the real part is used for the computation.

### 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 x 10^{-6}. For double-precision calculations, the default value is 1.0 x 10^{-12}.

### JMAX

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

## Version History

4.0 |
Introduced |

## Resources and References

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.