Страницы справочника OS X

Эта страница руководства для  версии 10.9 Mac OS X

Читать страницы руководства

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



____________________________________________________________________________________________________________

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_sqrtSingUp-per, romberg_sqrtSingUpper,
              per, romberg_powerLawLower, romberg_powerLawUpper, romberg_expLower, and romberg_expUpper pro-cedures, procedures,
              cedures,  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 eval-uations.) evaluations.)
              uations.)

       -degree d
              Requests that an extrapolating polynomial of degree d be used in Romberg integration; see sec-tion section
              tion 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 integra-tion integration
       tion 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 pos-sible possible
       sible 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 esti-mated estimated
       mated 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 func-tion function
       tion  (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  func-tion, function,
       tion,  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 func-tions functions
       tions that have singularities.  The trick is to make a change of variable to eliminate the  singular-ity, singularity,
       ity, and to put the singularity at one end or the other of the region of integration.  The math::cal-culus math::calculus
       culus 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 inte-gral integral
              gral 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_power-LawLower) (romberg_powerLawLower)
              LawLower) 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 sin-
              gularity 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 SF Trackers [http://source -
       forge.net/tracker/? group_id=12883].  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/>




math                                                 0.6                          math::calculus::romberg(n)