NAME¶
math::calculus::romberg - Romberg integration
SYNOPSIS¶
package require
Tcl 8.2
package require
math::calculus 0.6
::math::calculus::romberg f a b ?
-option
value...?
::math::calculus::romberg_infinity f a b ?
-option
value...?
::math::calculus::romberg_sqrtSingLower f a b
?
-option value...?
::math::calculus::romberg_sqrtSingUpper f a b
?
-option value...?
::math::calculus::romberg_powerLawLower gamma f a
b ?
-option value...?
::math::calculus::romberg_powerLawUpper gamma f a
b ?
-option value...?
::math::calculus::romberg_expLower f a b ?
-option
value...?
::math::calculus::romberg_expUpper f a b ?
-option
value...?
DESCRIPTION¶
The
romberg procedures in the
math::calculus package perform
numerical integration of a function of one variable. They are intended to be
of "production quality" in that they are robust, precise, and
reasonably efficient in terms of the number of function evaluations.
PROCEDURES¶
The following procedures are available for Romberg integration:
- ::math::calculus::romberg f a b ?-option
value...?
- Integrates an analytic function over a given interval.
- ::math::calculus::romberg_infinity f a b
?-option value...?
- Integrates an analytic function over a half-infinite interval.
- ::math::calculus::romberg_sqrtSingLower f a b
? -option value...?
- Integrates a function that is expected to be analytic over an interval
except for the presence of an inverse square root singularity at the lower
limit.
- ::math::calculus::romberg_sqrtSingUpper f a b
? -option value...?
- Integrates a function that is expected to be analytic over an interval
except for the presence of an inverse square root singularity at the upper
limit.
- ::math::calculus::romberg_powerLawLower gamma f
a b ?-option value...?
- Integrates a function that is expected to be analytic over an interval
except for the presence of a power law singularity at the lower
limit.
- ::math::calculus::romberg_powerLawUpper gamma f
a b ?-option value...?
- Integrates a function that is expected to be analytic over an interval
except for the presence of a power law singularity at the upper
limit.
- ::math::calculus::romberg_expLower f a b
?-option value...?
- Integrates an exponentially growing function; the lower limit of the
region of integration may be arbitrarily large and negative.
- ::math::calculus::romberg_expUpper f a b
?-option value...?
- Integrates an exponentially decaying function; the upper limit of the
region of integration may be arbitrarily large.
PARAMETERS¶
- f
- Function to integrate. Must be expressed as a single Tcl command, to which
will be appended a single argument, specifically, the abscissa at which
the function is to be evaluated. The first word of the command will be
processed with namespace which in the caller's scope prior to any
evaluation. Given this processing, the command may local to the calling
namespace rather than needing to be global.
- a
- Lower limit of the region of integration.
- b
- Upper limit of the region of integration. For the
romberg_sqrtSingLower, romberg_sqrtSingUpper,
romberg_powerLawLower, romberg_powerLawUpper,
romberg_expLower, and romberg_expUpper procedures, the lower
limit must be strictly less than the upper. For the other procedures, the
limits may appear in either order.
- gamma
- Power to use for a power law singularity; see section IMPROPER
INTEGRALS for details.
OPTIONS¶
- -abserror epsilon
- Requests that the integration machinery proceed at most until the
estimated absolute error of the integral is less than epsilon. The
error may be seriously over- or underestimated if the function (or any of
its derivatives) contains singularities; see section IMPROPER
INTEGRALS for details. Default is 1.0e-08.
- -relerror epsilon
- Requests that the integration machinery proceed at most until the
estimated relative error of the integral is less than epsilon. The
error may be seriously over- or underestimated if the function (or any of
its derivatives) contains singularities; see section IMPROPER
INTEGRALS for details. Default is 1.0e-06.
- -maxiter m
- Requests that integration terminate after at most n triplings of
the number of evaluations performed. In other words, given n for
-maxiter, the integration machinery will make at most 3** n
evaluations of the function. Default is 14, corresponding to a limit
approximately 4.8 million evaluations. (Well-behaved functions will seldom
require more than a few hundred evaluations.)
- -degree d
- Requests that an extrapolating polynomial of degree d be used in
Romberg integration; see section DESCRIPTION for details. Default
is 4. Can be at most m-1.
DESCRIPTION¶
The
romberg procedure performs Romberg integration using the modified
midpoint rule. Romberg integration is an iterative process. At the first step,
the function is evaluated at the midpoint of the region of integration, and
the value is multiplied by the width of the interval for the coarsest possible
estimate. At the second step, the interval is divided into three parts, and
the function is evaluated at the midpoint of each part; the sum of the values
is multiplied by three. At the third step, nine parts are used, at the fourth
twenty-seven, and so on, tripling the number of subdivisions at each step.
Once the interval has been divided at least
d times, a polynomial is
fitted to the integrals estimated in the last
d+1 divisions. The
integrals are considered to be a function of the square of the width of the
subintervals (any good numerical analysis text will discuss this process under
"Romberg integration"). The polynomial is extrapolated to a step
size of zero, computing a value for the integral and an estimate of the error.
This process will be well-behaved only if the function is analytic over the
region of integration; there may be removable singularities at either end of
the region provided that the limit of the function (and of all its
derivatives) exists as the ends are approached. Thus,
romberg may be
used to integrate a function like f(x)=sin(x)/x over an interval beginning or
ending at zero.
Note that
romberg will either fail to converge or else return incorrect
error estimates if the function, or any of its derivatives, has a singularity
anywhere in the region of integration (except for the case mentioned above).
Care must be used, therefore, in integrating a function like 1/(1-x**2) to
avoid the places where the derivative is singular.
IMPROPER INTEGRALS¶
Romberg integration is also useful for integrating functions over half-infinite
intervals or functions that have singularities. The trick is to make a change
of variable to eliminate the singularity, and to put the singularity at one
end or the other of the region of integration. The
math::calculus
package supplies a number of
romberg procedures to deal with the
commoner cases.
- romberg_infinity
- Integrates a function over a half-infinite interval; either a or
b may be infinite. a and b must be of the same sign;
if you need to integrate across the axis, say, from a negative value to
positive infinity, use romberg to integrate from the negative value
to a small positive value, and then romberg_infinity to integrate
from the positive value to positive infinity. The romberg_infinity
procedure works by making the change of variable u=1/x, so that the
integral from a to b of f(x) is evaluated as the integral from 1/a to 1/b
of f(1/u)/u**2.
- romberg_powerLawLower and romberg_powerLawUpper
- Integrate a function that has an integrable power law singularity at
either the lower or upper bound of the region of integration (or has a
derivative with a power law singularity there). These procedures take a
first parameter, gamma, which gives the power law. The function or
its first derivative are presumed to diverge as (x-
a)**(-gamma) or ( b-x)**(-gamma). gamma
must be greater than zero and less than 1.
These procedures are useful not only in integrating functions that go to
infinity at one end of the region of integration, but also functions whose
derivatives do not exist at the end of the region. For instance,
integrating f(x)=pow(x,0.25) with the origin as one end of the region will
result in the romberg procedure greatly underestimating the error
in the integral. The problem can be fixed by observing that the first
derivative of f(x), f'(x)=x**(-3/4)/4, goes to infinity at the origin.
Integrating using romberg_powerLawLower with gamma set to
0.75 gives much more orderly convergence.
These procedures operate by making the change of variable u=(x-a)**(1-gamma)
( romberg_powerLawLower) or u=(b-x)**(1-gamma) (
romberg_powerLawUpper).
To summarize the meaning of gamma:
- •
- If f(x) ~ x**(-a) (0 < a < 1), use gamma = a
- •
- If f'(x) ~ x**(-b) (0 < b < 1), use gamma = b
- romberg_sqrtSingLower and romberg_sqrtSingUpper
- These procedures behave identically to romberg_powerLawLower and
romberg_powerLawUpper for the common case of gamma=0.5; that
is, they integrate a function with an inverse square root singularity at
one end of the interval. They have a simpler implementation involving
square roots rather than arbitrary powers.
- romberg_expLower and romberg_expUpper
- These procedures are for integrating a function that grows or decreases
exponentially over a half-infinite interval. romberg_expLower
handles exponentially growing functions, and allows the lower limit of
integration to be an arbitrarily large negative number.
romberg_expUpper handles exponentially decaying functions and
allows the upper limit of integration to be an arbitrary large positive
number. The functions make the change of variable u=exp(-x) and u=exp(x)
respectively.
OTHER CHANGES OF VARIABLE¶
If you need an improper integral other than the ones listed here, a change of
variable can be written in very few lines of Tcl. Because the Tcl coding that
does it is somewhat arcane, we offer a worked example here.
Let's say that the function that we want to integrate is f(x)=exp(x)/sqrt(1-x*x)
(not a very natural function, but a good example), and we want to integrate it
over the interval (-1,1). The denominator falls to zero at both ends of the
interval. We wish to make a change of variable from x to u so that
dx/sqrt(1-x**2) maps to du. Choosing x=sin(u), we can find that dx=cos(u)*du,
and sqrt(1-x**2)=cos(u). The integral from a to b of f(x) is the integral from
asin(a) to asin(b) of f(sin(u))*cos(u).
We can make a function
g that accepts an arbitrary function
f and
the parameter u, and computes this new integrand.
proc g { f u } {
set x [expr { sin($u) }]
set cmd $f; lappend cmd $x; set y [eval $cmd]
return [expr { $y / cos($u) }]
}
Now integrating
f from
a to
b is the same as integrating
g from
asin(a) to
asin(b). It's a little tricky to get
f consistently evaluated in the caller's scope; the following procedure
does it.
proc romberg_sine { f a b args } {
set f [lreplace $f 0 0 [uplevel 1 [list namespace which [lindex $f 0]]]]
set f [list g $f]
return [eval [linsert $args 0 romberg $f [expr { asin($a) }] [expr { asin($b) }]]]
}
This
romberg_sine procedure will do any function with sqrt(1-x*x) in the
denominator. Our sample function is f(x)=exp(x)/sqrt(1-x*x):
proc f { x } {
expr { exp($x) / sqrt( 1. - $x*$x ) }
}
Integrating it is a matter of applying
romberg_sine as we would any of
the other
romberg procedures:
foreach { value error } [romberg_sine f -1.0 1.0] break
puts [format "integral is %.6g +/- %.6g" $value $error]
integral is 3.97746 +/- 2.3557e-010
BUGS, IDEAS, FEEDBACK¶
This document, and the package it describes, will undoubtedly contain bugs and
other problems. Please report such in the category
math :: calculus of
the
Tcllib Trackers [
http://core.tcl.tk/tcllib/reportlist]. Please also
report any ideas for enhancements you may have for either package and/or
documentation.
SEE ALSO¶
math::calculus, math::interpolate
CATEGORY¶
Mathematics
COPYRIGHT¶
Copyright (c) 2004 Kevin B. Kenny <kennykb@acm.org>. All rights reserved. Redistribution permitted under the terms of the Open Publication License <http://www.opencontent.org/openpub/>