numeric::quadrature

Numerical integration ( 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

We demonstrate multi-dimensional quadrature:

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

for the following quadrature:

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)

An arithmetical expression in 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.