|
ActiveTcl User Guide
|
|
|
[ Main table Of Contents | Tcllib Table Of Contents | Tcllib Index ]
math::calculus::romberg(n) 0.6 "Math"
math::calculus::romberg - Romberg integration
TABLE OF
CONTENTS
SYNOPSIS
DESCRIPTION
PROCEDURES
PARAMETERS
OPTIONS
DESCRIPTION
IMPROPER
INTEGRALS
OTHER
CHANGES OF VARIABLE
SEE ALSO
COPYRIGHT
package require Tcl 8.2
package require math::calculus 0.6
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.
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.
- 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.
- -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.
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.
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.
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
|
math::calculus , math::interpolate
Copyright © 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/>