calculus.tcl

Go to the documentation of this file.

/*  calculus.tcl --*/
/*     Package that implements several basic numerical methods, such*/
/*     as the integration of a one-dimensional function and the*/
/*     solution of a system of first-order differential equations.*/
/* */
/*  Copyright (c) 2002, 2003, 2004, 2006 by Arjen Markus.*/
/*  Copyright (c) 2004 by Kevin B. Kenny.  All rights reserved.*/
/*  See the file "license.terms" for information on usage and redistribution*/
/*  of this file, and for a DISCLAIMER OF ALL WARRANTIES.*/
/* */
/*  RCS: @(#) $Id: calculus.tcl,v 1.13 2006/03/28 20:36:00 arjenmarkus Exp $*/

package require Tcl 8.4
package require math::interpolate
package provide math::calculus 0.7

/*  math::calculus --*/
/*     Namespace for the commands*/

namespace ::math::calculus {

    namespace import ::math::interpolate::neville

    namespace import ::math::expectDouble ::math::expectInteger

    namespace export \
    integral integralExpr integral2D integral3D \
    eulerStep heunStep rungeKuttaStep           \
    boundaryValueSecondOrder solveTriDiagonal   \
    newtonRaphson newtonRaphsonParameters
    namespace export \
        integral2D_2accurate integral3D_accurate

    namespace export romberg romberg_infinity
    namespace export romberg_sqrtSingLower romberg_sqrtSingUpper
    namespace export romberg_powerLawLower romberg_powerLawUpper
    namespace export romberg_expLower romberg_expUpper

    namespace export regula_falsi

    variable nr_maxiter    20
    variable nr_tolerance   0.001

}

/*  integral --*/
/*     Integrate a function over a given interval using the Simpson rule*/
/* */
/*  Arguments:*/
/*     begin       Start of the interval*/
/*     end         End of the interval*/
/*     nosteps     Number of steps in which to divide the interval*/
/*     func        Name of the function to be integrated (takes one*/
/*                 argument)*/
/*  Return value:*/
/*     Computed integral*/
/* */
ret  ::math::calculus::integral ( type begin , type end , type nosteps , type func ) {

   set delta    [expr {($end-$begin)/double($nosteps)}]
   set hdelta   [expr {$delta/2.0}]
   set result   0.0
   set xval     $begin
   set func_end [uplevel 1 $func $xval]
   for { set i 1 } { $i <= $nosteps } { incr i } {
      set func_begin  $func_end
      set func_middle [uplevel 1 $func [expr {$xval+$hdelta}]]
      set func_end    [uplevel 1 $func [expr {$xval+$delta}]]
      set result      [expr  {$result+$func_begin+4.0*$func_middle+$func_end}]

      set xval        [expr {$begin+double($i)*$delta}]
   }

   return [expr {$result*$delta/6.0}]
}

/*  integralExpr --*/
/*     Integrate an expression with "x" as the integrate according to the*/
/*     Simpson rule*/
/* */
/*  Arguments:*/
/*     begin       Start of the interval*/
/*     end         End of the interval*/
/*     nosteps     Number of steps in which to divide the interval*/
/*     expression  Expression with "x" as the integrate*/
/*  Return value:*/
/*     Computed integral*/
/* */
ret  ::math::calculus::integralExpr ( type begin , type end , type nosteps , type expression ) {

   set delta    [expr {($end-$begin)/double($nosteps)}]
   set hdelta   [expr {$delta/2.0}]
   set result   0.0
   set x        $begin
   # FRINK: nocheck
   set func_end [expr $expression]
   for { set i 1 } { $i <= $nosteps } { incr i } {
      set func_begin  $func_end
      set x           [expr {$x+$hdelta}]
       # FRINK: nocheck
      set func_middle [expr $expression]
      set x           [expr {$x+$hdelta}]
       # FRINK: nocheck
      set func_end    [expr $expression]
      set result      [expr {$result+$func_begin+4.0*$func_middle+$func_end}]

      set x           [expr {$begin+double($i)*$delta}]
   }

   return [expr {$result*$delta/6.0}]
}

/*  integral2D --*/
/*     Integrate a given fucntion of two variables over a block,*/
/*     using bilinear interpolation (for this moment: block function)*/
/* */
/*  Arguments:*/
/*     xinterval   Start, stop and number of steps of the "x" interval*/
/*     yinterval   Start, stop and number of steps of the "y" interval*/
/*     func        Function of the two variables to be integrated*/
/*  Return value:*/
/*     Computed integral*/
/* */
ret  ::math::calculus::integral2D ( type xinterval , type yinterval , type func ) {

   foreach { xbegin xend xnumber } $xinterval { break }
   foreach { ybegin yend ynumber } $yinterval { break }

   set xdelta    [expr {($xend-$xbegin)/double($xnumber)}]
   set ydelta    [expr {($yend-$ybegin)/double($ynumber)}]
   set hxdelta   [expr {$xdelta/2.0}]
   set hydelta   [expr {$ydelta/2.0}]
   set result   0.0
   set dxdy      [expr {$xdelta*$ydelta}]
   for { set j 0 } { $j < $ynumber } { incr j } {
      set y [expr {$ybegin+$hydelta+double($j)*$ydelta}]
      for { set i 0 } { $i < $xnumber } { incr i } {
         set x           [expr {$xbegin+$hxdelta+double($i)*$xdelta}]
         set func_value  [uplevel 1 $func $x $y]
         set result      [expr {$result+$func_value}]
      }
   }

   return [expr {$result*$dxdy}]
}

/*  integral3D --*/
/*     Integrate a given fucntion of two variables over a block,*/
/*     using trilinear interpolation (for this moment: block function)*/
/* */
/*  Arguments:*/
/*     xinterval   Start, stop and number of steps of the "x" interval*/
/*     yinterval   Start, stop and number of steps of the "y" interval*/
/*     zinterval   Start, stop and number of steps of the "z" interval*/
/*     func        Function of the three variables to be integrated*/
/*  Return value:*/
/*     Computed integral*/
/* */
ret  ::math::calculus::integral3D ( type xinterval , type yinterval , type zinterval , type func ) {

   foreach { xbegin xend xnumber } $xinterval { break }
   foreach { ybegin yend ynumber } $yinterval { break }
   foreach { zbegin zend znumber } $zinterval { break }

   set xdelta    [expr {($xend-$xbegin)/double($xnumber)}]
   set ydelta    [expr {($yend-$ybegin)/double($ynumber)}]
   set zdelta    [expr {($zend-$zbegin)/double($znumber)}]
   set hxdelta   [expr {$xdelta/2.0}]
   set hydelta   [expr {$ydelta/2.0}]
   set hzdelta   [expr {$zdelta/2.0}]
   set result   0.0
   set dxdydz    [expr {$xdelta*$ydelta*$zdelta}]
   for { set k 0 } { $k < $znumber } { incr k } {
      set z [expr {$zbegin+$hzdelta+double($k)*$zdelta}]
      for { set j 0 } { $j < $ynumber } { incr j } {
         set y [expr {$ybegin+$hydelta+double($j)*$ydelta}]
         for { set i 0 } { $i < $xnumber } { incr i } {
            set x           [expr {$xbegin+$hxdelta+double($i)*$xdelta}]
            set func_value  [uplevel 1 $func $x $y $z]
            set result      [expr {$result+$func_value}]
         }
      }
   }

   return [expr {$result*$dxdydz}]
}

