1[manpage_begin math::calculus::romberg n 0.6] 2[copyright "2004 Kevin B. Kenny <kennykb@acm.org>. All rights\ 3reserved. Redistribution permitted under the terms of the Open\ 4Publication License <http://www.opencontent.org/openpub/>"] 5[moddesc {Tcl Math Library}] 6[titledesc {Romberg integration}] 7[category Mathematics] 8[require Tcl 8.2] 9[require math::calculus 0.6] 10[description] 11[para] 12The [cmd romberg] procedures in the [cmd math::calculus] package 13perform numerical integration of a function of one variable. They 14are intended to be of "production quality" in that they are robust, 15precise, and reasonably efficient in terms of the number of function 16evaluations. 17[section "PROCEDURES"] 18 19The following procedures are available for Romberg integration: 20 21[list_begin definitions] 22[call [cmd ::math::calculus::romberg] [arg f] [arg a] [arg b] [opt [arg "-option value"]...]] 23 24Integrates an analytic function over a given interval. 25 26[call [cmd ::math::calculus::romberg_infinity] [arg f] [arg a] [arg b] [opt [arg "-option value"]...]] 27 28Integrates an analytic function over a half-infinite interval. 29 30[call [cmd ::math::calculus::romberg_sqrtSingLower] [arg f] [arg a] [arg b] [opt [arg "-option value"]...]] 31 32Integrates a function that is expected to be analytic over an interval 33except for the presence of an inverse square root singularity at the 34lower limit. 35 36[call [cmd ::math::calculus::romberg_sqrtSingUpper] [arg f] [arg a] [arg b] [opt [arg "-option value"]...]] 37 38Integrates a function that is expected to be analytic over an interval 39except for the presence of an inverse square root singularity at the 40upper limit. 41 42[call [cmd ::math::calculus::romberg_powerLawLower] [arg gamma] [arg f] [arg a] [arg b] [opt [arg "-option value"]...]] 43 44Integrates a function that is expected to be analytic over an interval 45except for the presence of a power law singularity at the 46lower limit. 47 48[call [cmd ::math::calculus::romberg_powerLawUpper] [arg gamma] [arg f] [arg a] [arg b] [opt [arg "-option value"]...]] 49 50Integrates a function that is expected to be analytic over an interval 51except for the presence of a power law singularity at the 52upper limit. 53 54[call [cmd ::math::calculus::romberg_expLower] [arg f] [arg a] [arg b] [opt [arg "-option value"]...]] 55 56Integrates an exponentially growing function; the lower limit of the 57region of integration may be arbitrarily large and negative. 58 59[call [cmd ::math::calculus::romberg_expUpper] [arg f] [arg a] [arg b] [opt [arg "-option value"]...]] 60 61Integrates an exponentially decaying function; the upper limit of the 62region of integration may be arbitrarily large. 63 64[list_end] 65 66[section PARAMETERS] 67 68[list_begin definitions] 69 70[def [arg f]] 71 72Function to integrate. Must be expressed as a single Tcl command, 73to which will be appended a single argument, specifically, the 74abscissa at which the function is to be evaluated. The first word 75of the command will be processed with [cmd "namespace which"] in the 76caller's scope prior to any evaluation. Given this processing, the 77command may local to the calling namespace rather than needing to be 78global. 79 80[def [arg a]] 81 82Lower limit of the region of integration. 83 84[def [arg b]] 85 86Upper limit of the region of integration. For the 87[cmd romberg_sqrtSingLower], [cmd romberg_sqrtSingUpper], 88[cmd romberg_powerLawLower], [cmd romberg_powerLawUpper], 89[cmd romberg_expLower], and [cmd romberg_expUpper] procedures, 90the lower limit must be strictly less than the upper. For 91the other procedures, the limits may appear in either order. 92 93[def [arg gamma]] 94 95Power to use for a power law singularity; see section 96[sectref "IMPROPER INTEGRALS"] for details. 97 98[list_end] 99 100[section OPTIONS] 101 102[list_begin definitions] 103 104[def "[option -abserror] [arg epsilon]"] 105 106Requests that the integration machinery proceed at most until 107the estimated absolute error of the integral is less than 108[arg epsilon]. The error may be seriously over- or underestimated 109if the function (or any of its derivatives) contains singularities; 110see section [sectref "IMPROPER INTEGRALS"] for details. Default 111is 1.0e-08. 112 113[def "[option -relerror] [arg epsilon]"] 114 115Requests that the integration machinery proceed at most until 116the estimated relative error of the integral is less than 117[arg epsilon]. The error may be seriously over- or underestimated 118if the function (or any of its derivatives) contains singularities; 119see section [sectref "IMPROPER INTEGRALS"] for details. Default is 1201.0e-06. 121 122[def "[option -maxiter] [arg m]"] 123 124Requests that integration terminate after at most [arg n] triplings of 125the number of evaluations performed. In other words, given [arg n] 126for [option -maxiter], the integration machinery will make at most 1273**[arg n] evaluations of the function. Default is 14, corresponding 128to a limit approximately 4.8 million evaluations. (Well-behaved 129functions will seldom require more than a few hundred evaluations.) 130 131[def "[option -degree] [arg d]"] 132 133Requests that an extrapolating polynomial of degree [arg d] be used 134in Romberg integration; see section [sectref "DESCRIPTION"] for 135details. Default is 4. Can be at most [arg m]-1. 136 137[list_end] 138 139[section DESCRIPTION] 140 141The [cmd romberg] procedure performs Romberg integration using 142the modified midpoint rule. Romberg integration is an iterative 143process. At the first step, the function is evaluated at the 144midpoint of the region of integration, and the value is multiplied by 145the width of the interval for the coarsest possible estimate. 146At the second step, the interval is divided into three parts, 147and the function is evaluated at the midpoint of each part; the 148sum of the values is multiplied by three. At the third step, 149nine parts are used, at the fourth twenty-seven, and so on, 150tripling the number of subdivisions at each step. 151 152[para] 153 154Once the interval has been divided at least [arg d] times, 155a polynomial is fitted to the integrals estimated in the last 156[arg d]+1 divisions. The integrals are considered to be a 157function of the square of the width of the subintervals 158(any good numerical analysis text will discuss this process 159under "Romberg integration"). The polynomial is extrapolated 160to a step size of zero, computing a value for the integral and 161an estimate of the error. 162 163[para] 164 165This process will be well-behaved only if the function is analytic 166over the region of integration; there may be removable singularities 167at either end of the region provided that the limit of the function 168(and of all its derivatives) exists as the ends are approached. 169Thus, [cmd romberg] may be used to integrate a function like 170f(x)=sin(x)/x over an interval beginning or ending at zero. 171 172[para] 173 174Note that [cmd romberg] will either fail to converge or else return 175incorrect error estimates if the function, or any of its derivatives, 176has a singularity anywhere in the region of integration (except for 177the case mentioned above). Care must be used, therefore, in 178integrating a function like 1/(1-x**2) to avoid the places 179where the derivative is singular. 180 181[section "IMPROPER INTEGRALS"] 182 183Romberg integration is also useful for integrating functions over 184half-infinite intervals or functions that have singularities. 185The trick is to make a change of variable to eliminate the 186singularity, and to put the singularity at one end or the other 187of the region of integration. The [cmd math::calculus] package 188supplies a number of [cmd romberg] procedures to deal with the 189commoner cases. 190 191[list_begin definitions] 192 193[def [cmd romberg_infinity]] 194 195Integrates a function over a half-infinite interval; either 196[arg a] or [arg b] may be infinite. [arg a] and [arg b] must be 197of the same sign; if you need to integrate across the axis, 198say, from a negative value to positive infinity, 199use [cmd romberg] to integrate from the negative 200value to a small positive value, and then [cmd romberg_infinity] 201to integrate from the positive value to positive infinity. The 202[cmd romberg_infinity] procedure works by making the change of 203variable u=1/x, so that the integral from a to b of f(x) is 204evaluated as the integral from 1/a to 1/b of f(1/u)/u**2. 205 206[def "[cmd romberg_powerLawLower] and [cmd romberg_powerLawUpper]"] 207 208Integrate a function that has an integrable power law singularity 209at either the lower or upper bound of the region of integration 210(or has a derivative with a power law singularity there). 211These procedures take a first parameter, [arg gamma], which gives 212the power law. The function or its first derivative are presumed to diverge as 213(x-[arg a])**(-[arg gamma]) or ([arg b]-x)**(-[arg gamma]). [arg gamma] 214must be greater than zero and less than 1. 215 216[para] 217 218These procedures are useful not only in integrating functions 219that go to infinity at one end of the region of integration, but 220also functions whose derivatives do not exist at the end of 221the region. For instance, integrating f(x)=pow(x,0.25) with the 222origin as one end of the region will result in the [cmd romberg] 223procedure greatly underestimating the error in the integral. 224The problem can be fixed by observing that the first derivative 225of f(x), f'(x)=x**(-3/4)/4, goes to infinity at the origin. Integrating 226using [cmd romberg_powerLawLower] with [arg gamma] set to 0.75 227gives much more orderly convergence. 228 229[para] 230 231These procedures operate by making the change of variable 232u=(x-a)**(1-gamma) ([cmd romberg_powerLawLower]) or 233u=(b-x)**(1-gamma) ([cmd romberg_powerLawUpper]). 234 235[para] 236 237To summarize the meaning of gamma: 238 239[list_begin itemized] 240[item] 241If f(x) ~ x**(-a) (0 < a < 1), use gamma = a 242[item] 243If f'(x) ~ x**(-b) (0 < b < 1), use gamma = b 244[list_end] 245 246[def "[cmd romberg_sqrtSingLower] and [cmd romberg_sqrtSingUpper]"] 247 248These procedures behave identically to [cmd romberg_powerLawLower] and 249[cmd romberg_powerLawUpper] for the common case of [arg gamma]=0.5; 250that is, they integrate a function with an inverse square root 251singularity at one end of the interval. They have a simpler 252implementation involving square roots rather than arbitrary powers. 253 254[def "[cmd romberg_expLower] and [cmd romberg_expUpper]"] 255 256These procedures are for integrating a function that grows or 257decreases exponentially over a half-infinite interval. 258[cmd romberg_expLower] handles exponentially growing functions, and 259allows the lower limit of integration to be an arbitrarily large 260negative number. [cmd romberg_expUpper] handles exponentially 261decaying functions and allows the upper limit of integration to 262be an arbitrary large positive number. The functions make the 263change of variable u=exp(-x) and u=exp(x) respectively. 264 265[list_end] 266 267[section "OTHER CHANGES OF VARIABLE"] 268 269If you need an improper integral other than the ones listed here, 270a change of variable can be written in very few lines of Tcl. 271Because the Tcl coding that does it is somewhat arcane, 272we offer a worked example here. 273 274[para] 275 276Let's say that the function that we want to integrate is 277f(x)=exp(x)/sqrt(1-x*x) (not a very natural 278function, but a good example), and we want to integrate 279it over the interval (-1,1). The denominator falls to zero 280at both ends of the interval. We wish to make a change of variable 281from x to u 282so that dx/sqrt(1-x**2) maps to du. Choosing x=sin(u), we can 283find that dx=cos(u)*du, and sqrt(1-x**2)=cos(u). The integral 284from a to b of f(x) is the integral from asin(a) to asin(b) 285of f(sin(u))*cos(u). 286 287[para] 288 289We can make a function [cmd g] that accepts an arbitrary function 290[cmd f] and the parameter u, and computes this new integrand. 291 292[example { 293proc g { f u } { 294 set x [expr { sin($u) }] 295 set cmd $f; lappend cmd $x; set y [eval $cmd] 296 return [expr { $y / cos($u) }] 297} 298}] 299 300Now integrating [cmd f] from [arg a] to [arg b] is the same 301as integrating [cmd g] from [arg asin(a)] to [arg asin(b)]. 302It's a little tricky to get [cmd f] consistently evaluated in 303the caller's scope; the following procedure does it. 304 305[example { 306proc romberg_sine { f a b args } { 307 set f [lreplace $f 0 0\ 308 [uplevel 1 [list namespace which [lindex $f 0]]]] 309 set f [list g $f] 310 return [eval [linsert $args 0\ 311 romberg $f\ 312 [expr { asin($a) }] [expr { asin($b) }]]] 313} 314}] 315 316This [cmd romberg_sine] procedure will do any function with 317sqrt(1-x*x) in the denominator. Our sample function is 318f(x)=exp(x)/sqrt(1-x*x): 319 320[example { 321proc f { x } { 322 expr { exp($x) / sqrt( 1. - $x*$x ) } 323} 324}] 325 326Integrating it is a matter of applying [cmd romberg_sine] 327as we would any of the other [cmd romberg] procedures: 328 329[example { 330foreach { value error } [romberg_sine f -1.0 1.0] break 331puts [format "integral is %.6g +/- %.6g" $value $error] 332 333integral is 3.97746 +/- 2.3557e-010 334}] 335 336 337[section {BUGS, IDEAS, FEEDBACK}] 338 339This document, and the package it describes, will undoubtedly contain 340bugs and other problems. 341 342Please report such in the category [emph {math :: calculus}] of the 343[uri {http://sourceforge.net/tracker/?group_id=12883} {Tcllib SF Trackers}]. 344 345Please also report any ideas for enhancements you may have for either 346package and/or documentation. 347 348 349 350[see_also math::calculus math::interpolate] 351[manpage_end] 352