numeric::quadrature
--
numerical integrationnumeric::quadrature
(f(x), x=a..b, ..)
computes a numerical approximation of int(f(x), x=a..b).
numeric::quadrature(f(x), x=a..b, <, method=n> <, Adaptive=v> <, MaxCalls=m> )
f(x) |
- | expression in x |
x |
- | identifier or indexed identifier |
a,b |
- | real or complex numerical values or
+-infinity |
|
- | 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. |
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.
The function is sensitive to the environment variable DIGITS
, which determines the
numerical working precision.
int
, numeric::gldata
, numeric::gtdata
, numeric::int
, numeric::ncdata
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.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.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.
quadrature(f(x), x=a..b) = - quadrature(f(x), x=b..a).
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
.
numeric::quadrature
(
numeric::quadrature
(f(x,y), y=A(x)..B(x)), x=a..b)
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.numeric::quadrature
is purely numerical:
no symbolic check for singularities etc. is carried out.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.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.
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).
GaussLegendre, Gauss, GL,
GaussTschebyscheff, GT, GaussChebyshev, GC,
NewtonCotes, NC.
Adaptive=TRUE
. An adaptive
quadrature scheme with automatic control of the quadrature error is
used.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
.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.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
.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.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
.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
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:
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:
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:
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:
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:
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
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
m
was
introduced to set the maximal number of calls.