/*  integral2D_accurate --*/
/*     Integrate a given function of two variables over a block,*/
/*     using a four-point quadrature formula*/
/* */
/*  Arguments:*/
/*     xinterval   Start, stop and number of steps of the "x" interval*/
/*     yinterval   Start, stop and number of steps of the "y" interval*/
/*     func        Function of the two variables to be integrated*/
/*  Return value:*/
/*     Computed integral*/
/* */
ret  ::math::calculus::integral2D_accurate ( type xinterval , type yinterval , type func ) {

    foreach { xbegin xend xnumber } $xinterval { break }
    foreach { ybegin yend ynumber } $yinterval { break }

    set alpha     [expr {sqrt(2.0/3.0)}]
    set minalpha  [expr {-$alpha}]
    set dpoints   [list $alpha 0.0 $minalpha 0.0 0.0 $alpha 0.0 $minalpha]

    set xdelta    [expr {($xend-$xbegin)/double($xnumber)}]
    set ydelta    [expr {($yend-$ybegin)/double($ynumber)}]
    set hxdelta   [expr {$xdelta/2.0}]
    set hydelta   [expr {$ydelta/2.0}]
    set result   0.0
    set dxdy      [expr {0.25*$xdelta*$ydelta}]

    for { set j 0 } { $j < $ynumber } { incr j } {
        set y [expr {$ybegin+$hydelta+double($j)*$ydelta}]
        for { set i 0 } { $i < $xnumber } { incr i } {
            set x [expr {$xbegin+$hxdelta+double($i)*$xdelta}]

            foreach {dx dy} $dpoints {
                set x1          [expr {$x+$dx}]
                set y1          [expr {$y+$dy}]
                set func_value  [uplevel 1 $func $x1 $y1]
                set result      [expr {$result+$func_value}]
            }
        }
    }

    return [expr {$result*$dxdy}]
}

/*  integral3D_accurate --*/
/*     Integrate a given function of three variables over a block,*/
/*     using an 8-point quadrature formula*/
/* */
/*  Arguments:*/
/*     xinterval   Start, stop and number of steps of the "x" interval*/
/*     yinterval   Start, stop and number of steps of the "y" interval*/
/*     zinterval   Start, stop and number of steps of the "z" interval*/
/*     func        Function of the three variables to be integrated*/
/*  Return value:*/
/*     Computed integral*/
/* */
ret  ::math::calculus::integral3D_accurate ( type xinterval , type yinterval , type zinterval , type func ) {

    foreach { xbegin xend xnumber } $xinterval { break }
    foreach { ybegin yend ynumber } $yinterval { break }
    foreach { zbegin zend znumber } $zinterval { break }

    set alpha     [expr {sqrt(1.0/3.0)}]
    set minalpha  [expr {-$alpha}]

    set dpoints   [list $alpha    $alpha    $alpha    \
                        $alpha    $alpha    $minalpha \
                        $alpha    $minalpha $alpha    \
                        $alpha    $minalpha $minalpha \
                        $minalpha $alpha    $alpha    \
                        $minalpha $alpha    $minalpha \
                        $minalpha $minalpha $alpha    \
                        $minalpha $minalpha $minalpha ]

    set xdelta    [expr {($xend-$xbegin)/double($xnumber)}]
    set ydelta    [expr {($yend-$ybegin)/double($ynumber)}]
    set zdelta    [expr {($zend-$zbegin)/double($znumber)}]
    set hxdelta   [expr {$xdelta/2.0}]
    set hydelta   [expr {$ydelta/2.0}]
    set hzdelta   [expr {$zdelta/2.0}]
    set result    0.0
    set dxdydz    [expr {0.125*$xdelta*$ydelta*$zdelta}]

    for { set k 0 } { $k < $znumber } { incr k } {
        set z [expr {$zbegin+$hzdelta+double($k)*$zdelta}]
        for { set j 0 } { $j < $ynumber } { incr j } {
            set y [expr {$ybegin+$hydelta+double($j)*$ydelta}]
            for { set i 0 } { $i < $xnumber } { incr i } {
                set x [expr {$xbegin+$hxdelta+double($i)*$xdelta}]

                foreach {dx dy dz} $dpoints {
                    set x1 [expr {$x+$dx}]
                    set y1 [expr {$y+$dy}]
                    set z1 [expr {$z+$dz}]
                    set func_value  [uplevel 1 $func $x1 $y1 $z1]
                    set result      [expr {$result+$func_value}]
                }
            }
        }
    }

    return [expr {$result*$dxdydz}]
}

/*  eulerStep --*/
/*     Integrate a system of ordinary differential equations of the type*/
/*     x' = f(x,t), where x is a vector of quantities. Integration is*/
/*     done over a single step according to Euler's method.*/
/* */
/*  Arguments:*/
/*     t           Start value of independent variable (time for instance)*/
/*     tstep       Step size of interval*/
/*     xvec        Vector of dependent values at the start*/
/*     func        Function taking the arguments t and xvec to return*/
/*                 the derivative of each dependent variable.*/
/*  Return value:*/
/*     List of values at the end of the step*/
/* */
ret  ::math::calculus::eulerStep ( type t , type tstep , type xvec , type func ) {

   set xderiv   [uplevel 1 $func $t [list $xvec]]
   set result   {}
   foreach xv $xvec dx $xderiv {
      set xnew [expr {$xv+$tstep*$dx}]
      lappend result $xnew
   }

   return $result
}

