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.
numeric::quadrature(f(x
), x = a .. b
, <GaussLegendre = n  GaussTschebyscheff = n  NewtonCotes = n
>, <Adaptive = v
>, <MaxCalls = m
>)
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 userdefined 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.
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
.
Multidimensional 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.
The function is sensitive to the environment variable DIGITS
,
which determines the numerical working precision.
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)
The user should make sure that the integrand is well defined and sufficiently regular. The following fails, because NewtonCotes 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 (GaussLegendre 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 + x^{4})^{2}  1 suffers from severe numerical cancellation and is dominated by roundoff. A numerically stable representation is (1 + x^{4})^{2}  1 = x^{4} (x^{4} + 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:
We demonstrate multidimensional 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^2t*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 GaussLegendre
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:
We demonstrate how integrands given by userdefined 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 multidimensional quadrature of userdefined procedures.
delete f:
The following integrand
f := proc(x, y) begin if x<y then xy else (xy) + (xy)^5 end_if end_proc:
can only be called with numerical arguments and must be delayed twice for 2dimensional 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 higherdimensional
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:
The following example uses nonadaptive GaussTschebyscheff 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:
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)
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.910318545e148'. Increase 'MAXDEPTH' and try again for a more accurate result. [adaptiveQuad]

An arithmetical expression in x 

An identifier or an indexed identifier 

Real or complex numerical values or 

Options, specified as 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 With For DIGITS ≤ 200,
the weights and abscissae of Gaussian quadrature with n For DIGITS > 200,
the routine For With NoteWith Following alternative names for the methods are accepted:


Option, specified as
The default setting is The call
Usually there is no need to switch off the internal adaptive
quadrature via 

Option, specified as
When called interactively, The default value is The default value of 
Floating point number is returned, unless nonnumerical 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.