Previous Page Next Page Contents

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

Option: method = n

Option: Adaptive= v

Option: MaxCalls= m

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




Do you have questions or comments?


Copyright © SciFace Software GmbH & Co. KG 2000