/*  heunStep --*/
/*     Integrate a system of ordinary differential equations of the type*/
/*     x' = f(x,t), where x is a vector of quantities. Integration is*/
/*     done over a single step according to Heun's method.*/
/* */
/*  Arguments:*/
/*     t           Start value of independent variable (time for instance)*/
/*     tstep       Step size of interval*/
/*     xvec        Vector of dependent values at the start*/
/*     func        Function taking the arguments t and xvec to return*/
/*                 the derivative of each dependent variable.*/
/*  Return value:*/
/*     List of values at the end of the step*/
/* */
ret  ::math::calculus::heunStep ( type t , type tstep , type xvec , type func ) {

   #
   # Predictor step
   #
   set funcq    [uplevel 1 namespace which -command $func]
   set xpred    [eulerStep $t $tstep $xvec $funcq]

   #
   # Corrector step
   #
   set tcorr    [expr {$t+$tstep}]
   set xcorr    [eulerStep $t     $tstep $xpred $funcq]

   set result   {}
   foreach xv $xvec xc $xcorr {
      set xnew [expr {0.5*($xv+$xc)}]
      lappend result $xnew
   }

   return $result
}

/*  rungeKuttaStep --*/
/*     Integrate a system of ordinary differential equations of the type*/
/*     x' = f(x,t), where x is a vector of quantities. Integration is*/
/*     done over a single step according to Runge-Kutta 4th order.*/
/* */
/*  Arguments:*/
/*     t           Start value of independent variable (time for instance)*/
/*     tstep       Step size of interval*/
/*     xvec        Vector of dependent values at the start*/
/*     func        Function taking the arguments t and xvec to return*/
/*                 the derivative of each dependent variable.*/
/*  Return value:*/
/*     List of values at the end of the step*/
/* */
ret  ::math::calculus::rungeKuttaStep ( type t , type tstep , type xvec , type func ) {

   set funcq    [uplevel 1 namespace which -command $func]

   #
   # Four steps:
   # - k1 = tstep*func(t,x0)
   # - k2 = tstep*func(t+0.5*tstep,x0+0.5*k1)
   # - k3 = tstep*func(t+0.5*tstep,x0+0.5*k2)
   # - k4 = tstep*func(t+    tstep,x0+    k3)
   # - x1 = x0 + (k1+2*k2+2*k3+k4)/6
   #
   set tstep2   [expr {$tstep/2.0}]
   set tstep6   [expr {$tstep/6.0}]

   set xk1      [$funcq $t $xvec]
   set xvec2    {}
   foreach x1 $xvec xv $xk1 {
      lappend xvec2 [expr {$x1+$tstep2*$xv}]
   }
   set xk2      [$funcq [expr {$t+$tstep2}] $xvec2]

   set xvec3    {}
   foreach x1 $xvec xv $xk2 {
      lappend xvec3 [expr {$x1+$tstep2*$xv}]
   }
   set xk3      [$funcq [expr {$t+$tstep2}] $xvec3]

   set xvec4    {}
   foreach x1 $xvec xv $xk3 {
      lappend xvec4 [expr {$x1+$tstep*$xv}]
   }
   set xk4      [$funcq [expr {$t+$tstep}] $xvec4]

   set result   {}
   foreach x0 $xvec k1 $xk1 k2 $xk2 k3 $xk3 k4 $xk4 {
      set dx [expr {$k1+2.0*$k2+2.0*$k3+$k4}]
      lappend result [expr {$x0+$dx*$tstep6}]
   }

   return $result
}

/*  boundaryValueSecondOrder --*/
/*     Integrate a second-order differential equation and solve for*/
/*     given boundary values.*/
/* */
/*     The equation is (see the documentation):*/
/*        d       dy   d*/
/*        -- A(x) -- + -- B(x) y + C(x) y = D(x)*/
/*        dx      dx   dx*/
/* */
/*     The procedure uses finite differences and tridiagonal matrices to*/
/*     solve the equation. The boundary values are put in the matrix*/
/*     directly.*/
/* */
/*  Arguments:*/
/*     coeff_func  Name of triple-valued function for coefficients A, B, C*/
/*     force_func  Name of the function providing the force term D(x)*/
/*     leftbnd     Left boundary condition (list of: xvalue, boundary*/
/*                 value or keyword zero-flux, zero-derivative)*/
/*     rightbnd    Right boundary condition (ditto)*/
/*     nostep      Number of steps*/
/*  Return value:*/
/*     List of x-values and calculated values (x1, y1, x2, y2, ...)*/
/* */
ret  ::math::calculus::boundaryValueSecondOrder (
   type coeff_, type func , type force_, type func , type leftbnd , type rightbnd , type nostep ) {

   set coeffq    [uplevel 1 namespace which -command $coeff_func]
   set forceq    [uplevel 1 namespace which -command $force_func]

   if { [llength $leftbnd] != 2 || [llength $rightbnd] != 2 } {
      error "Boundary condition(s) incorrect"
   }
   if { $nostep < 1 } {
      error "Number of steps must be larger/equal 1"
   }

   #
   # Set up the matrix, as three different lists and the
   # righthand side as the fourth
   #
   set xleft  [lindex $leftbnd 0]
   set xright [lindex $rightbnd 0]
   set xstep  [expr {($xright-$xleft)/double($nostep)}]

   set acoeff {}
   set bcoeff {}
   set ccoeff {}
   set dvalue {}

   set x $xleft
   foreach {A B C} [$coeffq $x] { break }

   set A1 [expr {$A/$xstep-0.5*$B}]
   set B1 [expr {$A/$xstep+0.5*$B+0.5*$C*$xstep}]
   set C1 0.0

   for { set i 1 } { $i <= $nostep } { incr i } {
      set x [expr {$xleft+double($i)*$xstep}]
      if { [expr {abs($x)-0.5*abs($xstep)}] < 0.0 } {
         set x 0.0
      }
      foreach {A B C} [$coeffq $x] { break }

      set A2 0.0
      set B2 [expr {$A/$xstep-0.5*$B+0.5*$C*$xstep}]
      set C2 [expr {$A/$xstep+0.5*$B}]
      lappend acoeff [expr {$A1+$A2}]
      lappend bcoeff [expr {-$B1-$B2}]
      lappend ccoeff [expr {$C1+$C2}]
      set A1 [expr {$A/$xstep-0.5*$B}]
      set B1 [expr {$A/$xstep+0.5*$B+0.5*$C*$xstep}]
      set C1 0.0
   }
   set xvec {}
   for { set i 0 } { $i < $nostep } { incr i } {
      set x [expr {$xleft+(0.5+double($i))*$xstep}]
      if { [expr {abs($x)-0.25*abs($xstep)}] < 0.0 } {
         set x 0.0
      }
      lappend xvec   $x
      lappend dvalue [expr {$xstep*[$forceq $x]}]
   }

   #
   # Substitute the boundary values
   #
   set A  [lindex $acoeff 0]
   set D  [lindex $dvalue 0]
   set D1 [expr {$D-$A*[lindex $leftbnd 1]}]
   set C  [lindex $ccoeff end]
   set D  [lindex $dvalue end]
   set D2 [expr {$D-$C*[lindex $rightbnd 1]}]
   set dvalue [concat $D1 [lrange $dvalue 1 end-1] $D2]

   set yvec [solveTriDiagonal $acoeff $bcoeff $ccoeff $dvalue]

   foreach x $xvec y $yvec {
      lappend result $x $y
   }
   return $result
}

