Tcl Library Source Code

Documentation
Login
EuroTcl/OpenACS 11 - 12 JULY 2024, VIENNA


[ Main Table Of Contents | Table Of Contents | Keyword Index | Categories | Modules | Applications ]

NAME

math::calculus::romberg - Romberg integration

Table Of Contents

SYNOPSIS

package require Tcl 8.5 9
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:

PARAMETERS

OPTIONS

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.

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. 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 © 2004 Kevin B. Kenny . All rights reserved. Redistribution permitted under the terms of the Open Publication License http://www\.opencontent\.org/openpub/