Provided by: tcllib_1.21+dfsg-1_all 
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.
When proposing code changes, please provide unified diffs, i.e the output of diff -u.
Note further that attachments are strongly preferred over inlined patches. Attachments can
be made by going to the Edit form of the ticket immediately after its creation, and then
using the left-most button in the secondary navigation bar.
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/>