/*  solveTriDiagonal --*/
/*     Solve a system of equations Ax = b where A is a tridiagonal matrix*/
/* */
/*  Arguments:*/
/*     acoeff      Values on lower diagonal*/
/*     bcoeff      Values on main diagonal*/
/*     ccoeff      Values on upper diagonal*/
/*     dvalue      Values on righthand side*/
/*  Return value:*/
/*     List of values forming the solution*/
/* */
ret  ::math::calculus::solveTriDiagonal ( type acoeff , type bcoeff , type ccoeff , type dvalue ) {

   set nostep [llength $acoeff]
   #
   # First step: Gauss-elimination
   #
   set B [lindex $bcoeff 0]
   set C [lindex $ccoeff 0]
   set D [lindex $dvalue 0]
   set bcoeff2 [list $B]
   set dvalue2 [list $D]
   for { set i 1 } { $i < $nostep } { incr i } {
      set A2    [lindex $acoeff $i]
      set B2    [lindex $bcoeff $i]
      set C2    [lindex $ccoeff $i]
      set D2    [lindex $dvalue $i]
      set ratab [expr {$A2/$B}]
      set B2    [expr {$B2-$ratab*$C}]
      set D2    [expr {$D2-$ratab*$D}]
      lappend bcoeff2 $B2
      lappend dvalue2 $D2
      set B     $B2
      set D     $D2
   }

   #
   # Second step: substitution
   #
   set yvec {}
   set B [lindex $bcoeff2 end]
   set D [lindex $dvalue2 end]
   set y [expr {$D/$B}]
   for { set i [expr {$nostep-2}] } { $i >= 0 } { incr i -1 } {
      set yvec  [concat $y $yvec]
      set B     [lindex $bcoeff2 $i]
      set C     [lindex $ccoeff  $i]
      set D     [lindex $dvalue2 $i]
      set y     [expr {($D-$C*$y)/$B}]
   }
   set yvec [concat $y $yvec]

   return $yvec
}

/*  newtonRaphson --*/
/*     Determine the root of an equation via the Newton-Raphson method*/
/* */
/*  Arguments:*/
/*     func        Function (proc) in x*/
/*     deriv       Derivative (proc) of func w.r.t. x*/
/*     initval     Initial value for x*/
/*  Return value:*/
/*     Estimate of root*/
/* */
ret  ::math::calculus::newtonRaphson ( type func , type deriv , type initval ) {
   variable nr_maxiter
   variable nr_tolerance

   set funcq  [uplevel 1 namespace which -command $func]
   set derivq [uplevel 1 namespace which -command $deriv]

   set value $initval
   set diff  [expr {10.0*$nr_tolerance}]

   for { set i 0 } { $i < $nr_maxiter } { incr i } {
      if { $diff < $nr_tolerance } {
         break
      }

      set newval [expr {$value-[$funcq $value]/[$derivq $value]}]
      if { $value != 0.0 } {
         set diff   [expr {abs($newval-$value)/abs($value)}]
      } else {
         set diff   [expr {abs($newval-$value)}]
      }
      set value $newval
   }

   return $newval
}

/*  newtonRaphsonParameters --*/
/*     Set the parameters for the Newton-Raphson method*/
/* */
/*  Arguments:*/
/*     maxiter     Maximum number of iterations*/
/*     tolerance   Relative precisiion of the result*/
/*  Return value:*/
/*     None*/
/* */
ret  ::math::calculus::newtonRaphsonParameters ( type maxiter , type tolerance ) {
   variable nr_maxiter
   variable nr_tolerance

   if { $maxiter > 0 } {
      set nr_maxiter $maxiter
   }
   if { $tolerance > 0 } {
      set nr_tolerance $tolerance
   }
}

/* ----------------------------------------------------------------------*/
/* */
/*  midpoint --*/
/* */
/*  Perform one set of steps in evaluating an integral using the*/
/*  midpoint method.*/
/* */
/*  Usage:*/
/*  midpoint f a b s ?n?*/
/* */
/*  Parameters:*/
/*  f - function to integrate*/
/*  a - One limit of integration*/
/*  b - Other limit of integration.  a and b need not be in ascending*/
/*      order.*/
/*  s - Value returned from a previous call to midpoint (see below)*/
/*  n - Step number (see below)*/
/* */
/*  Results:*/
/*  Returns an estimate of the integral obtained by dividing the*/
/*  interval into 3**n equal intervals and using the midpoint rule.*/
/* */
/*  Side effects:*/
/*  f is evaluated 2*3**(n-1) times and may have side effects.*/
/* */
/*  The 'midpoint' procedure is designed for successive approximations.*/
/*  It should be called initially with n==0.  On this initial call, s*/
/*  is ignored.  The function is evaluated at the midpoint of the interval, and*/
/*  the value is multiplied by the width of the interval to give the*/
/*  coarsest possible estimate of the integral.*/
/* */
/*  On each iteration except the first, n should be incremented by one,*/
/*  and the previous value returned from [midpoint] should be supplied*/
/*  as 's'.  The function will be evaluated at additional points*/
/*  to give a total of 3**n equally spaced points, and the estimate*/
/*  of the integral will be updated and returned*/
/* */
/*  Under normal circumstances, user code will not call this function*/
/*  directly. Instead, it will use ::math::calculus::romberg to*/
/*  do error control and extrapolation to a zero step size.*/
/* */
/* ----------------------------------------------------------------------*/

ret  ::math::calculus::midpoint ( type f , type a , type b , optional n =0  , optional s =0.  ) {

    if { $n == 0 } {

    # First iteration.  Simply evaluate the function at the midpoint
    # of the interval.

    set cmd $f; lappend cmd [expr { 0.5 * ( $a + $b ) }]; set v [eval $cmd]
    return [expr { ( $b - $a ) * $v }]

    } else {

    # Subsequent iterations. We've divided the interval into
    # $it subintervals.  Evaluate the function at the 1/3 and
    # 2/3 points of each subinterval.  Then update the estimate
    # of the integral that we produced on the last step with
    # the new sum.

    set it [expr { pow( 3, $n-1 ) }]
    set h [expr { ( $b - $a ) / ( 3. * $it ) }]
    set h2 [expr { $h + $h }]
    set x [expr { $a + 0.5 * $h }]
    set sum 0
    for { set j 0 } { $j < $it } { incr j } {
        set cmd $f; lappend cmd $x; set y [eval $cmd]
        set sum [expr { $sum + $y }]
        set x [expr { $x + $h2 }]
        set cmd $f; lappend cmd $x; set y [eval $cmd]
        set sum [expr { $sum + $y }]
        set x [expr { $x + $h}]
    }
    return [expr { ( $s + ( $b - $a ) * $sum / $it ) / 3. }]

    }
}

