Documentation

# `numeric`::`quadrature`

MuPAD® notebooks will be removed in a future release. Use MATLAB® live scripts instead.

MATLAB live scripts support most MuPAD functionality, though there are some differences. For more information, see Convert MuPAD Notebooks to MATLAB Live Scripts.

## Syntax

```numeric::quadrature(f(`x`), `x = a .. b`, <`GaussLegendre = n | GaussTschebyscheff = n | NewtonCotes = n`>, <`Adaptive = v`>, <`MaxCalls = m`>)
```

## Description

`numeric::quadrature(f(x), x = a..b)` computes a numerical approximation of .

`numeric::quadrature` returns itself symbolically if the integrand `f(x)` contains symbolic objects apart from the integration variable `x` that cannot be converted to numerical values via `float`. Symbolic objects such as π or etc. are accepted.

The integrand `f(x)` should be integrable in the Riemann sense. In particular, `f(x)` should be bounded on the integration interval `x = a..b`. Certain types of mild singularities are handled, but numerical convergence is not guaranteed and will be slow in most cases. Also discontinuities and singularities of (higher) derivatives of `f(x)` slow down numerical convergence. For integrands with irregular points, it is recommended to split the integration into several parts, using subintervals on which the integrand is smooth. Cf. Example 4.

Integrands given by user-defined procedures can be handled. See Example 4 and Example 5.

`numeric::quadrature` returns itself symbolically if the boundaries `a,b` contain symbolic objects that cannot be converted to numerical values via `float`. Symbolic objects such as π or etc. as well as `infinity` and `-infinity` are accepted.

### Note

For infinite ranges, the user should make sure that the integral exists! If `f(x)` does not decay as fast as at infinity, then convergence may be slow.

Boundaries a > b are accepted, using .

For complex values of `a,b`, the integration is to be understood as a contour integral along a straight line from `a` to `b`. Complex boundaries must not involve `infinity`.

Multi-dimensional integration such as

```numeric::quadrature ( numeric::quadrature(f(x,y), y = A(x)..B(x)), x = a..b)```

is possible. See Example 3 and Example 5.

Internally, an adaptive mechanism based on a fixed quadrature rule specified by `method = n` is used. It tries to keep the relative quadrature error of the result below . The last digit(s) of the result may be incorrect due to roundoff effects.

The routine `numeric::quadrature` is purely numerical: no symbolic check for singularities etc. is carried out.

## Environment Interactions

The function is sensitive to the environment variable `DIGITS`, which determines the numerical working precision.

## Examples

### Example 1

We demonstrate the standard calls for adaptive numerical integration:

`numeric::quadrature(exp(x^2), x = -1..1)`
` `
`numeric::quadrature(max(1/10, cos(PI*x)), x = -2..0.0123)`
` `
`numeric::quadrature(exp(-x^2), x = -2..infinity)`
` `

The precision goal is set by `DIGITS`:

```DIGITS := 50: numeric::quadrature(besselJ(0, x), x = 0..PI)```
` `

Note that due to the internal adaptive mechanism, the choice of different methods should not have any significant effect on the quadrature result:

```DIGITS := 10: numeric::quadrature(sin(x)/x, x = -1..10, GaussLegendre = 5), numeric::quadrature(sin(x)/x, x = -1..10, GaussLegendre = 160), numeric::quadrature(sin(x)/x, x = -1..10, NewtonCotes = 8)```
` `

### Example 2

The user should make sure that the integrand is well defined and sufficiently regular. The following fails, because Newton-Cotes quadrature tries to evaluate the integrand at the boundaries:

`numeric::quadrature(sin(x)/x, x = 0..1, NewtonCotes = 8)`
```Error: Division by zero. [_power] Evaluating: Quadsum ```

One may cure this problem be assigning a value to `f(0)`. The integrand is passed to the integrator as `hold(f)` to prevent premature evaluation of `f(x)` to `sin(x)/x`. Internally, `numeric::quadrature` replaces `x` by numerical values and then evaluates the integrand. Note that one has to define `f(0.0) := 1` rather than ```f(0) := 1```:

```f := x -> sin(x)/x: f(0.0) := 1: numeric::quadrature(hold(f)(x), x = 0..1, NewtonCotes = 8)```
` `

The default method (Gauss-Legendre quadrature) does not evaluate the integrand at the end points:

`numeric::quadrature(sin(x)/x, x = 0..1)`
` `

Nevertheless, problems may still arise for improper integrals:

`numeric::quadrature(ln((1 + x^4)^2 - 1), x = 0 .. 1)`
```Warning: Precision goal not achieved after 10000 function calls. Increase 'MaxCalls' and try again for a more accurate result. [numeric::quadrature] ```
` `

In this example, the integrand is evaluated close to 0. The expression (1 + x4)2 - 1 suffers from severe numerical cancellation and is dominated by round-off. A numerically stable representation is (1 + x4)2 - 1 = x4 (x4 + 2):

`numeric::quadrature(ln(x^4*(x^4 + 2)), x = 0..1)`
` `

Note that convergence is rather slow, because the integrand is not bounded.

