math::calculus::romberg(n) Tcl Math Library math::calculus::romberg(n) 







math::calculus::romberg - Romberg integration








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...?












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.











  

  
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.

  

  
Lower limit of the region of integration.

  

  
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.

  

  
Power to use for a power law singularity; see section IMPROPER
      INTEGRALS for details.











  

  
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.

  

  
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.

  

  
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.)

  

  
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.



  

  
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.

  

  
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
    • 






  

  
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.

  

  
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









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 SF Trackers
    [http://sourceforge.net/tracker/?group_id=12883]. Please also report any
    ideas for enhancements you may have for either package and/or
  documentation.








math::calculus, math::interpolate








Mathematics








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/>








  

    0.6
    math