1math::calculus::romberg(n) Tcl Math Library math::calculus::romberg(n)
2______________________________________________________________________________
6
8 math::calculus::romberg - Romberg integration
9
11 package require Tcl 8.2
12 package require math::calculus 0.6
14 ::math::calculus::romberg f a b ?-option value...?
16 ::math::calculus::romberg_infinity f a b ?-option value...?
18 ::math::calculus::romberg_sqrtSingLower f a b ?-option value...?
20 ::math::calculus::romberg_sqrtSingUpper f a b ?-option value...?
22 ::math::calculus::romberg_powerLawLower gamma f a b ?-option value...?
24 ::math::calculus::romberg_powerLawUpper gamma f a b ?-option value...?
26 ::math::calculus::romberg_expLower f a b ?-option value...?
28 ::math::calculus::romberg_expUpper f a b ?-option value...?
30_________________________________________________________________
32
34 The romberg procedures in the math::calculus package perform numerical
35 integration of a function of one variable. They are intended to be of
36 "production quality" in that they are robust, precise, and reasonably
37 efficient in terms of the number of function evaluations.
38
40 The following procedures are available for Romberg integration:
41 ::math::calculus::romberg f a b ?-option value...?
43 Integrates an analytic function over a given interval.
44 ::math::calculus::romberg_infinity f a b ?-option value...?
46 Integrates an analytic function over a half-infinite interval.
47 ::math::calculus::romberg_sqrtSingLower f a b ?-option value...?
49 Integrates a function that is expected to be analytic over an
50 interval except for the presence of an inverse square root sin‐
51 gularity at the lower limit.
52 ::math::calculus::romberg_sqrtSingUpper f a b ?-option value...?
54 Integrates a function that is expected to be analytic over an
55 interval except for the presence of an inverse square root sin‐
56 gularity at the upper limit.
57 ::math::calculus::romberg_powerLawLower gamma f a b ?-option value...?
59 Integrates a function that is expected to be analytic over an
60 interval except for the presence of a power law singularity at
61 the lower limit.
62 ::math::calculus::romberg_powerLawUpper gamma f a b ?-option value...?
64 Integrates a function that is expected to be analytic over an
65 interval except for the presence of a power law singularity at
66 the upper limit.
67 ::math::calculus::romberg_expLower f a b ?-option value...?
69 Integrates an exponentially growing function; the lower limit of
70 the region of integration may be arbitrarily large and negative.
71 ::math::calculus::romberg_expUpper f a b ?-option value...?
73 Integrates an exponentially decaying function; the upper limit
74 of the region of integration may be arbitrarily large.
75
77 f Function to integrate. Must be expressed as a single Tcl com‐
78 mand, to which will be appended a single argument, specifically,
79 the abscissa at which the function is to be evaluated. The first
80 word of the command will be processed with namespace which in
81 the caller's scope prior to any evaluation. Given this process‐
82 ing, the command may local to the calling namespace rather than
83 needing to be global.
84 a Lower limit of the region of integration.
86 b Upper limit of the region of integration. For the
88 romberg_sqrtSingLower, romberg_sqrtSingUpper, romberg_power‐
89 LawLower, romberg_powerLawUpper, romberg_expLower, and
90 romberg_expUpper procedures, the lower limit must be strictly
91 less than the upper. For the other procedures, the limits may
92 appear in either order.
93 gamma Power to use for a power law singularity; see section IMPROPER
95 INTEGRALS for details.
96
98 -abserror epsilon
99 Requests that the integration machinery proceed at most until
100 the estimated absolute error of the integral is less than
101 epsilon. The error may be seriously over- or underestimated if
102 the function (or any of its derivatives) contains singularities;
103 see section IMPROPER INTEGRALS for details. Default is 1.0e-08.
104 -relerror epsilon
106 Requests that the integration machinery proceed at most until
107 the estimated relative error of the integral is less than
108 epsilon. The error may be seriously over- or underestimated if
109 the function (or any of its derivatives) contains singularities;
110 see section IMPROPER INTEGRALS for details. Default is 1.0e-06.
111 -maxiter m
113 Requests that integration terminate after at most n triplings of
114 the number of evaluations performed. In other words, given n
115 for -maxiter, the integration machinery will make at most 3**n
116 evaluations of the function. Default is 14, corresponding to a
117 limit approximately 4.8 million evaluations. (Well-behaved func‐
118 tions will seldom require more than a few hundred evaluations.)
119 -degree d
121 Requests that an extrapolating polynomial of degree d be used in
122 Romberg integration; see section DESCRIPTION for details.
123 Default is 4. Can be at most m-1.
124
126 The romberg procedure performs Romberg integration using the modified
127 midpoint rule. Romberg integration is an iterative process. At the
128 first step, the function is evaluated at the midpoint of the region of
129 integration, and the value is multiplied by the width of the interval
130 for the coarsest possible estimate. At the second step, the interval
131 is divided into three parts, and the function is evaluated at the mid‐
132 point of each part; the sum of the values is multiplied by three. At
133 the third step, nine parts are used, at the fourth twenty-seven, and so
134 on, tripling the number of subdivisions at each step.
135 Once the interval has been divided at least d times, a polynomial is
137 fitted to the integrals estimated in the last d+1 divisions. The inte‐
138 grals are considered to be a function of the square of the width of the
139 subintervals (any good numerical analysis text will discuss this
140 process under "Romberg integration"). The polynomial is extrapolated
141 to a step size of zero, computing a value for the integral and an esti‐
142 mate of the error.
143 This process will be well-behaved only if the function is analytic over
145 the region of integration; there may be removable singularities at
146 either end of the region provided that the limit of the function (and
147 of all its derivatives) exists as the ends are approached. Thus,
148 romberg may be used to integrate a function like f(x)=sin(x)/x over an
149 interval beginning or ending at zero.
150 Note that romberg will either fail to converge or else return incorrect
152 error estimates if the function, or any of its derivatives, has a sin‐
153 gularity anywhere in the region of integration (except for the case
154 mentioned above). Care must be used, therefore, in integrating a func‐
155 tion like 1/(1-x**2) to avoid the places where the derivative is singu‐
156 lar.
157
159 Romberg integration is also useful for integrating functions over half-
160 infinite intervals or functions that have singularities. The trick is
161 to make a change of variable to eliminate the singularity, and to put
162 the singularity at one end or the other of the region of integration.
163 The math::calculus package supplies a number of romberg procedures to
164 deal with the commoner cases.
165 romberg_infinity
167 Integrates a function over a half-infinite interval; either a or
168 b may be infinite. a and b must be of the same sign; if you
169 need to integrate across the axis, say, from a negative value to
170 positive infinity, use romberg to integrate from the negative
171 value to a small positive value, and then romberg_infinity to
172 integrate from the positive value to positive infinity. The
173 romberg_infinity procedure works by making the change of vari‐
174 able u=1/x, so that the integral from a to b of f(x) is evalu‐
175 ated as the integral from 1/a to 1/b of f(1/u)/u**2.
176 romberg_powerLawLower and romberg_powerLawUpper
178 Integrate a function that has an integrable power law singular‐
179 ity at either the lower or upper bound of the region of integra‐
180 tion (or has a derivative with a power law singularity there).
181 These procedures take a first parameter, gamma, which gives the
182 power law. The function or its first derivative are presumed to
183 diverge as (x-a)**(-gamma) or (b-x)**(-gamma). gamma must be
184 greater than zero and less than 1.
185 These procedures are useful not only in integrating functions
187 that go to infinity at one end of the region of integration, but
188 also functions whose derivatives do not exist at the end of the
189 region. For instance, integrating f(x)=pow(x,0.25) with the
190 origin as one end of the region will result in the romberg pro‐
191 cedure greatly underestimating the error in the integral. The
192 problem can be fixed by observing that the first derivative of
193 f(x), f'(x)=x**(-3/4)/4, goes to infinity at the origin. Inte‐
194 grating using romberg_powerLawLower with gamma set to 0.75 gives
195 much more orderly convergence.
196 These procedures operate by making the change of variable u=(x-
198 a)**(1-gamma) (romberg_powerLawLower) or u=(b-x)**(1-gamma)
199 (romberg_powerLawUpper).
200 To summarize the meaning of gamma:
202 · If f(x) ~ x**(-a) (0 < a < 1), use gamma = a
204 · If f'(x) ~ x**(-b) (0 < b < 1), use gamma = b
206 romberg_sqrtSingLower and romberg_sqrtSingUpper
208 These procedures behave identically to romberg_powerLawLower and
209 romberg_powerLawUpper for the common case of gamma=0.5; that is,
210 they integrate a function with an inverse square root singular‐
211 ity at one end of the interval. They have a simpler implementa‐
212 tion involving square roots rather than arbitrary powers.
213 romberg_expLower and romberg_expUpper
215 These procedures are for integrating a function that grows or
216 decreases exponentially over a half-infinite interval.
217 romberg_expLower handles exponentially growing functions, and
218 allows the lower limit of integration to be an arbitrarily large
219 negative number. romberg_expUpper handles exponentially decay‐
220 ing functions and allows the upper limit of integration to be an
221 arbitrary large positive number. The functions make the change
222 of variable u=exp(-x) and u=exp(x) respectively.
223
225 If you need an improper integral other than the ones listed here, a
226 change of variable can be written in very few lines of Tcl. Because
227 the Tcl coding that does it is somewhat arcane, we offer a worked exam‐
228 ple here.
229 Let's say that the function that we want to integrate is
231 f(x)=exp(x)/sqrt(1-x*x) (not a very natural function, but a good exam‐
232 ple), and we want to integrate it over the interval (-1,1). The denom‐
233 inator falls to zero at both ends of the interval. We wish to make a
234 change of variable from x to u so that dx/sqrt(1-x**2) maps to du.
235 Choosing x=sin(u), we can find that dx=cos(u)*du, and
236 sqrt(1-x**2)=cos(u). The integral from a to b of f(x) is the integral
237 from asin(a) to asin(b) of f(sin(u))*cos(u).
238 We can make a function g that accepts an arbitrary function f and the
240 parameter u, and computes this new integrand.
241 proc g { f u } {
243 set x [expr { sin($u) }]
244 set cmd $f; lappend cmd $x; set y [eval $cmd]
245 return [expr { $y / cos($u) }]
246 }
247 Now integrating f from a to b is the same as integrating g from asin(a)
249 to asin(b). It's a little tricky to get f consistently evaluated in
250 the caller's scope; the following procedure does it.
251 proc romberg_sine { f a b args } {
253 set f [lreplace $f 0 0 [uplevel 1 [list namespace which [lindex $f 0]]]]
254 set f [list g $f]
255 return [eval [linsert $args 0 romberg $f [expr { asin($a) }] [expr { asin($b) }]]]
256 }
257 This romberg_sine procedure will do any function with sqrt(1-x*x) in
259 the denominator. Our sample function is f(x)=exp(x)/sqrt(1-x*x):
260 proc f { x } {
262 expr { exp($x) / sqrt( 1. - $x*$x ) }
263 }
264 Integrating it is a matter of applying romberg_sine as we would any of
266 the other romberg procedures:
267 foreach { value error } [romberg_sine f -1.0 1.0] break
269 puts [format "integral is %.6g +/- %.6g" $value $error]
270 integral is 3.97746 +/- 2.3557e-010
272
275 This document, and the package it describes, will undoubtedly contain
276 bugs and other problems. Please report such in the category math ::
277 calculus of the Tcllib SF Trackers [http://source‐
278 forge.net/tracker/?group_id=12883]. Please also report any ideas for
279 enhancements you may have for either package and/or documentation.
280
282 math::calculus, math::interpolate
283
285 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/>
286math 0.6 math::calculus::romberg(n)