`delete f:`

### Example 3

```Q := numeric::quadrature: Q(Q(exp(x*y), x = 0..y), y = 0..1)```
` `

Also more complex types of nested calls are possible. We use numerically defined functions

`b := y -> Q(exp(-t^2-t*y), t = y..infinity):`

and

`f := y -> cos(y^2) + Q(exp(x*y), x = 0..b(y)):`

`Q(f(y), y = 0..1)`
` `

Multi dimensional quadrature is computationally expensive. Note that a call to `numeric::quadrature` evaluates the integrand at least 3 n times, where n is the number of nodes of the internal quadrature rule (by default, n = 20 for DIGITS ≤ 10). The following triple quadrature would call the exp function no less than (3 20)3 = 216000 times!

`Q(Q(Q(exp(x*y*z), x = 0..y+z), y = 0..z), z = 0..1)`

For low precision goals, low order quadrature rules suffice. In the following, we reduce the computational costs by using Gauss-Legendre quadrature with 5 nodes. We use the shorthand notation `GL` to specify the `GaussLegendre` method:

```DIGITS := 4: Q(Q(Q(exp(x*y*z), x=0..y+z, GL=5), y=0..z, GL=5), z=0..1, GL=5)```
` `
`delete Q, b, f, DIGITS:`

### Example 4

We demonstrate how integrands given by user-defined procedures should be handled. The following integrand

```f := proc(x) begin if x<1 then sin(x^2) else cos(x^5) end_if end_proc:```

cannot be called with a symbolic argument:

`f(x)`
```Error: Unable to evaluate to Boolean. [_less] Evaluating: f ```

Consequently, one must use `hold` to prevent premature evaluation of `f(x)`:

`numeric::quadrature(hold(f)(x), x = -1..PI/2)`
` `

Note that the above integrand is discontinuous at x = 1, causing slow convergence of the numerical quadrature. It is much more efficient to split the integral into two subquadratures with smooth integrands:

```numeric::quadrature(sin(x^2), x = -1..1) + numeric::quadrature(cos(x^5), x = 1..PI/2)```
` `

See Example 5 for multi-dimensional quadrature of user-defined procedures.

`delete f:`

### Example 5

The following integrand

```f := proc(x, y) begin if x<y then x-y else (x-y) + (x-y)^5 end_if end_proc:```

can only be called with numerical arguments and must be delayed twice for 2-dimensional quadrature:

`Q := numeric::quadrature:`
`Q(Q(hold(hold(f))(x, y), x = 0..1), y = 0..1)`
` `

Note that delaying the integrand via `hold` will not work for triple or higher-dimensional quadrature! The user can handle this by making sure that the integrand can also be evaluated for symbolic arguments:

```f := proc(x, y, z) begin if not testtype([args()], Type::ListOf(Type::Numeric)) then return(procname(args())) end_if; if x^2 + y^2 + z^2 <= 1 then return(1) else return(0) end_if end_proc:```

Note that this function is not continuous, implying slow convergence of the numerical quadrature. For this reason, we use a low precision goal of only 2 digits and reduce the costs by using a low order quadrature rule:

```DIGITS := 2: Q(Q(Q(f(x, y, z), x=0..1, GL=5), y=0..1, GL=5), z=0..1, GL=5)```
` `
`delete f, Q, DIGITS:`

### Example 6

The following example uses non-adaptive Gauss-Tschebyscheff quadrature with an increasing number of nodes. The results converge quickly to the exact value:

```a := exp(x)/sqrt(1 - x^2), x = -1..1: numeric::quadrature(a, Adaptive = FALSE, GT = n) \$ n = 3..7```
` `
`delete a:`

### Example 7

The improper integral exists. Numerical convergence, however, is rather slow because of the singularity at x = 0:

`numeric::quadrature(x^(-9/10), x = 0..1)`
```Warning: Precision goal not achieved after 10000 function calls. Increase 'MaxCalls' and try again for a more accurate result. [numeric::quadrature] ```
` `

We remove the limit for the number of function calls and let `numeric::quadrature` grind along until a result is found. The time for the computation grows accordingly, the last digit is incorrect due to roundoff effects:

`numeric::quadrature(x^(-9/10), x = 0..1, MaxCalls = infinity)`
` `

### Example 8

The following integral does not exist in the Riemann sense, because the integrand is not bounded:

`numeric::quadrature(1/x, x = -1..2)`
```Warning: Precision goal not achieved after 10000 function calls. Increase 'MaxCalls' and try again for a more accurate result. [numeric::quadrature] ```
` `

We increase `MaxCalls`. There is no convergence of the numerical algorithm, because the integral does not exist. Consequently, some internal problem must arise: `numeric::quadrature` reaches its maximal recursive depth:

`numeric::quadrature(1/x, x = -1..2, MaxCalls = infinity)`
```Warning: Precision goal not achieved after 'MAXDEPTH=500' recursive calls. There might be a singularity of '1/x' close to 'x=3.910318545e-148'. Increase 'MAXDEPTH' and try again for a more accurate result. [adaptiveQuad] ```
` `