/* ----------------------------------------------------------------------*/
/* */
/*  romberg --*/
/*  */
/*  Compute the integral of a function over an interval using*/
/*  Romberg's method.*/
/* */
/*  Usage:*/
/*  romberg f a b ?-option value?...*/
/* */
/*  Parameters:*/
/*  f - Function to integrate.  Must be a single Tcl command,*/
/*      to which will be appended the abscissa at which the function*/
/*      should be evaluated.  f should be analytic over the*/
/*      region of integration, but may have a removable singularity*/
/*      at either endpoint.*/
/*  a - One bound of the interval*/
/*  b - The other bound of the interval.  a and b need not be in*/
/*      ascending order.*/
/* */
/*  Options:*/
/*  -abserror ABSERROR*/
/*      Requests that the integration be performed to make*/
/*      the estimated absolute error of the integral less than*/
/*      the given value.  Default is 1.e-10.*/
/*  -relerror RELERROR*/
/*      Requests that the integration be performed to make*/
/*      the estimated absolute error of the integral less than*/
/*      the given value.  Default is 1.e-6.*/
/*  -degree N*/
/*      Specifies the degree of the polynomial that will be*/
/*      used to extrapolate to a zero step size.  -degree 0*/
/*      requests integration with the midpoint rule; -degree 1*/
/*      is equivalent to Simpson's 3/8 rule; higher degrees*/
/*      are difficult to describe but (within reason) give*/
/*      faster convergence for smooth functions.  Default is*/
/*      -degree 4.*/
/*  -maxiter N*/
/*      Specifies the maximum number of triplings of the*/
/*      number of steps to take in integration.  At most*/
/*      3**N function evaluations will be performed in*/
/*      integrating with -maxiter N.  The integration*/
/*      will terminate at that time, even if the result*/
/*      satisfies neither the -relerror nor -abserror tests.*/
/* */
/*  Results:*/
/*  Returns a two-element list.  The first element is the estimated*/
/*  value of the integral; the second is the estimated absolute*/
/*  error of the value.*/
/* */
/* ----------------------------------------------------------------------*/

ret  ::math::calculus::romberg ( type f , type a , type b , type args ) {

    # Replace f with a context-independent version

    set f [lreplace $f 0 0 [uplevel 1 [list namespace which [lindex $f 0]]]]

    # Assign default parameters

    array set params {
    -abserror 1.0e-10
    -degree 4
    -relerror 1.0e-6
    -maxiter 14
    }

    # Extract parameters

    if { ( [llength $args] % 2 ) != 0 } {
        return -code error -errorcode [list romberg wrongNumArgs] \
            "wrong \# args, should be\
                 \"[lreplace [info level 0] 1 end \
                         f x1 x2 ?-option value?...]\""
    }
    foreach { key value } $args {
        if { ![info exists params($key)] } {
            return -code error -errorcode [list romberg badoption $key] \
                "unknown option \"$key\",\
                     should be -abserror, -degree, -relerror, or -maxiter"
        }
        set params($key) $value
    }

    # Check params

    if { ![string is double -strict $a] } {
    return -code error [expectDouble $a]
    }
    if { ![string is double -strict $b] } {
    return -code error [expectDouble $b]
    }
    if { ![string is double -strict $params(-abserror)] } {
    return -code error [expectDouble $params(-abserror)]
    }
    if { ![string is integer -strict $params(-degree)] } {
    return -code error [expectInteger $params(-degree)]
    }
    if { ![string is integer -strict $params(-maxiter)] } {
    return -code error [expectInteger $params(-maxiter)]
    }
    if { ![string is double -strict $params(-relerror)] } {
    return -code error [expectDouble $params(-relerror)]
    }
    foreach key {-abserror -degree -maxiter -relerror} {
    if { $params($key) <= 0 } {
        return -code error -errorcode [list romberg notPositive $key] \
        "$key must be positive"
    }
    }
    if { $params(-maxiter) <= $params(-degree) } {
    return -code error -errorcode [list romberg tooFewIter] \
        "-maxiter must be greater than -degree"
    }

    # Create lists of step size and sum with the given number of steps.

    set x [list]
    set y [list]
    set s 0;                # Current best estimate of integral
    set indx end-$params(-degree)
    set pow3 1.;            # Current step size (times b-a)

    # Perform successive integrations, tripling the number of steps each time

    for { set i 0 } { $i < $params(-maxiter) } { incr i } {
    set s [midpoint $f $a $b $i $s]
    lappend x $pow3
    lappend y $s
    set pow3 [expr { $pow3 / 9. }]

    # Once $degree steps have been done, start Richardson extrapolation
    # to a zero step size.

    if { $i >= $params(-degree) } {
        set x [lrange $x $indx end]
        set y [lrange $y $indx end]
        foreach {estimate err} [neville $x $y 0.] break
        if { $err < $params(-abserror)
         || $err < $params(-relerror) * abs($estimate) } {
        return [list $estimate $err]
        }
    }
    }

    # If -maxiter iterations have been done, give up, and return
    # with the current error estimate.

    return [list $estimate $err]
}

/* ----------------------------------------------------------------------*/
/* */
/*  u_infinity --*/
/*  Change of variable for integrating over a half-infinite*/
/*  interval*/
/* */
/*  Parameters:*/
/*  f - Function being integrated*/
/*  u - 1/x, where x is the abscissa where f is to be evaluated*/
/* */
/*  Results:*/
/*  Returns f(1/u)/(u**2)*/
/* */
/*  Side effects:*/
/*  Whatever f does.*/
/* */
/* ----------------------------------------------------------------------*/

ret  ::math::calculus::u_infinity ( type f , type u ) {
    set cmd $f
    lappend cmd [expr { 1.0 / $u }]
    set y [eval $cmd]
    return [expr { $y / ( $u * $u ) }]
}

/* ----------------------------------------------------------------------*/
/* */
/*  romberg_infinity --*/
/*  Evaluate a function on a half-open interval*/
/* */
/*  Usage:*/
/*  Same as 'romberg'*/
/* */
/*  The 'romberg_infinity' procedure performs Romberg integration on*/
/*  an interval [a,b] where an infinite a or b may be represented by*/
/*  a large number (e.g. 1.e30).  It operates by a change of variable;*/
/*  instead of integrating f(x) from a to b, it makes a change*/
/*  of variable u = 1/x, and integrates from 1/b to 1/a f(1/u)/u**2 du.*/
/* */
/* ----------------------------------------------------------------------*/

