source: tol/tolsh/tolsh.vfs/lib/tcllib1.10/math/calculus.tcl @ 8597867d

1941-docker-image-to-testdevel
Last change on this file since 8597867d was 8597867d, checked in by Jorge Pérez <josp.jorge@…>, 16 years ago

closes #976, se incluye tcllib1.10, base64 debe eliminarse del raiz de lib ya que esta incluido en tcllib

git-svn-id: ​https://www.tol-project.org/svn/tolp/trunk@2556 be4846ae-f051-483a-831f-b64e7d08527f
  • Property mode set to 100644
File size: 44.5 KB
Line 
1# calculus.tcl --
2#    Package that implements several basic numerical methods, such
3#    as the integration of a one-dimensional function and the
4#    solution of a system of first-order differential equations.
5#
6# Copyright (c) 2002, 2003, 2004, 2006 by Arjen Markus.
7# Copyright (c) 2004 by Kevin B. Kenny.  All rights reserved.
8# See the file "license.terms" for information on usage and redistribution
9# of this file, and for a DISCLAIMER OF ALL WARRANTIES.
10#
11# RCS: @(#) $Id: calculus.tcl,v 1.13 2006/03/28 20:36:00 arjenmarkus Exp $
12
13package require Tcl 8.4
14package require math::interpolate
15package provide math::calculus 0.7
16
17# math::calculus --
18#    Namespace for the commands
19
20namespace eval ::math::calculus {
21
22    namespace import ::math::interpolate::neville
23
24    namespace import ::math::expectDouble ::math::expectInteger
25
26    namespace export \
27        integral integralExpr integral2D integral3D \
28        eulerStep heunStep rungeKuttaStep           \
29        boundaryValueSecondOrder solveTriDiagonal   \
30        newtonRaphson newtonRaphsonParameters
31    namespace export \
32        integral2D_2accurate integral3D_accurate
33
34    namespace export romberg romberg_infinity
35    namespace export romberg_sqrtSingLower romberg_sqrtSingUpper
36    namespace export romberg_powerLawLower romberg_powerLawUpper
37    namespace export romberg_expLower romberg_expUpper
38
39    namespace export regula_falsi
40
41    variable nr_maxiter    20
42    variable nr_tolerance   0.001
43
44}
45
46# integral --
47#    Integrate a function over a given interval using the Simpson rule
48#
49# Arguments:
50#    begin       Start of the interval
51#    end         End of the interval
52#    nosteps     Number of steps in which to divide the interval
53#    func        Name of the function to be integrated (takes one
54#                argument)
55# Return value:
56#    Computed integral
57#
58proc ::math::calculus::integral { begin end nosteps func } {
59
60   set delta    [expr {($end-$begin)/double($nosteps)}]
61   set hdelta   [expr {$delta/2.0}]
62   set result   0.0
63   set xval     $begin
64   set func_end [uplevel 1 $func $xval]
65   for { set i 1 } { $i <= $nosteps } { incr i } {
66      set func_begin  $func_end
67      set func_middle [uplevel 1 $func [expr {$xval+$hdelta}]]
68      set func_end    [uplevel 1 $func [expr {$xval+$delta}]]
69      set result      [expr  {$result+$func_begin+4.0*$func_middle+$func_end}]
70
71      set xval        [expr {$begin+double($i)*$delta}]
72   }
73
74   return [expr {$result*$delta/6.0}]
75}
76
77# integralExpr --
78#    Integrate an expression with "x" as the integrate according to the
79#    Simpson rule
80#
81# Arguments:
82#    begin       Start of the interval
83#    end         End of the interval
84#    nosteps     Number of steps in which to divide the interval
85#    expression  Expression with "x" as the integrate
86# Return value:
87#    Computed integral
88#
89proc ::math::calculus::integralExpr { begin end nosteps expression } {
90
91   set delta    [expr {($end-$begin)/double($nosteps)}]
92   set hdelta   [expr {$delta/2.0}]
93   set result   0.0
94   set x        $begin
95   # FRINK: nocheck
96   set func_end [expr $expression]
97   for { set i 1 } { $i <= $nosteps } { incr i } {
98      set func_begin  $func_end
99      set x           [expr {$x+$hdelta}]
100       # FRINK: nocheck
101      set func_middle [expr $expression]
102      set x           [expr {$x+$hdelta}]
103       # FRINK: nocheck
104      set func_end    [expr $expression]
105      set result      [expr {$result+$func_begin+4.0*$func_middle+$func_end}]
106
107      set x           [expr {$begin+double($i)*$delta}]
108   }
109
110   return [expr {$result*$delta/6.0}]
111}
112
113# integral2D --
114#    Integrate a given fucntion of two variables over a block,
115#    using bilinear interpolation (for this moment: block function)
116#
117# Arguments:
118#    xinterval   Start, stop and number of steps of the "x" interval
119#    yinterval   Start, stop and number of steps of the "y" interval
120#    func        Function of the two variables to be integrated
121# Return value:
122#    Computed integral
123#
124proc ::math::calculus::integral2D { xinterval yinterval func } {
125
126   foreach { xbegin xend xnumber } $xinterval { break }
127   foreach { ybegin yend ynumber } $yinterval { break }
128
129   set xdelta    [expr {($xend-$xbegin)/double($xnumber)}]
130   set ydelta    [expr {($yend-$ybegin)/double($ynumber)}]
131   set hxdelta   [expr {$xdelta/2.0}]
132   set hydelta   [expr {$ydelta/2.0}]
133   set result   0.0
134   set dxdy      [expr {$xdelta*$ydelta}]
135   for { set j 0 } { $j < $ynumber } { incr j } {
136      set y [expr {$ybegin+$hydelta+double($j)*$ydelta}]
137      for { set i 0 } { $i < $xnumber } { incr i } {
138         set x           [expr {$xbegin+$hxdelta+double($i)*$xdelta}]
139         set func_value  [uplevel 1 $func $x $y]
140         set result      [expr {$result+$func_value}]
141      }
142   }
143
144   return [expr {$result*$dxdy}]
145}
146
147# integral3D --
148#    Integrate a given fucntion of two variables over a block,
149#    using trilinear interpolation (for this moment: block function)
150#
151# Arguments:
152#    xinterval   Start, stop and number of steps of the "x" interval
153#    yinterval   Start, stop and number of steps of the "y" interval
154#    zinterval   Start, stop and number of steps of the "z" interval
155#    func        Function of the three variables to be integrated
156# Return value:
157#    Computed integral
158#
159proc ::math::calculus::integral3D { xinterval yinterval zinterval func } {
160
161   foreach { xbegin xend xnumber } $xinterval { break }
162   foreach { ybegin yend ynumber } $yinterval { break }
163   foreach { zbegin zend znumber } $zinterval { break }
164
165   set xdelta    [expr {($xend-$xbegin)/double($xnumber)}]
166   set ydelta    [expr {($yend-$ybegin)/double($ynumber)}]
167   set zdelta    [expr {($zend-$zbegin)/double($znumber)}]
168   set hxdelta   [expr {$xdelta/2.0}]
169   set hydelta   [expr {$ydelta/2.0}]
170   set hzdelta   [expr {$zdelta/2.0}]
171   set result   0.0
172   set dxdydz    [expr {$xdelta*$ydelta*$zdelta}]
173   for { set k 0 } { $k < $znumber } { incr k } {
174      set z [expr {$zbegin+$hzdelta+double($k)*$zdelta}]
175      for { set j 0 } { $j < $ynumber } { incr j } {
176         set y [expr {$ybegin+$hydelta+double($j)*$ydelta}]
177         for { set i 0 } { $i < $xnumber } { incr i } {
178            set x           [expr {$xbegin+$hxdelta+double($i)*$xdelta}]
179            set func_value  [uplevel 1 $func $x $y $z]
180            set result      [expr {$result+$func_value}]
181         }
182      }
183   }
184
185   return [expr {$result*$dxdydz}]
186}
187
188# integral2D_accurate --
189#    Integrate a given function of two variables over a block,
190#    using a four-point quadrature formula
191#
192# Arguments:
193#    xinterval   Start, stop and number of steps of the "x" interval
194#    yinterval   Start, stop and number of steps of the "y" interval
195#    func        Function of the two variables to be integrated
196# Return value:
197#    Computed integral
198#
199proc ::math::calculus::integral2D_accurate { xinterval yinterval func } {
200
201    foreach { xbegin xend xnumber } $xinterval { break }
202    foreach { ybegin yend ynumber } $yinterval { break }
203
204    set alpha     [expr {sqrt(2.0/3.0)}]
205    set minalpha  [expr {-$alpha}]
206    set dpoints   [list $alpha 0.0 $minalpha 0.0 0.0 $alpha 0.0 $minalpha]
207
208    set xdelta    [expr {($xend-$xbegin)/double($xnumber)}]
209    set ydelta    [expr {($yend-$ybegin)/double($ynumber)}]
210    set hxdelta   [expr {$xdelta/2.0}]
211    set hydelta   [expr {$ydelta/2.0}]
212    set result   0.0
213    set dxdy      [expr {0.25*$xdelta*$ydelta}]
214
215    for { set j 0 } { $j < $ynumber } { incr j } {
216        set y [expr {$ybegin+$hydelta+double($j)*$ydelta}]
217        for { set i 0 } { $i < $xnumber } { incr i } {
218            set x [expr {$xbegin+$hxdelta+double($i)*$xdelta}]
219
220            foreach {dx dy} $dpoints {
221                set x1          [expr {$x+$dx}]
222                set y1          [expr {$y+$dy}]
223                set func_value  [uplevel 1 $func $x1 $y1]
224                set result      [expr {$result+$func_value}]
225            }
226        }
227    }
228
229    return [expr {$result*$dxdy}]
230}
231
232# integral3D_accurate --
233#    Integrate a given function of three variables over a block,
234#    using an 8-point quadrature formula
235#
236# Arguments:
237#    xinterval   Start, stop and number of steps of the "x" interval
238#    yinterval   Start, stop and number of steps of the "y" interval
239#    zinterval   Start, stop and number of steps of the "z" interval
240#    func        Function of the three variables to be integrated
241# Return value:
242#    Computed integral
243#
244proc ::math::calculus::integral3D_accurate { xinterval yinterval zinterval func } {
245
246    foreach { xbegin xend xnumber } $xinterval { break }
247    foreach { ybegin yend ynumber } $yinterval { break }
248    foreach { zbegin zend znumber } $zinterval { break }
249
250    set alpha     [expr {sqrt(1.0/3.0)}]
251    set minalpha  [expr {-$alpha}]
252
253    set dpoints   [list $alpha    $alpha    $alpha    \
254                        $alpha    $alpha    $minalpha \
255                        $alpha    $minalpha $alpha    \
256                        $alpha    $minalpha $minalpha \
257                        $minalpha $alpha    $alpha    \
258                        $minalpha $alpha    $minalpha \
259                        $minalpha $minalpha $alpha    \
260                        $minalpha $minalpha $minalpha ]
261
262    set xdelta    [expr {($xend-$xbegin)/double($xnumber)}]
263    set ydelta    [expr {($yend-$ybegin)/double($ynumber)}]
264    set zdelta    [expr {($zend-$zbegin)/double($znumber)}]
265    set hxdelta   [expr {$xdelta/2.0}]
266    set hydelta   [expr {$ydelta/2.0}]
267    set hzdelta   [expr {$zdelta/2.0}]
268    set result    0.0
269    set dxdydz    [expr {0.125*$xdelta*$ydelta*$zdelta}]
270
271    for { set k 0 } { $k < $znumber } { incr k } {
272        set z [expr {$zbegin+$hzdelta+double($k)*$zdelta}]
273        for { set j 0 } { $j < $ynumber } { incr j } {
274            set y [expr {$ybegin+$hydelta+double($j)*$ydelta}]
275            for { set i 0 } { $i < $xnumber } { incr i } {
276                set x [expr {$xbegin+$hxdelta+double($i)*$xdelta}]
277
278                foreach {dx dy dz} $dpoints {
279                    set x1 [expr {$x+$dx}]
280                    set y1 [expr {$y+$dy}]
281                    set z1 [expr {$z+$dz}]
282                    set func_value  [uplevel 1 $func $x1 $y1 $z1]
283                    set result      [expr {$result+$func_value}]
284                }
285            }
286        }
287    }
288
289    return [expr {$result*$dxdydz}]
290}
291
292# eulerStep --
293#    Integrate a system of ordinary differential equations of the type
294#    x' = f(x,t), where x is a vector of quantities. Integration is
295#    done over a single step according to Euler's method.
296#
297# Arguments:
298#    t           Start value of independent variable (time for instance)
299#    tstep       Step size of interval
300#    xvec        Vector of dependent values at the start
301#    func        Function taking the arguments t and xvec to return
302#                the derivative of each dependent variable.
303# Return value:
304#    List of values at the end of the step
305#
306proc ::math::calculus::eulerStep { t tstep xvec func } {
307
308   set xderiv   [uplevel 1 $func $t [list $xvec]]
309   set result   {}
310   foreach xv $xvec dx $xderiv {
311      set xnew [expr {$xv+$tstep*$dx}]
312      lappend result $xnew
313   }
314
315   return $result
316}
317
318# heunStep --
319#    Integrate a system of ordinary differential equations of the type
320#    x' = f(x,t), where x is a vector of quantities. Integration is
321#    done over a single step according to Heun's method.
322#
323# Arguments:
324#    t           Start value of independent variable (time for instance)
325#    tstep       Step size of interval
326#    xvec        Vector of dependent values at the start
327#    func        Function taking the arguments t and xvec to return
328#                the derivative of each dependent variable.
329# Return value:
330#    List of values at the end of the step
331#
332proc ::math::calculus::heunStep { t tstep xvec func } {
333
334   #
335   # Predictor step
336   #
337   set funcq    [uplevel 1 namespace which -command $func]
338   set xpred    [eulerStep $t $tstep $xvec $funcq]
339
340   #
341   # Corrector step
342   #
343   set tcorr    [expr {$t+$tstep}]
344   set xcorr    [eulerStep $t     $tstep $xpred $funcq]
345
346   set result   {}
347   foreach xv $xvec xc $xcorr {
348      set xnew [expr {0.5*($xv+$xc)}]
349      lappend result $xnew
350   }
351
352   return $result
353}
354
355# rungeKuttaStep --
356#    Integrate a system of ordinary differential equations of the type
357#    x' = f(x,t), where x is a vector of quantities. Integration is
358#    done over a single step according to Runge-Kutta 4th order.
359#
360# Arguments:
361#    t           Start value of independent variable (time for instance)
362#    tstep       Step size of interval
363#    xvec        Vector of dependent values at the start
364#    func        Function taking the arguments t and xvec to return
365#                the derivative of each dependent variable.
366# Return value:
367#    List of values at the end of the step
368#
369proc ::math::calculus::rungeKuttaStep { t tstep xvec func } {
370
371   set funcq    [uplevel 1 namespace which -command $func]
372
373   #
374   # Four steps:
375   # - k1 = tstep*func(t,x0)
376   # - k2 = tstep*func(t+0.5*tstep,x0+0.5*k1)
377   # - k3 = tstep*func(t+0.5*tstep,x0+0.5*k2)
378   # - k4 = tstep*func(t+    tstep,x0+    k3)
379   # - x1 = x0 + (k1+2*k2+2*k3+k4)/6
380   #
381   set tstep2   [expr {$tstep/2.0}]
382   set tstep6   [expr {$tstep/6.0}]
383
384   set xk1      [$funcq $t $xvec]
385   set xvec2    {}
386   foreach x1 $xvec xv $xk1 {
387      lappend xvec2 [expr {$x1+$tstep2*$xv}]
388   }
389   set xk2      [$funcq [expr {$t+$tstep2}] $xvec2]
390
391   set xvec3    {}
392   foreach x1 $xvec xv $xk2 {
393      lappend xvec3 [expr {$x1+$tstep2*$xv}]
394   }
395   set xk3      [$funcq [expr {$t+$tstep2}] $xvec3]
396
397   set xvec4    {}
398   foreach x1 $xvec xv $xk3 {
399      lappend xvec4 [expr {$x1+$tstep*$xv}]
400   }
401   set xk4      [$funcq [expr {$t+$tstep}] $xvec4]
402
403   set result   {}
404   foreach x0 $xvec k1 $xk1 k2 $xk2 k3 $xk3 k4 $xk4 {
405      set dx [expr {$k1+2.0*$k2+2.0*$k3+$k4}]
406      lappend result [expr {$x0+$dx*$tstep6}]
407   }
408
409   return $result
410}
411
412# boundaryValueSecondOrder --
413#    Integrate a second-order differential equation and solve for
414#    given boundary values.
415#
416#    The equation is (see the documentation):
417#       d       dy   d
418#       -- A(x) -- + -- B(x) y + C(x) y = D(x)
419#       dx      dx   dx
420#
421#    The procedure uses finite differences and tridiagonal matrices to
422#    solve the equation. The boundary values are put in the matrix
423#    directly.
424#
425# Arguments:
426#    coeff_func  Name of triple-valued function for coefficients A, B, C
427#    force_func  Name of the function providing the force term D(x)
428#    leftbnd     Left boundary condition (list of: xvalue, boundary
429#                value or keyword zero-flux, zero-derivative)
430#    rightbnd    Right boundary condition (ditto)
431#    nostep      Number of steps
432# Return value:
433#    List of x-values and calculated values (x1, y1, x2, y2, ...)
434#
435proc ::math::calculus::boundaryValueSecondOrder {
436   coeff_func force_func leftbnd rightbnd nostep } {
437
438   set coeffq    [uplevel 1 namespace which -command $coeff_func]
439   set forceq    [uplevel 1 namespace which -command $force_func]
440
441   if { [llength $leftbnd] != 2 || [llength $rightbnd] != 2 } {
442      error "Boundary condition(s) incorrect"
443   }
444   if { $nostep < 1 } {
445      error "Number of steps must be larger/equal 1"
446   }
447
448   #
449   # Set up the matrix, as three different lists and the
450   # righthand side as the fourth
451   #
452   set xleft  [lindex $leftbnd 0]
453   set xright [lindex $rightbnd 0]
454   set xstep  [expr {($xright-$xleft)/double($nostep)}]
455
456   set acoeff {}
457   set bcoeff {}
458   set ccoeff {}
459   set dvalue {}
460
461   set x $xleft
462   foreach {A B C} [$coeffq $x] { break }
463
464   set A1 [expr {$A/$xstep-0.5*$B}]
465   set B1 [expr {$A/$xstep+0.5*$B+0.5*$C*$xstep}]
466   set C1 0.0
467
468   for { set i 1 } { $i <= $nostep } { incr i } {
469      set x [expr {$xleft+double($i)*$xstep}]
470      if { [expr {abs($x)-0.5*abs($xstep)}] < 0.0 } {
471         set x 0.0
472      }
473      foreach {A B C} [$coeffq $x] { break }
474
475      set A2 0.0
476      set B2 [expr {$A/$xstep-0.5*$B+0.5*$C*$xstep}]
477      set C2 [expr {$A/$xstep+0.5*$B}]
478      lappend acoeff [expr {$A1+$A2}]
479      lappend bcoeff [expr {-$B1-$B2}]
480      lappend ccoeff [expr {$C1+$C2}]
481      set A1 [expr {$A/$xstep-0.5*$B}]
482      set B1 [expr {$A/$xstep+0.5*$B+0.5*$C*$xstep}]
483      set C1 0.0
484   }
485   set xvec {}
486   for { set i 0 } { $i < $nostep } { incr i } {
487      set x [expr {$xleft+(0.5+double($i))*$xstep}]
488      if { [expr {abs($x)-0.25*abs($xstep)}] < 0.0 } {
489         set x 0.0
490      }
491      lappend xvec   $x
492      lappend dvalue [expr {$xstep*[$forceq $x]}]
493   }
494
495   #
496   # Substitute the boundary values
497   #
498   set A  [lindex $acoeff 0]
499   set D  [lindex $dvalue 0]
500   set D1 [expr {$D-$A*[lindex $leftbnd 1]}]
501   set C  [lindex $ccoeff end]
502   set D  [lindex $dvalue end]
503   set D2 [expr {$D-$C*[lindex $rightbnd 1]}]
504   set dvalue [concat $D1 [lrange $dvalue 1 end-1] $D2]
505
506   set yvec [solveTriDiagonal $acoeff $bcoeff $ccoeff $dvalue]
507
508   foreach x $xvec y $yvec {
509      lappend result $x $y
510   }
511   return $result
512}
513
514# solveTriDiagonal --
515#    Solve a system of equations Ax = b where A is a tridiagonal matrix
516#
517# Arguments:
518#    acoeff      Values on lower diagonal
519#    bcoeff      Values on main diagonal
520#    ccoeff      Values on upper diagonal
521#    dvalue      Values on righthand side
522# Return value:
523#    List of values forming the solution
524#
525proc ::math::calculus::solveTriDiagonal { acoeff bcoeff ccoeff dvalue } {
526
527   set nostep [llength $acoeff]
528   #
529   # First step: Gauss-elimination
530   #
531   set B [lindex $bcoeff 0]
532   set C [lindex $ccoeff 0]
533   set D [lindex $dvalue 0]
534   set bcoeff2 [list $B]
535   set dvalue2 [list $D]
536   for { set i 1 } { $i < $nostep } { incr i } {
537      set A2    [lindex $acoeff $i]
538      set B2    [lindex $bcoeff $i]
539      set C2    [lindex $ccoeff $i]
540      set D2    [lindex $dvalue $i]
541      set ratab [expr {$A2/$B}]
542      set B2    [expr {$B2-$ratab*$C}]
543      set D2    [expr {$D2-$ratab*$D}]
544      lappend bcoeff2 $B2
545      lappend dvalue2 $D2
546      set B     $B2
547      set D     $D2
548   }
549
550   #
551   # Second step: substitution
552   #
553   set yvec {}
554   set B [lindex $bcoeff2 end]
555   set D [lindex $dvalue2 end]
556   set y [expr {$D/$B}]
557   for { set i [expr {$nostep-2}] } { $i >= 0 } { incr i -1 } {
558      set yvec  [concat $y $yvec]
559      set B     [lindex $bcoeff2 $i]
560      set C     [lindex $ccoeff  $i]
561      set D     [lindex $dvalue2 $i]
562      set y     [expr {($D-$C*$y)/$B}]
563   }
564   set yvec [concat $y $yvec]
565
566   return $yvec
567}
568
569# newtonRaphson --
570#    Determine the root of an equation via the Newton-Raphson method
571#
572# Arguments:
573#    func        Function (proc) in x
574#    deriv       Derivative (proc) of func w.r.t. x
575#    initval     Initial value for x
576# Return value:
577#    Estimate of root
578#
579proc ::math::calculus::newtonRaphson { func deriv initval } {
580   variable nr_maxiter
581   variable nr_tolerance
582
583   set funcq  [uplevel 1 namespace which -command $func]
584   set derivq [uplevel 1 namespace which -command $deriv]
585
586   set value $initval
587   set diff  [expr {10.0*$nr_tolerance}]
588
589   for { set i 0 } { $i < $nr_maxiter } { incr i } {
590      if { $diff < $nr_tolerance } {
591         break
592      }
593
594      set newval [expr {$value-[$funcq $value]/[$derivq $value]}]
595      if { $value != 0.0 } {
596         set diff   [expr {abs($newval-$value)/abs($value)}]
597      } else {
598         set diff   [expr {abs($newval-$value)}]
599      }
600      set value $newval
601   }
602
603   return $newval
604}
605
606# newtonRaphsonParameters --
607#    Set the parameters for the Newton-Raphson method
608#
609# Arguments:
610#    maxiter     Maximum number of iterations
611#    tolerance   Relative precisiion of the result
612# Return value:
613#    None
614#
615proc ::math::calculus::newtonRaphsonParameters { maxiter tolerance } {
616   variable nr_maxiter
617   variable nr_tolerance
618
619   if { $maxiter > 0 } {
620      set nr_maxiter $maxiter
621   }
622   if { $tolerance > 0 } {
623      set nr_tolerance $tolerance
624   }
625}
626
627#----------------------------------------------------------------------
628#
629# midpoint --
630#
631#       Perform one set of steps in evaluating an integral using the
632#       midpoint method.
633#
634# Usage:
635#       midpoint f a b s ?n?
636#
637# Parameters:
638#       f - function to integrate
639#       a - One limit of integration
640#       b - Other limit of integration.  a and b need not be in ascending
641#           order.
642#       s - Value returned from a previous call to midpoint (see below)
643#       n - Step number (see below)
644#
645# Results:
646#       Returns an estimate of the integral obtained by dividing the
647#       interval into 3**n equal intervals and using the midpoint rule.
648#
649# Side effects:
650#       f is evaluated 2*3**(n-1) times and may have side effects.
651#
652# The 'midpoint' procedure is designed for successive approximations.
653# It should be called initially with n==0.  On this initial call, s
654# is ignored.  The function is evaluated at the midpoint of the interval, and
655# the value is multiplied by the width of the interval to give the
656# coarsest possible estimate of the integral.
657#
658# On each iteration except the first, n should be incremented by one,
659# and the previous value returned from [midpoint] should be supplied
660# as 's'.  The function will be evaluated at additional points
661# to give a total of 3**n equally spaced points, and the estimate
662# of the integral will be updated and returned
663#
664# Under normal circumstances, user code will not call this function
665# directly. Instead, it will use ::math::calculus::romberg to
666# do error control and extrapolation to a zero step size.
667#
668#----------------------------------------------------------------------
669
670proc ::math::calculus::midpoint { f a b { n 0 } { s 0. } } {
671
672    if { $n == 0 } {
673
674        # First iteration.  Simply evaluate the function at the midpoint
675        # of the interval.
676
677        set cmd $f; lappend cmd [expr { 0.5 * ( $a + $b ) }]; set v [eval $cmd]
678        return [expr { ( $b - $a ) * $v }]
679
680    } else {
681
682        # Subsequent iterations. We've divided the interval into
683        # $it subintervals.  Evaluate the function at the 1/3 and
684        # 2/3 points of each subinterval.  Then update the estimate
685        # of the integral that we produced on the last step with
686        # the new sum.
687
688        set it [expr { pow( 3, $n-1 ) }]
689        set h [expr { ( $b - $a ) / ( 3. * $it ) }]
690        set h2 [expr { $h + $h }]
691        set x [expr { $a + 0.5 * $h }]
692        set sum 0
693        for { set j 0 } { $j < $it } { incr j } {
694            set cmd $f; lappend cmd $x; set y [eval $cmd]
695            set sum [expr { $sum + $y }]
696            set x [expr { $x + $h2 }]
697            set cmd $f; lappend cmd $x; set y [eval $cmd]
698            set sum [expr { $sum + $y }]
699            set x [expr { $x + $h}]
700        }
701        return [expr { ( $s + ( $b - $a ) * $sum / $it ) / 3. }]
702
703    }
704}
705
706#----------------------------------------------------------------------
707#
708# romberg --
709#       
710#       Compute the integral of a function over an interval using
711#       Romberg's method.
712#
713# Usage:
714#       romberg f a b ?-option value?...
715#
716# Parameters:
717#       f - Function to integrate.  Must be a single Tcl command,
718#           to which will be appended the abscissa at which the function
719#           should be evaluated.  f should be analytic over the
720#           region of integration, but may have a removable singularity
721#           at either endpoint.
722#       a - One bound of the interval
723#       b - The other bound of the interval.  a and b need not be in
724#           ascending order.
725#
726# Options:
727#       -abserror ABSERROR
728#               Requests that the integration be performed to make
729#               the estimated absolute error of the integral less than
730#               the given value.  Default is 1.e-10.
731#       -relerror RELERROR
732#               Requests that the integration be performed to make
733#               the estimated absolute error of the integral less than
734#               the given value.  Default is 1.e-6.
735#       -degree N
736#               Specifies the degree of the polynomial that will be
737#               used to extrapolate to a zero step size.  -degree 0
738#               requests integration with the midpoint rule; -degree 1
739#               is equivalent to Simpson's 3/8 rule; higher degrees
740#               are difficult to describe but (within reason) give
741#               faster convergence for smooth functions.  Default is
742#               -degree 4.
743#       -maxiter N
744#               Specifies the maximum number of triplings of the
745#               number of steps to take in integration.  At most
746#               3**N function evaluations will be performed in
747#               integrating with -maxiter N.  The integration
748#               will terminate at that time, even if the result
749#               satisfies neither the -relerror nor -abserror tests.
750#
751# Results:
752#       Returns a two-element list.  The first element is the estimated
753#       value of the integral; the second is the estimated absolute
754#       error of the value.
755#
756#----------------------------------------------------------------------
757
758proc ::math::calculus::romberg { f a b args } {
759
760    # Replace f with a context-independent version
761
762    set f [lreplace $f 0 0 [uplevel 1 [list namespace which [lindex $f 0]]]]
763
764    # Assign default parameters
765
766    array set params {
767        -abserror 1.0e-10
768        -degree 4
769        -relerror 1.0e-6
770        -maxiter 14
771    }
772
773    # Extract parameters
774
775    if { ( [llength $args] % 2 ) != 0 } {
776        return -code error -errorcode [list romberg wrongNumArgs] \
777            "wrong \# args, should be\
778                 \"[lreplace [info level 0] 1 end \
779                         f x1 x2 ?-option value?...]\""
780    }
781    foreach { key value } $args {
782        if { ![info exists params($key)] } {
783            return -code error -errorcode [list romberg badoption $key] \
784                "unknown option \"$key\",\
785                     should be -abserror, -degree, -relerror, or -maxiter"
786        }
787        set params($key) $value
788    }
789
790    # Check params
791
792    if { ![string is double -strict $a] } {
793        return -code error [expectDouble $a]
794    }
795    if { ![string is double -strict $b] } {
796        return -code error [expectDouble $b]
797    }
798    if { ![string is double -strict $params(-abserror)] } {
799        return -code error [expectDouble $params(-abserror)]
800    }
801    if { ![string is integer -strict $params(-degree)] } {
802        return -code error [expectInteger $params(-degree)]
803    }
804    if { ![string is integer -strict $params(-maxiter)] } {
805        return -code error [expectInteger $params(-maxiter)]
806    }
807    if { ![string is double -strict $params(-relerror)] } {
808        return -code error [expectDouble $params(-relerror)]
809    }
810    foreach key {-abserror -degree -maxiter -relerror} {
811        if { $params($key) <= 0 } {
812            return -code error -errorcode [list romberg notPositive $key] \
813                "$key must be positive"
814        }
815    }
816    if { $params(-maxiter) <= $params(-degree) } {
817        return -code error -errorcode [list romberg tooFewIter] \
818            "-maxiter must be greater than -degree"
819    }
820
821    # Create lists of step size and sum with the given number of steps.
822
823    set x [list]
824    set y [list]
825    set s 0;                            # Current best estimate of integral
826    set indx end-$params(-degree)
827    set pow3 1.;                        # Current step size (times b-a)
828
829    # Perform successive integrations, tripling the number of steps each time
830
831    for { set i 0 } { $i < $params(-maxiter) } { incr i } {
832        set s [midpoint $f $a $b $i $s]
833        lappend x $pow3
834        lappend y $s
835        set pow3 [expr { $pow3 / 9. }]
836
837        # Once $degree steps have been done, start Richardson extrapolation
838        # to a zero step size.
839
840        if { $i >= $params(-degree) } {
841            set x [lrange $x $indx end]
842            set y [lrange $y $indx end]
843            foreach {estimate err} [neville $x $y 0.] break
844            if { $err < $params(-abserror)
845                 || $err < $params(-relerror) * abs($estimate) } {
846                return [list $estimate $err]
847            }
848        }
849    }
850
851    # If -maxiter iterations have been done, give up, and return
852    # with the current error estimate.
853
854    return [list $estimate $err]
855}
856
857#----------------------------------------------------------------------
858#
859# u_infinity --
860#       Change of variable for integrating over a half-infinite
861#       interval
862#
863# Parameters:
864#       f - Function being integrated
865#       u - 1/x, where x is the abscissa where f is to be evaluated
866#
867# Results:
868#       Returns f(1/u)/(u**2)
869#
870# Side effects:
871#       Whatever f does.
872#
873#----------------------------------------------------------------------
874
875proc ::math::calculus::u_infinity { f u } {
876    set cmd $f
877    lappend cmd [expr { 1.0 / $u }]
878    set y [eval $cmd]
879    return [expr { $y / ( $u * $u ) }]
880}
881
882#----------------------------------------------------------------------
883#
884# romberg_infinity --
885#       Evaluate a function on a half-open interval
886#
887# Usage:
888#       Same as 'romberg'
889#
890# The 'romberg_infinity' procedure performs Romberg integration on
891# an interval [a,b] where an infinite a or b may be represented by
892# a large number (e.g. 1.e30).  It operates by a change of variable;
893# instead of integrating f(x) from a to b, it makes a change
894# of variable u = 1/x, and integrates from 1/b to 1/a f(1/u)/u**2 du.
895#
896#----------------------------------------------------------------------
897
898proc ::math::calculus::romberg_infinity { f a b args } {
899    if { ![string is double -strict $a] } {
900        return -code error [expectDouble $a]
901    }
902    if { ![string is double -strict $b] } {
903        return -code error [expectDouble $b]
904    }
905    if { $a * $b <= 0. } {
906        return -code error -errorcode {romberg_infinity cross-axis} \
907            "limits of integration have opposite sign"
908    }
909    set f [lreplace $f 0 0 [uplevel 1 [list namespace which [lindex $f 0]]]]
910    set f [list u_infinity $f]
911    return [eval [linsert $args 0 \
912                      romberg $f [expr { 1.0 / $b }] [expr { 1.0 / $a }]]]
913}
914
915#----------------------------------------------------------------------
916#
917# u_sqrtSingLower --
918#       Change of variable for integrating over an interval with
919#       an inverse square root singularity at the lower bound.
920#
921# Parameters:
922#       f - Function being integrated
923#       a - Lower bound
924#       u - sqrt(x-a), where x is the abscissa where f is to be evaluated
925#
926# Results:
927#       Returns 2 * u * f( a + u**2 )
928#
929# Side effects:
930#       Whatever f does.
931#
932#----------------------------------------------------------------------
933
934proc ::math::calculus::u_sqrtSingLower { f a u } {
935    set cmd $f
936    lappend cmd [expr { $a + $u * $u }]
937    set y [eval $cmd]
938    return [expr { 2. * $u * $y }]
939}
940
941#----------------------------------------------------------------------
942#
943# u_sqrtSingUpper --
944#       Change of variable for integrating over an interval with
945#       an inverse square root singularity at the upper bound.
946#
947# Parameters:
948#       f - Function being integrated
949#       b - Upper bound
950#       u - sqrt(b-x), where x is the abscissa where f is to be evaluated
951#
952# Results:
953#       Returns 2 * u * f( b - u**2 )
954#
955# Side effects:
956#       Whatever f does.
957#
958#----------------------------------------------------------------------
959
960proc ::math::calculus::u_sqrtSingUpper { f b u } {
961    set cmd $f
962    lappend cmd [expr { $b - $u * $u }]
963    set y [eval $cmd]
964    return [expr { 2. * $u * $y }]
965}
966
967#----------------------------------------------------------------------
968#
969# math::calculus::romberg_sqrtSingLower --
970#       Integrate a function with an inverse square root singularity
971#       at the lower bound
972#
973# Usage:
974#       Same as 'romberg'
975#
976# The 'romberg_sqrtSingLower' procedure is a wrapper for 'romberg'
977# for integrating a function with an inverse square root singularity
978# at the lower bound of the interval.  It works by making the change
979# of variable u = sqrt( x-a ).
980#
981#----------------------------------------------------------------------
982
983proc ::math::calculus::romberg_sqrtSingLower { f a b args } {
984    if { ![string is double -strict $a] } {
985        return -code error [expectDouble $a]
986    }
987    if { ![string is double -strict $b] } {
988        return -code error [expectDouble $b]
989    }
990    if { $a >= $b } {
991        return -code error "limits of integration out of order"
992    }
993    set f [lreplace $f 0 0 [uplevel 1 [list namespace which [lindex $f 0]]]]
994    set f [list u_sqrtSingLower $f $a]
995    return [eval [linsert $args 0 \
996                             romberg $f 0 [expr { sqrt( $b - $a ) }]]]
997}
998
999#----------------------------------------------------------------------
1000#
1001# math::calculus::romberg_sqrtSingUpper --
1002#       Integrate a function with an inverse square root singularity
1003#       at the upper bound
1004#
1005# Usage:
1006#       Same as 'romberg'
1007#
1008# The 'romberg_sqrtSingUpper' procedure is a wrapper for 'romberg'
1009# for integrating a function with an inverse square root singularity
1010# at the upper bound of the interval.  It works by making the change
1011# of variable u = sqrt( b-x ).
1012#
1013#----------------------------------------------------------------------
1014
1015proc ::math::calculus::romberg_sqrtSingUpper { f a b args } {
1016    if { ![string is double -strict $a] } {
1017        return -code error [expectDouble $a]
1018    }
1019    if { ![string is double -strict $b] } {
1020        return -code error [expectDouble $b]
1021    }
1022    if { $a >= $b } {
1023        return -code error "limits of integration out of order"
1024    }
1025    set f [lreplace $f 0 0 [uplevel 1 [list namespace which [lindex $f 0]]]]
1026    set f [list u_sqrtSingUpper $f $b]
1027    return [eval [linsert $args 0 \
1028                      romberg $f 0. [expr { sqrt( $b - $a ) }]]]
1029}
1030
1031#----------------------------------------------------------------------
1032#
1033# u_powerLawLower --
1034#       Change of variable for integrating over an interval with
1035#       an integrable power law singularity at the lower bound.
1036#
1037# Parameters:
1038#       f - Function being integrated
1039#       gammaover1mgamma - gamma / (1 - gamma), where gamma is the power
1040#       oneover1mgamma - 1 / (1 - gamma), where gamma is the power
1041#       a - Lower limit of integration
1042#       u - Changed variable u == (x-a)**(1-gamma)
1043#
1044# Results:
1045#       Returns u**(1/1-gamma) * f(a + u**(1/1-gamma) ).
1046#
1047# Side effects:
1048#       Whatever f does.
1049#
1050#----------------------------------------------------------------------
1051
1052proc ::math::calculus::u_powerLawLower { f gammaover1mgamma oneover1mgamma
1053                                         a u } {
1054    set cmd $f
1055    lappend cmd [expr { $a + pow( $u, $oneover1mgamma ) }]
1056    set y [eval $cmd]
1057    return [expr { $y * pow( $u, $gammaover1mgamma ) }]
1058}
1059
1060#----------------------------------------------------------------------
1061#
1062# math::calculus::romberg_powerLawLower --
1063#       Integrate a function with an integrable power law singularity
1064#       at the lower bound
1065#
1066# Usage:
1067#       romberg_powerLawLower gamma f a b ?-option value...?
1068#
1069# Parameters:
1070#       gamma - Power (0<gamma<1) of the singularity
1071#       f - Function to integrate.  Must be a single Tcl command,
1072#           to which will be appended the abscissa at which the function
1073#           should be evaluated.  f is expected to have an integrable
1074#           power law singularity at the lower endpoint; that is, the
1075#           integrand is expected to diverge as (x-a)**gamma.
1076#       a - One bound of the interval
1077#       b - The other bound of the interval.  a and b need not be in
1078#           ascending order.
1079#
1080# Options:
1081#       -abserror ABSERROR
1082#               Requests that the integration be performed to make
1083#               the estimated absolute error of the integral less than
1084#               the given value.  Default is 1.e-10.
1085#       -relerror RELERROR
1086#               Requests that the integration be performed to make
1087#               the estimated absolute error of the integral less than
1088#               the given value.  Default is 1.e-6.
1089#       -degree N
1090#               Specifies the degree of the polynomial that will be
1091#               used to extrapolate to a zero step size.  -degree 0
1092#               requests integration with the midpoint rule; -degree 1
1093#               is equivalent to Simpson's 3/8 rule; higher degrees
1094#               are difficult to describe but (within reason) give
1095#               faster convergence for smooth functions.  Default is
1096#               -degree 4.
1097#       -maxiter N
1098#               Specifies the maximum number of triplings of the
1099#               number of steps to take in integration.  At most
1100#               3**N function evaluations will be performed in
1101#               integrating with -maxiter N.  The integration
1102#               will terminate at that time, even if the result
1103#               satisfies neither the -relerror nor -abserror tests.
1104#
1105# Results:
1106#       Returns a two-element list.  The first element is the estimated
1107#       value of the integral; the second is the estimated absolute
1108#       error of the value.
1109#
1110# The 'romberg_sqrtSingLower' procedure is a wrapper for 'romberg'
1111# for integrating a function with an integrable power law singularity
1112# at the lower bound of the interval.  It works by making the change
1113# of variable u = (x-a)**(1-gamma).
1114#
1115#----------------------------------------------------------------------
1116
1117proc ::math::calculus::romberg_powerLawLower { gamma f a b args } {
1118    if { ![string is double -strict $gamma] } {
1119        return -code error [expectDouble $gamma]
1120    }
1121    if { $gamma <= 0.0 || $gamma >= 1.0 } {
1122        return -code error -errorcode [list romberg gammaTooBig] \
1123            "gamma must lie in the interval (0,1)"
1124    }
1125    if { ![string is double -strict $a] } {
1126        return -code error [expectDouble $a]
1127    }
1128    if { ![string is double -strict $b] } {
1129        return -code error [expectDouble $b]
1130    }
1131    if { $a >= $b } {
1132        return -code error "limits of integration out of order"
1133    }
1134    set f [lreplace $f 0 0 [uplevel 1 [list namespace which [lindex $f 0]]]]
1135    set onemgamma [expr { 1. - $gamma }]
1136    set f [list u_powerLawLower $f \
1137               [expr { $gamma / $onemgamma }] \
1138               [expr { 1 / $onemgamma }] \
1139               $a]
1140       
1141    set limit [expr { pow( $b - $a, $onemgamma ) }]
1142    set result {}
1143    foreach v [eval [linsert $args 0 romberg $f 0 $limit]] {
1144        lappend result [expr { $v / $onemgamma }]
1145    }
1146    return $result
1147
1148}
1149
1150#----------------------------------------------------------------------
1151#
1152# u_powerLawLower --
1153#       Change of variable for integrating over an interval with
1154#       an integrable power law singularity at the upper bound.
1155#
1156# Parameters:
1157#       f - Function being integrated
1158#       gammaover1mgamma - gamma / (1 - gamma), where gamma is the power
1159#       oneover1mgamma - 1 / (1 - gamma), where gamma is the power
1160#       b - Upper limit of integration
1161#       u - Changed variable u == (b-x)**(1-gamma)
1162#
1163# Results:
1164#       Returns u**(1/1-gamma) * f(b-u**(1/1-gamma) ).
1165#
1166# Side effects:
1167#       Whatever f does.
1168#
1169#----------------------------------------------------------------------
1170
1171proc ::math::calculus::u_powerLawUpper { f gammaover1mgamma oneover1mgamma
1172                                         b u } {
1173    set cmd $f
1174    lappend cmd [expr { $b - pow( $u, $oneover1mgamma ) }]
1175    set y [eval $cmd]
1176    return [expr { $y * pow( $u, $gammaover1mgamma ) }]
1177}
1178
1179#----------------------------------------------------------------------
1180#
1181# math::calculus::romberg_powerLawUpper --
1182#       Integrate a function with an integrable power law singularity
1183#       at the upper bound
1184#
1185# Usage:
1186#       romberg_powerLawLower gamma f a b ?-option value...?
1187#
1188# Parameters:
1189#       gamma - Power (0<gamma<1) of the singularity
1190#       f - Function to integrate.  Must be a single Tcl command,
1191#           to which will be appended the abscissa at which the function
1192#           should be evaluated.  f is expected to have an integrable
1193#           power law singularity at the upper endpoint; that is, the
1194#           integrand is expected to diverge as (b-x)**gamma.
1195#       a - One bound of the interval
1196#       b - The other bound of the interval.  a and b need not be in
1197#           ascending order.
1198#
1199# Options:
1200#       -abserror ABSERROR
1201#               Requests that the integration be performed to make
1202#               the estimated absolute error of the integral less than
1203#               the given value.  Default is 1.e-10.
1204#       -relerror RELERROR
1205#               Requests that the integration be performed to make
1206#               the estimated absolute error of the integral less than
1207#               the given value.  Default is 1.e-6.
1208#       -degree N
1209#               Specifies the degree of the polynomial that will be
1210#               used to extrapolate to a zero step size.  -degree 0
1211#               requests integration with the midpoint rule; -degree 1
1212#               is equivalent to Simpson's 3/8 rule; higher degrees
1213#               are difficult to describe but (within reason) give
1214#               faster convergence for smooth functions.  Default is
1215#               -degree 4.
1216#       -maxiter N
1217#               Specifies the maximum number of triplings of the
1218#               number of steps to take in integration.  At most
1219#               3**N function evaluations will be performed in
1220#               integrating with -maxiter N.  The integration
1221#               will terminate at that time, even if the result
1222#               satisfies neither the -relerror nor -abserror tests.
1223#
1224# Results:
1225#       Returns a two-element list.  The first element is the estimated
1226#       value of the integral; the second is the estimated absolute
1227#       error of the value.
1228#
1229# The 'romberg_PowerLawUpper' procedure is a wrapper for 'romberg'
1230# for integrating a function with an integrable power law singularity
1231# at the upper bound of the interval.  It works by making the change
1232# of variable u = (b-x)**(1-gamma).
1233#
1234#----------------------------------------------------------------------
1235
1236proc ::math::calculus::romberg_powerLawUpper { gamma f a b args } {
1237    if { ![string is double -strict $gamma] } {
1238        return -code error [expectDouble $gamma]
1239    }
1240    if { $gamma <= 0.0 || $gamma >= 1.0 } {
1241        return -code error -errorcode [list romberg gammaTooBig] \
1242            "gamma must lie in the interval (0,1)"
1243    }
1244    if { ![string is double -strict $a] } {
1245        return -code error [expectDouble $a]
1246    }
1247    if { ![string is double -strict $b] } {
1248        return -code error [expectDouble $b]
1249    }
1250    if { $a >= $b } {
1251        return -code error "limits of integration out of order"
1252    }
1253    set f [lreplace $f 0 0 [uplevel 1 [list namespace which [lindex $f 0]]]]
1254    set onemgamma [expr { 1. - $gamma }]
1255    set f [list u_powerLawUpper $f \
1256               [expr { $gamma / $onemgamma }] \
1257               [expr { 1. / $onemgamma }] \
1258               $b]
1259       
1260    set limit [expr { pow( $b - $a, $onemgamma ) }]
1261    set result {}
1262    foreach v [eval [linsert $args 0 romberg $f 0 $limit]] {
1263        lappend result [expr { $v / $onemgamma }]
1264    }
1265    return $result
1266
1267}
1268
1269#----------------------------------------------------------------------
1270#
1271# u_expUpper --
1272#
1273#       Change of variable to integrate a function that decays
1274#       exponentially.
1275#
1276# Parameters:
1277#       f - Function to integrate
1278#       u - Changed variable u = exp(-x)
1279#
1280# Results:
1281#       Returns (1/u)*f(-log(u))
1282#
1283# Side effects:
1284#       Whatever f does.
1285#
1286#----------------------------------------------------------------------
1287
1288proc ::math::calculus::u_expUpper { f u } {
1289    set cmd $f
1290    lappend cmd [expr { -log($u) }]
1291    set y [eval $cmd]
1292    return [expr { $y / $u }]
1293}
1294
1295#----------------------------------------------------------------------
1296#
1297# romberg_expUpper --
1298#
1299#       Integrate a function that decays exponentially over a
1300#       half-infinite interval.
1301#
1302# Parameters:
1303#       Same as romberg.  The upper limit of integration, 'b',
1304#       is expected to be very large.
1305#
1306# Results:
1307#       Same as romberg.
1308#
1309# The romberg_expUpper function operates by making the change of
1310# variable, u = exp(-x).
1311#
1312#----------------------------------------------------------------------
1313
1314proc ::math::calculus::romberg_expUpper { f a b args } {
1315    if { ![string is double -strict $a] } {
1316        return -code error [expectDouble $a]
1317    }
1318    if { ![string is double -strict $b] } {
1319        return -code error [expectDouble $b]
1320    }
1321    if { $a >= $b } {
1322        return -code error "limits of integration out of order"
1323    }
1324    set f [lreplace $f 0 0 [uplevel 1 [list namespace which [lindex $f 0]]]]
1325    set f [list u_expUpper $f]
1326    return [eval [linsert $args 0 \
1327                      romberg $f [expr {exp(-$b)}] [expr {exp(-$a)}]]]
1328}
1329
1330#----------------------------------------------------------------------
1331#
1332# u_expLower --
1333#
1334#       Change of variable to integrate a function that grows
1335#       exponentially.
1336#
1337# Parameters:
1338#       f - Function to integrate
1339#       u - Changed variable u = exp(x)
1340#
1341# Results:
1342#       Returns (1/u)*f(log(u))
1343#
1344# Side effects:
1345#       Whatever f does.
1346#
1347#----------------------------------------------------------------------
1348
1349proc ::math::calculus::u_expLower { f u } {
1350    set cmd $f
1351    lappend cmd [expr { log($u) }]
1352    set y [eval $cmd]
1353    return [expr { $y / $u }]
1354}
1355
1356#----------------------------------------------------------------------
1357#
1358# romberg_expLower --
1359#
1360#       Integrate a function that grows exponentially over a
1361#       half-infinite interval.
1362#
1363# Parameters:
1364#       Same as romberg.  The lower limit of integration, 'a',
1365#       is expected to be very large and negative.
1366#
1367# Results:
1368#       Same as romberg.
1369#
1370# The romberg_expUpper function operates by making the change of
1371# variable, u = exp(x).
1372#
1373#----------------------------------------------------------------------
1374
1375proc ::math::calculus::romberg_expLower { f a b args } {
1376    if { ![string is double -strict $a] } {
1377        return -code error [expectDouble $a]
1378    }
1379    if { ![string is double -strict $b] } {
1380        return -code error [expectDouble $b]
1381    }
1382    if { $a >= $b } {
1383        return -code error "limits of integration out of order"
1384    }
1385    set f [lreplace $f 0 0 [uplevel 1 [list namespace which [lindex $f 0]]]]
1386    set f [list u_expLower $f]
1387    return [eval [linsert $args 0 \
1388                      romberg $f [expr {exp($a)}] [expr {exp($b)}]]]
1389}
1390
1391
1392# regula_falsi --
1393#    Compute the zero of a function via regula falsi
1394# Arguments:
1395#    f       Name of the procedure/command that evaluates the function
1396#    xb      Start of the interval that brackets the zero
1397#    xe      End of the interval that brackets the zero
1398#    eps     Relative error that is allowed (default: 1.0e-4)
1399# Result:
1400#    Estimate of the zero, such that the estimated (!)
1401#    error < eps * abs(xe-xb)
1402# Note:
1403#    f(xb)*f(xe) must be negative and eps must be positive
1404#
1405proc ::math::calculus::regula_falsi { f xb xe {eps 1.0e-4} } {
1406    if { $eps <= 0.0 } {
1407       return -code error "Relative error must be positive"
1408    }
1409
1410    set fb [$f $xb]
1411    set fe [$f $xe]
1412
1413    if { $fb * $fe > 0.0 } {
1414       return -code error "Interval must be chosen such that the \
1415function has a different sign at the beginning than at the end"
1416    }
1417
1418    set max_error [expr {$eps * abs($xe-$xb)}]
1419    set interval  [expr {abs($xe-$xb)}]
1420
1421    while { $interval > $max_error } {
1422       set coeff [expr {($fe-$fb)/($xe-$xb)}]
1423       set xi    [expr {$xb-$fb/$coeff}]
1424       set fi    [$f $xi]
1425
1426       if { $fi == 0.0 } {
1427          break
1428       }
1429       set diff1 [expr {abs($xe-$xi)}]
1430       set diff2 [expr {abs($xb-$xi)}]
1431       if { $diff1 > $diff2 } {
1432          set interval $diff2
1433       } else {
1434          set interval $diff1
1435       }
1436
1437       if { $fb*$fi < 0.0 } {
1438          set xe $xi
1439          set fe $fi
1440       } else {
1441          set xb $xi
1442          set fb $fi
1443       }
1444    }
1445
1446    return $xi
1447}
Note: See TracBrowser for help on using the repository browser.