## Parameters

 `f(x)` `x` An identifier or an indexed identifier `a`, `b` Real or complex numerical values or ## Options

`GaussLegendre`, `GaussTschebyscheff`, `NewtonCotes`

Options, specified as `GaussLegendre = n`, ```GaussTschebyscheff = n```, `NewtonCotes = n`

The name of the underlying quadrature scheme. Each quadrature rule can be used with an arbitrary number of nodes specified by the positive integer n.

Usually there is no need to use this option to change the default method `GaussLegendre = n` with `n = 20,40,80` or `160`, depending on the precision goal determined by the environment variable `DIGITS`. Due to the corresponding high quadrature orders 40, 80, 160 or 320, respectively, the default settings are suitable for high precision computations.

With `GaussLegendre = n`, an adaptive version of Gauss-Legendre quadrature with n nodes is used.

For DIGITS ≤ 200, the weights and abscissae of Gaussian quadrature with n `=` 20, 40, 80, and 160 nodes are available and the integration starts immediately.

For DIGITS > 200, the routine `numeric::gldata` is called to compute the Gaussian data before the actual integration starts. This will be noted by some delay in the first call of `numeric::quadrature`.

For `DIGITS` much larger than 200, it is recommended not to use the default setting but to use ```GaussLegendre = n``` with sufficiently high `n` instead. A reasonable choice is nDIGITS.

With `GaussTschebyscheff = n`, non-adaptive Gauss-Tschebyscheff quadrature with n nodes is used. This method may only be used in conjunction with ```Adaptive = FALSE```.

### Note

With `NewtonCotes = n`, an adaptive version of Newton-Cotes quadrature with n nodes is used. Note that these quadrature rules become numerically unstable for large n (n much larger than 10).

Following alternative names for the methods are accepted:

`GaussLegendre`, `Gauss`, `GL`,

`GaussTschebyscheff`, `GT`, `GaussChebyshev`, `GC`,

`NewtonCotes`, `NC`.

`Adaptive`

Option, specified as `Adaptive = v`

`v` may be `TRUE` or `FALSE`. With ```Adaptive = FALSE```, the internal error control is switched off.

The default setting is `Adaptive = TRUE`. An adaptive quadrature scheme with automatic control of the quadrature error is used.

The call ```numeric::quadrature(f(x), x = a..b, method = n, Adaptive = FALSE)``` returns the quadrature sum approximating without any control of the quadrature error. The weights Bi and abscissae Ci are determined by the quadrature rule given by `method = n`. For the methods `GaussLegendre`, `GaussTschebyscheff` and `NewtonCotes`, these data are available via `numeric::gldata`, `numeric::gtdata` and `numeric::ncdata`, respectively.

`Adaptive = FALSE` may only be used in conjunction with `method = n`.

Usually there is no need to switch off the internal adaptive quadrature via `Adaptive = FALSE`. This option is meant to give easy access to standard non-adaptive quadrature rules of Gauss-Legendre, Gauss-Tschebyscheff and Newton-Cotes type, respectively. The user may want to build his own adaptive quadrature based on these non-adaptive rules. Cf. Example 6.

`MaxCalls`

Option, specified as `MaxCalls = m`

`m` must be a (large) positive integer or `infinity`. It is the maximal number of evaluations of the integrand, before `numeric::quadrature` gives up.

When called interactively, `numeric::quadrature` returns the approximation it has computed so far and issues a warning. See Example 7. When called from inside a procedure, `numeric::quadrature` returns `FAIL`.

The default value is `m = MAXDEPTH*n`. The environment variable `MAXDEPTH` (default value 500) represents the maximal recursive depth of a function. `n` is the number of nodes of the basic quadrature rule specified by the optional argument `method = n`. If no such method is specified by the user, then Gauss-Legendre quadrature is used with n = 20 for DIGITS ≤ 10, n = 40 for 10 < DIGITS ≤ 50, n = 80 for 50 < DIGITS ≤ 100, n = 160 for 100 < DIGITS. Thus, for DIGITS = 10, the default setting is `MaxCalls = 10000`.

The default value of `m` ensures that the `MaxCalls` limit is reached before `numeric::quadrature` reaches its maximal internal recursive depth. Specifying `MaxCalls = infinity` removes this restriction and `numeric::quadrature` computes until it obtains an approximation with about `DIGITS` correct digits or until it runs into an internal error. The typical error that may occur is the evaluation of the integrand at a singularity. There also is the danger of `numeric::quadrature` reaching its maximal internal recursive depth. When called interactively, `numeric::quadrature` returns the approximation it has computed so far and issues a warning. See Example 8. When called from inside a procedure, `numeric::quadrature` returns `FAIL`.

## Return Values

Floating point number is returned, unless non-numerical symbolic objects in the integrand `f(x)` or in the boundaries `a,b` prevent numerical evaluation. In this case, `numeric::quadrature` returns itself symbolically. If numerical problems occur, then `FAIL` is returned.

#### Mathematical Modeling with Symbolic Math Toolbox

Get examples and videos