ret  ::math::calculus::romberg_infinity ( type f , type a , type b , type args ) {
    if { ![string is double -strict $a] } {
    return -code error [expectDouble $a]
    }
    if { ![string is double -strict $b] } {
    return -code error [expectDouble $b]
    }
    if { $a * $b <= 0. } {
        return -code error -errorcode {romberg_infinity cross-axis} \
            "limits of integration have opposite sign"
    }
    set f [lreplace $f 0 0 [uplevel 1 [list namespace which [lindex $f 0]]]]
    set f [list u_infinity $f]
    return [eval [linsert $args 0 \
                      romberg $f [expr { 1.0 / $b }] [expr { 1.0 / $a }]]]
}

/* ----------------------------------------------------------------------*/
/* */
/*  u_sqrtSingLower --*/
/*  Change of variable for integrating over an interval with*/
/*  an inverse square root singularity at the lower bound.*/
/* */
/*  Parameters:*/
/*  f - Function being integrated*/
/*  a - Lower bound*/
/*  u - sqrt(x-a), where x is the abscissa where f is to be evaluated*/
/* */
/*  Results:*/
/*  Returns 2 * u * f( a + u**2 )*/
/* */
/*  Side effects:*/
/*  Whatever f does.*/
/* */
/* ----------------------------------------------------------------------*/

ret  ::math::calculus::u_sqrtSingLower ( type f , type a , type u ) {
    set cmd $f
    lappend cmd [expr { $a + $u * $u }]
    set y [eval $cmd]
    return [expr { 2. * $u * $y }]
}

/* ----------------------------------------------------------------------*/
/* */
/*  u_sqrtSingUpper --*/
/*  Change of variable for integrating over an interval with*/
/*  an inverse square root singularity at the upper bound.*/
/* */
/*  Parameters:*/
/*  f - Function being integrated*/
/*  b - Upper bound*/
/*  u - sqrt(b-x), where x is the abscissa where f is to be evaluated*/
/* */
/*  Results:*/
/*  Returns 2 * u * f( b - u**2 )*/
/* */
/*  Side effects:*/
/*  Whatever f does.*/
/* */
/* ----------------------------------------------------------------------*/

ret  ::math::calculus::u_sqrtSingUpper ( type f , type b , type u ) {
    set cmd $f
    lappend cmd [expr { $b - $u * $u }]
    set y [eval $cmd]
    return [expr { 2. * $u * $y }]
}

/* ----------------------------------------------------------------------*/
/* */
/*  math::calculus::romberg_sqrtSingLower --*/
/*  Integrate a function with an inverse square root singularity*/
/*  at the lower bound*/
/* */
/*  Usage:*/
/*  Same as 'romberg'*/
/* */
/*  The 'romberg_sqrtSingLower' procedure is a wrapper for 'romberg'*/
/*  for integrating a function with an inverse square root singularity*/
/*  at the lower bound of the interval.  It works by making the change*/
/*  of variable u = sqrt( x-a ).*/
/* */
/* ----------------------------------------------------------------------*/

ret  ::math::calculus::romberg_sqrtSingLower ( type f , type a , type b , type args ) {
    if { ![string is double -strict $a] } {
    return -code error [expectDouble $a]
    }
    if { ![string is double -strict $b] } {
    return -code error [expectDouble $b]
    }
    if { $a >= $b } {
    return -code error "limits of integration out of order"
    }
    set f [lreplace $f 0 0 [uplevel 1 [list namespace which [lindex $f 0]]]]
    set f [list u_sqrtSingLower $f $a]
    return [eval [linsert $args 0 \
                 romberg $f 0 [expr { sqrt( $b - $a ) }]]]
}

/* ----------------------------------------------------------------------*/
/* */
/*  math::calculus::romberg_sqrtSingUpper --*/
/*  Integrate a function with an inverse square root singularity*/
/*  at the upper bound*/
/* */
/*  Usage:*/
/*  Same as 'romberg'*/
/* */
/*  The 'romberg_sqrtSingUpper' procedure is a wrapper for 'romberg'*/
/*  for integrating a function with an inverse square root singularity*/
/*  at the upper bound of the interval.  It works by making the change*/
/*  of variable u = sqrt( b-x ).*/
/* */
/* ----------------------------------------------------------------------*/

ret  ::math::calculus::romberg_sqrtSingUpper ( type f , type a , type b , type args ) {
    if { ![string is double -strict $a] } {
    return -code error [expectDouble $a]
    }
    if { ![string is double -strict $b] } {
    return -code error [expectDouble $b]
    }
    if { $a >= $b } {
    return -code error "limits of integration out of order"
    }
    set f [lreplace $f 0 0 [uplevel 1 [list namespace which [lindex $f 0]]]]
    set f [list u_sqrtSingUpper $f $b]
    return [eval [linsert $args 0 \
              romberg $f 0. [expr { sqrt( $b - $a ) }]]]
}

/* ----------------------------------------------------------------------*/
/* */
/*  u_powerLawLower --*/
/*  Change of variable for integrating over an interval with*/
/*  an integrable power law singularity at the lower bound.*/
/* */
/*  Parameters:*/
/*  f - Function being integrated*/
/*  gammaover1mgamma - gamma / (1 - gamma), where gamma is the power*/
/*  oneover1mgamma - 1 / (1 - gamma), where gamma is the power*/
/*  a - Lower limit of integration*/
/*  u - Changed variable u == (x-a)**(1-gamma)*/
/* */
/*  Results:*/
/*  Returns u**(1/1-gamma) * f(a + u**(1/1-gamma) ).*/
/* */
/*  Side effects:*/
/*  Whatever f does.*/
/* */
/* ----------------------------------------------------------------------*/

ret  ::math::calculus::u_powerLawLower ( type f , type gammaover1mgamma , type oneover1mgamma
                     , type a , type u ) {
    set cmd $f
    lappend cmd [expr { $a + pow( $u, $oneover1mgamma ) }]
    set y [eval $cmd]
    return [expr { $y * pow( $u, $gammaover1mgamma ) }]
}

