## `numeric::quadrature` -- numerical integration

#### Introduction

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

#### Call(s)

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

#### Parameters

`f(x)` | - | expression in x
`x` | - | identifier or indexed identifier
`a,b` | - | real or complex numerical values or +-`infinity`

#### Options

``method=n`` | - | `method` is the name of the underlying quadrature scheme, either GaussLegendre, GaussTschebyscheff, or NewtonCotes. Each quadrature rule can be used with an arbitray number of nodes specified by the positive integer n.
`Adaptive=`v | - | `v` may be `TRUE` or `FALSE`. With `Adaptive=FALSE` the internal error control is switched off.
`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.

#### Returns

A 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 unevaluated. If numerical problems occur, then `FAIL` is returned.

#### Side Effects

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

#### Related Functions

`int`, `numeric::gldata`, `numeric::gtdata`, `numeric::int`, `numeric::ncdata`

#### Details

- `numeric::quadrature` returns itself unevaluated, 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 `PI` or `sqrt(2)` 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. Cf. Examples 4 and 5.
- `numeric::quadrature` returns itself unevaluated, if the boundaries `a,b` contain symbolic objects that cannot be converted to numerical values via `float`. Symbolic objects such as `PI` or `sqrt(2)` etc. as well as `infinity` and `-infinity` are accepted.

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

- Boundaries a>b are accepted. Note:

      quadrature(f(x), x=a..b) = - quadrature(f(x), x=b..a).

- 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. Cf. examples 3 and 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 10^(-DIGITS). 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.

#### Option: method = 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`>>200 it is recommended not to use the default setting but to use GaussLegendre=`n` with sufficiently high `n` (approximately `DIGITS`) instead.

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

- 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 (>>10).

- Following alternative names for the methods are accepted:

  GaussLegendre, Gauss, GL,

  GaussTschebyscheff, GT, GaussChebyshev, GC,

  NewtonCotes, NC.

#### Option: Adaptive= v

- 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 (b-a)*sum(B[i]*f(a+C[i]*(b-a)), i=1..n approximating int(f(x),x=a..b) without any control of the quadrature error. The weights B[i] and abscissae C[i] 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.

#### Option: MaxCalls= m

- `numeric::quadrature` gives up after `m` evaluations of the integrand. When called interactively, `numeric::quadrature` returns the approximation it has computed so far and issues a warning. Cf. 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. Cf. example 8. When called from inside a procedure, `numeric::quadrature` returns `FAIL`.

#### Example 1

We demonstrate the standard calls for adaptive numerical integration:

    >> numeric::quadrature(exp(x^2), x=-1..1)

                                    2.925303492

    >> numeric::quadrature(max(1/10, cos(PI*x)), x=-2..0.0123)

                                   0.7521024709

    >> numeric::quadrature(exp(-x^2), x=-2..infinity)

                                    1.768308316

The precision goal is set by `DIGITS`:

    >> DIGITS := 50:

    >> numeric::quadrature(besselJ(0, x), x=0..PI)

                1.3475263146739901712314731279612149642205400522774

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)

                       2.604430665, 2.604430665, 2.604430665

#### 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];
          during evaluation of '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)

                                   0.9460830704

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

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

                                   0.9460830704

Nevertheless, problems may still arise for improper integrals:

    >> numeric::quadrature(ln(1 - cos(x)), x=0..PI)

          Error: singularity [ln]

In this example the integrand is evaluated close to 0. Due to numerical cancellation 1-cos(x) may yield `0.0` for non-zero x. A numerically stable representation is 1-cos(x)=2*sin(x/2)^2:

    >> numeric::quadrature(ln(2*sin(x/2)^2), x=0..PI)

                                    -2.17758609

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

    >> delete f:

#### Example 3

We demonstrate multi-dimensional quadrature:

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

                                   0.6589510757

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)):

for the following quadrature:

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

                                    1.261578952

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 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)

                                       0.665

    >> 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: Can't evaluate to boolean [_less];
          during evaluation of 'f'

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

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

                                   0.5354101317

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)

                                   0.5354101318

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)

                                   0.0238095238

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)

                                       0.53

    >> 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

          3.97732196, 3.977462635, 3.977463259, 3.977463261, 3.977463261

    >> delete a:

#### Example 7

The improper integral int(x^(-9/10), x=0..1)=10 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 call\
          s!
          Increase MaxCalls and try again for a more accurate result. [n\
          umeric::quadrature]

                                    9.998221196

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)

                                    9.999999993

#### Example 8

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

    >> numeric::quadrature(1/x, x=-1..1)

          Warning: Precision goal not achieved after 10000 function call\
          s!
          Increase MaxCalls and try again for a more accurate result. [n\
          umeric::quadrature]

                                    93.14572971

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..1, MaxCalls=infinity)

          Warning: Precision goal not achieved after MAXDEPTH=500 recurs\
          ive calls!
          There may be a singularity of 1/x close to x=1.466369455e-149.
          Increase MAXDEPTH and try again for a more accurate result! [a\
          daptiveQuad]

                                    343.1078544

#### Changes

- Handling of infinite ranges was added. Quadrature now stops after m function calls.
- The option MaxCalls=`m` was introduced to set the maximal number of calls.

Do you have questions or comments?
Copyright © SciFace Software GmbH & Co. KG 2000