1math::calculus::romberg(n)     Tcl Math Library     math::calculus::romberg(n)
2
3
4
5______________________________________________________________________________
6

NAME

8       math::calculus::romberg - Romberg integration
9

SYNOPSIS

11       package require Tcl  8.2
12
13       package require math::calculus  0.6
14
15       ::math::calculus::romberg f a b ?-option value...?
16
17       ::math::calculus::romberg_infinity f a b ?-option value...?
18
19       ::math::calculus::romberg_sqrtSingLower f a b ?-option value...?
20
21       ::math::calculus::romberg_sqrtSingUpper f a b ?-option value...?
22
23       ::math::calculus::romberg_powerLawLower gamma f a b ?-option value...?
24
25       ::math::calculus::romberg_powerLawUpper gamma f a b ?-option value...?
26
27       ::math::calculus::romberg_expLower f a b ?-option value...?
28
29       ::math::calculus::romberg_expUpper f a b ?-option value...?
30
31______________________________________________________________________________
32

DESCRIPTION

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

PROCEDURES

40       The following procedures are available for Romberg integration:
41
42       ::math::calculus::romberg f a b ?-option value...?
43              Integrates an analytic function over a given interval.
44
45       ::math::calculus::romberg_infinity f a b ?-option value...?
46              Integrates an analytic function over a half-infinite interval.
47
48       ::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
53       ::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
58       ::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
63       ::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
68       ::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
72       ::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

PARAMETERS

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
85       a      Lower limit of the region of integration.
86
87       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
94       gamma  Power  to  use for a power law singularity; see section IMPROPER
95              INTEGRALS for details.
96

OPTIONS

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
105       -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
112       -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
120       -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

DESCRIPTION

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
136       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
144       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
151       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

IMPROPER INTEGRALS

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
166       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
177       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
186              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
197              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
201              To summarize the meaning of gamma:
202
203              ·      If f(x) ~ x**(-a) (0 < a < 1), use gamma = a
204
205              ·      If f'(x) ~ x**(-b) (0 < b < 1), use gamma = b
206
207       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
214       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

OTHER CHANGES OF VARIABLE

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
230       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
239       We can make a function g that accepts an arbitrary function f  and  the
240       parameter u, and computes this new integrand.
241
242
243              proc g { f u } {
244                  set x [expr { sin($u) }]
245                  set cmd $f; lappend cmd $x; set y [eval $cmd]
246                  return [expr { $y / cos($u) }]
247              }
248
249       Now integrating f from a to b is the same as integrating g from asin(a)
250       to asin(b).  It's a little tricky to get f  consistently  evaluated  in
251       the caller's scope; the following procedure does it.
252
253
254              proc romberg_sine { f a b args } {
255                  set f [lreplace $f 0 0 [uplevel 1 [list namespace which [lindex $f 0]]]]
256                  set f [list g $f]
257                  return [eval [linsert $args 0 romberg $f [expr { asin($a) }] [expr { asin($b) }]]]
258              }
259
260       This  romberg_sine  procedure  will do any function with sqrt(1-x*x) in
261       the denominator. Our sample function is f(x)=exp(x)/sqrt(1-x*x):
262
263
264              proc f { x } {
265                  expr { exp($x) / sqrt( 1. - $x*$x ) }
266              }
267
268       Integrating it is a matter of applying romberg_sine as we would any  of
269       the other romberg procedures:
270
271
272              foreach { value error } [romberg_sine f -1.0 1.0] break
273              puts [format "integral is %.6g +/- %.6g" $value $error]
274
275              integral is 3.97746 +/- 2.3557e-010
276
277

BUGS, IDEAS, FEEDBACK

279       This  document,  and the package it describes, will undoubtedly contain
280       bugs and other problems.  Please report such in the  category  math  ::
281       calculus of the Tcllib Trackers [http://core.tcl.tk/tcllib/reportlist].
282       Please also report any ideas for enhancements you may have  for  either
283       package and/or documentation.
284
285       When proposing code changes, please provide unified diffs, i.e the out‐
286       put of diff -u.
287
288       Note further that  attachments  are  strongly  preferred  over  inlined
289       patches.  Attachments  can  be  made  by  going to the Edit form of the
290       ticket immediately after its creation, and  then  using  the  left-most
291       button in the secondary navigation bar.
292

SEE ALSO

294       math::calculus, math::interpolate
295

CATEGORY

297       Mathematics
298
300       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/>
301
302
303
304
305tcllib                                0.6           math::calculus::romberg(n)
Impressum