/* ----------------------------------------------------------------------*/
/* */
/*  math::calculus::romberg_powerLawLower --*/
/*  Integrate a function with an integrable power law singularity*/
/*  at the lower bound*/
/* */
/*  Usage:*/
/*  romberg_powerLawLower gamma f a b ?-option value...?*/
/* */
/*  Parameters:*/
/*  gamma - Power (0<gamma<1) of the singularity*/
/*  f - Function to integrate.  Must be a single Tcl command,*/
/*      to which will be appended the abscissa at which the function*/
/*      should be evaluated.  f is expected to have an integrable*/
/*      power law singularity at the lower endpoint; that is, the*/
/*      integrand is expected to diverge as (x-a)**gamma.*/
/*  a - One bound of the interval*/
/*  b - The other bound of the interval.  a and b need not be in*/
/*      ascending order.*/
/* */
/*  Options:*/
/*  -abserror ABSERROR*/
/*      Requests that the integration be performed to make*/
/*      the estimated absolute error of the integral less than*/
/*      the given value.  Default is 1.e-10.*/
/*  -relerror RELERROR*/
/*      Requests that the integration be performed to make*/
/*      the estimated absolute error of the integral less than*/
/*      the given value.  Default is 1.e-6.*/
/*  -degree N*/
/*      Specifies the degree of the polynomial that will be*/
/*      used to extrapolate to a zero step size.  -degree 0*/
/*      requests integration with the midpoint rule; -degree 1*/
/*      is equivalent to Simpson's 3/8 rule; higher degrees*/
/*      are difficult to describe but (within reason) give*/
/*      faster convergence for smooth functions.  Default is*/
/*      -degree 4.*/
/*  -maxiter N*/
/*      Specifies the maximum number of triplings of the*/
/*      number of steps to take in integration.  At most*/
/*      3**N function evaluations will be performed in*/
/*      integrating with -maxiter N.  The integration*/
/*      will terminate at that time, even if the result*/
/*      satisfies neither the -relerror nor -abserror tests.*/
/* */
/*  Results:*/
/*  Returns a two-element list.  The first element is the estimated*/
/*  value of the integral; the second is the estimated absolute*/
/*  error of the value.*/
/* */
/*  The 'romberg_sqrtSingLower' procedure is a wrapper for 'romberg'*/
/*  for integrating a function with an integrable power law singularity*/
/*  at the lower bound of the interval.  It works by making the change*/
/*  of variable u = (x-a)**(1-gamma).*/
/* */
/* ----------------------------------------------------------------------*/

ret  ::math::calculus::romberg_powerLawLower ( type gamma , type f , type a , type b , type args ) {
    if { ![string is double -strict $gamma] } {
    return -code error [expectDouble $gamma]
    }
    if { $gamma <= 0.0 || $gamma >= 1.0 } {
    return -code error -errorcode [list romberg gammaTooBig] \
        "gamma must lie in the interval (0,1)"
    }
    if { ![string is double -strict $a] } {
    return -code error [expectDouble $a]
    }
    if { ![string is double -strict $b] } {
    return -code error [expectDouble $b]
    }
    if { $a >= $b } {
    return -code error "limits of integration out of order"
    }
    set f [lreplace $f 0 0 [uplevel 1 [list namespace which [lindex $f 0]]]]
    set onemgamma [expr { 1. - $gamma }]
    set f [list u_powerLawLower $f \
           [expr { $gamma / $onemgamma }] \
           [expr { 1 / $onemgamma }] \
           $a]
    
    set limit [expr { pow( $b - $a, $onemgamma ) }]
    set result {}
    foreach v [eval [linsert $args 0 romberg $f 0 $limit]] {
    lappend result [expr { $v / $onemgamma }]
    }
    return $result

}

/* ----------------------------------------------------------------------*/
/* */
/*  u_powerLawLower --*/
/*  Change of variable for integrating over an interval with*/
/*  an integrable power law singularity at the upper bound.*/
/* */
/*  Parameters:*/
/*  f - Function being integrated*/
/*  gammaover1mgamma - gamma / (1 - gamma), where gamma is the power*/
/*  oneover1mgamma - 1 / (1 - gamma), where gamma is the power*/
/*  b - Upper limit of integration*/
/*  u - Changed variable u == (b-x)**(1-gamma)*/
/* */
/*  Results:*/
/*  Returns u**(1/1-gamma) * f(b-u**(1/1-gamma) ).*/
/* */
/*  Side effects:*/
/*  Whatever f does.*/
/* */
/* ----------------------------------------------------------------------*/

ret  ::math::calculus::u_powerLawUpper ( type f , type gammaover1mgamma , type oneover1mgamma
                     , type b , type u ) {
    set cmd $f
    lappend cmd [expr { $b - pow( $u, $oneover1mgamma ) }]
    set y [eval $cmd]
    return [expr { $y * pow( $u, $gammaover1mgamma ) }]
}

/* ----------------------------------------------------------------------*/
/* */
/*  math::calculus::romberg_powerLawUpper --*/
/*  Integrate a function with an integrable power law singularity*/
/*  at the upper bound*/
/* */
/*  Usage:*/
/*  romberg_powerLawLower gamma f a b ?-option value...?*/
/* */
/*  Parameters:*/
/*  gamma - Power (0<gamma<1) of the singularity*/
/*  f - Function to integrate.  Must be a single Tcl command,*/
/*      to which will be appended the abscissa at which the function*/
/*      should be evaluated.  f is expected to have an integrable*/
/*      power law singularity at the upper endpoint; that is, the*/
/*      integrand is expected to diverge as (b-x)**gamma.*/
/*  a - One bound of the interval*/
/*  b - The other bound of the interval.  a and b need not be in*/
/*      ascending order.*/
/* */
/*  Options:*/
/*  -abserror ABSERROR*/
/*      Requests that the integration be performed to make*/
/*      the estimated absolute error of the integral less than*/
/*      the given value.  Default is 1.e-10.*/
/*  -relerror RELERROR*/
/*      Requests that the integration be performed to make*/
/*      the estimated absolute error of the integral less than*/
/*      the given value.  Default is 1.e-6.*/
/*  -degree N*/
/*      Specifies the degree of the polynomial that will be*/
/*      used to extrapolate to a zero step size.  -degree 0*/
/*      requests integration with the midpoint rule; -degree 1*/
/*      is equivalent to Simpson's 3/8 rule; higher degrees*/
/*      are difficult to describe but (within reason) give*/
/*      faster convergence for smooth functions.  Default is*/
/*      -degree 4.*/
/*  -maxiter N*/
/*      Specifies the maximum number of triplings of the*/
/*      number of steps to take in integration.  At most*/
/*      3**N function evaluations will be performed in*/
/*      integrating with -maxiter N.  The integration*/
/*      will terminate at that time, even if the result*/
/*      satisfies neither the -relerror nor -abserror tests.*/
/* */
/*  Results:*/
/*  Returns a two-element list.  The first element is the estimated*/
/*  value of the integral; the second is the estimated absolute*/
/*  error of the value.*/
/* */
/*  The 'romberg_PowerLawUpper' procedure is a wrapper for 'romberg'*/
/*  for integrating a function with an integrable power law singularity*/
/*  at the upper bound of the interval.  It works by making the change*/
/*  of variable u = (b-x)**(1-gamma).*/
/* */
/* ----------------------------------------------------------------------*/

ret  ::math::calculus::romberg_powerLawUpper ( type gamma , type f , type a , type b , type args ) {
    if { ![string is double -strict $gamma] } {
    return -code error [expectDouble $gamma]
    }
    if { $gamma <= 0.0 || $gamma >= 1.0 } {
    return -code error -errorcode [list romberg gammaTooBig] \
        "gamma must lie in the interval (0,1)"
    }
    if { ![string is double -strict $a] } {
    return -code error [expectDouble $a]
    }
    if { ![string is double -strict $b] } {
    return -code error [expectDouble $b]
    }
    if { $a >= $b } {
    return -code error "limits of integration out of order"
    }
    set f [lreplace $f 0 0 [uplevel 1 [list namespace which [lindex $f 0]]]]
    set onemgamma [expr { 1. - $gamma }]
    set f [list u_powerLawUpper $f \
           [expr { $gamma / $onemgamma }] \
           [expr { 1. / $onemgamma }] \
           $b]
    
    set limit [expr { pow( $b - $a, $onemgamma ) }]
    set result {}
    foreach v [eval [linsert $args 0 romberg $f 0 $limit]] {
    lappend result [expr { $v / $onemgamma }]
    }
    return $result

}

/* ----------------------------------------------------------------------*/
/* */
/*  u_expUpper --*/
/* */
/*  Change of variable to integrate a function that decays*/
/*  exponentially.*/
/* */
/*  Parameters:*/
/*  f - Function to integrate*/
/*  u - Changed variable u = exp(-x)*/
/* */
/*  Results:*/
/*  Returns (1/u)*f(-log(u))*/
/* */
/*  Side effects:*/
/*  Whatever f does.*/
/* */
/* ----------------------------------------------------------------------*/

ret  ::math::calculus::u_expUpper ( type f , type u ) {
    set cmd $f
    lappend cmd [expr { -log($u) }]
    set y [eval $cmd]
    return [expr { $y / $u }]
}

/* ----------------------------------------------------------------------*/
/* */
/*  romberg_expUpper --*/
/* */
/*  Integrate a function that decays exponentially over a*/
/*  half-infinite interval.*/
/* */
/*  Parameters:*/
/*  Same as romberg.  The upper limit of integration, 'b',*/
/*  is expected to be very large.*/
/* */
/*  Results:*/
/*  Same as romberg.*/
/* */
/*  The romberg_expUpper function operates by making the change of*/
/*  variable, u = exp(-x).*/
/* */
/* ----------------------------------------------------------------------*/

ret  ::math::calculus::romberg_expUpper ( type f , type a , type b , type args ) {
    if { ![string is double -strict $a] } {
    return -code error [expectDouble $a]
    }
    if { ![string is double -strict $b] } {
    return -code error [expectDouble $b]
    }
    if { $a >= $b } {
    return -code error "limits of integration out of order"
    }
    set f [lreplace $f 0 0 [uplevel 1 [list namespace which [lindex $f 0]]]]
    set f [list u_expUpper $f]
    return [eval [linsert $args 0 \
              romberg $f [expr {exp(-$b)}] [expr {exp(-$a)}]]]
}

/* ----------------------------------------------------------------------*/
/* */
/*  u_expLower --*/
/* */
/*  Change of variable to integrate a function that grows*/
/*  exponentially.*/
/* */
/*  Parameters:*/
/*  f - Function to integrate*/
/*  u - Changed variable u = exp(x)*/
/* */
/*  Results:*/
/*  Returns (1/u)*f(log(u))*/
/* */
/*  Side effects:*/
/*  Whatever f does.*/
/* */
/* ----------------------------------------------------------------------*/

ret  ::math::calculus::u_expLower ( type f , type u ) {
    set cmd $f
    lappend cmd [expr { log($u) }]
    set y [eval $cmd]
    return [expr { $y / $u }]
}

/* ----------------------------------------------------------------------*/
/* */
/*  romberg_expLower --*/
/* */
/*  Integrate a function that grows exponentially over a*/
/*  half-infinite interval.*/
/* */
/*  Parameters:*/
/*  Same as romberg.  The lower limit of integration, 'a',*/
/*  is expected to be very large and negative.*/
/* */
/*  Results:*/
/*  Same as romberg.*/
/* */
/*  The romberg_expUpper function operates by making the change of*/
/*  variable, u = exp(x).*/
/* */
/* ----------------------------------------------------------------------*/

ret  ::math::calculus::romberg_expLower ( type f , type a , type b , type args ) {
    if { ![string is double -strict $a] } {
    return -code error [expectDouble $a]
    }
    if { ![string is double -strict $b] } {
    return -code error [expectDouble $b]
    }
    if { $a >= $b } {
    return -code error "limits of integration out of order"
    }
    set f [lreplace $f 0 0 [uplevel 1 [list namespace which [lindex $f 0]]]]
    set f [list u_expLower $f]
    return [eval [linsert $args 0 \
              romberg $f [expr {exp($a)}] [expr {exp($b)}]]]
}


/*  regula_falsi --*/
/*     Compute the zero of a function via regula falsi*/
/*  Arguments:*/
/*     f       Name of the procedure/command that evaluates the function*/
/*     xb      Start of the interval that brackets the zero*/
/*     xe      End of the interval that brackets the zero*/
/*     eps     Relative error that is allowed (default: 1.0e-4)*/
/*  Result:*/
/*     Estimate of the zero, such that the estimated (!)*/
/*     error < eps * abs(xe-xb)*/
/*  Note:*/
/*     f(xb)*f(xe) must be negative and eps must be positive*/
/* */
ret  ::math::calculus::regula_falsi ( type f , type xb , type xe , optional eps =1.0e-4 ) {
    if { $eps <= 0.0 } {
       return -code error "Relative error must be positive"
    }

    set fb [$f $xb]
    set fe [$f $xe]

    if { $fb * $fe > 0.0 } {
       return -code error "Interval must be chosen such that the \
function has a different sign at the beginning than at the end"
    }

    set max_error [expr {$eps * abs($xe-$xb)}]
    set interval  [expr {abs($xe-$xb)}]

    while { $interval > $max_error } {
       set coeff [expr {($fe-$fb)/($xe-$xb)}]
       set xi    [expr {$xb-$fb/$coeff}]
       set fi    [$f $xi]

       if { $fi == 0.0 } {
          break
       }
       set diff1 [expr {abs($xe-$xi)}]
       set diff2 [expr {abs($xb-$xi)}]
       if { $diff1 > $diff2 } {
          set interval $diff2
       } else {
          set interval $diff1
       }

       if { $fb*$fi < 0.0 } {
          set xe $xi
          set fe $fi
       } else {
          set xb $xi
          set fb $fi
       }
    }

    return $xi
}