Index: modules/math/calculus.CHANGES
===================================================================
RCS file: /cvsroot/tcllib/tcllib/modules/math/calculus.CHANGES,v
retrieving revision 1.2
diff -u -r1.2 calculus.CHANGES
--- modules/math/calculus.CHANGES	15 Mar 2002 21:42:44 -0000	1.2
+++ modules/math/calculus.CHANGES	31 May 2002 20:59:52 -0000
@@ -10,6 +10,12 @@
 Version 0.2: november 2001
    Extended with Euler and Heun methods, 2D and 3D simple integration
 
-Version 0.3: november 2001
-   Implementation of Runge-Kutta method. Converted the documentation to
-   doctools format
+Version 0.3: march 2002
+   Implemented Runge-Kutta, converted documentation to doctools' 
+   man format
+
+Version 0.4: march 2002
+   Implemented Newton-Raphson method for finding roots of equations
+
+Version 0.5: may 2002
+   Fixed problem with namespaces
Index: modules/math/calculus.README
===================================================================
RCS file: /cvsroot/tcllib/tcllib/modules/math/calculus.README,v
retrieving revision 1.2
diff -u -r1.2 calculus.README
--- modules/math/calculus.README	15 Mar 2002 21:42:44 -0000	1.2
+++ modules/math/calculus.README	31 May 2002 20:59:52 -0000
@@ -1,6 +1,6 @@
-Package: Calculus
-----------------
-The Calculus package is an all-Tcl package that implements
+Package: math::calculus
+-----------------------
+The math::calculus package is an all-Tcl package that implements 
 several basic numerical algorithms, such as the integration
 of functions of one variable or the integration of ordinary
 differential equations.
@@ -12,7 +12,10 @@
 calculus.test  - Several simple tests
 calculus.html  - Documentation of the package
 
-The current version is: 0.3, the first public release, march 2002
+The current version is: 0.5, may 2002
+
+This package is available as part of Tcllib at:
+   http://tcllib.sourceforge.net
 
 Please contact Arjen Markus (arjen.markus@wldelft.nl) for questions,
 bug reports, enhancements and so on.
Index: modules/math/calculus.doc
===================================================================
RCS file: /cvsroot/tcllib/tcllib/modules/math/calculus.doc,v
retrieving revision 1.1
diff -u -r1.1 calculus.doc
--- modules/math/calculus.doc	18 Jan 2002 20:28:22 -0000	1.1
+++ modules/math/calculus.doc	31 May 2002 20:59:52 -0000
@@ -1,4 +1,18 @@
 [pageheader "Package: Calculus"]
+[synopsis \
+{package require Tcl 8.2
+package require math::calculus 0.5
+::math::calculus::integral begin end nosteps func
+::math::calculus::integralExpr begin end nosteps expression
+::math::calculus::integral2D xinterval yinterval func
+::math::calculus::integral3D xinterval yinterval zinterval func
+::math::calculus::eulerStep t tstep xvec func
+::math::calculus::heunStep t tstep xvec func
+::math::calculus::rungeKuttaStep tstep xvec func
+::math::calculus::boundaryValueSecondOrder coeff_func force_func leftbnd rightbnd nostep}]
+::math::calculus::newtonRaphson func deriv initval
+::math::calculus::newtonRaphsonParameters maxiter tolerance
+
 [section "Introduction"]
 The package Calculus implements several simple mathematical algorithms,
 such as the integration of a function over an interval and the numerical
@@ -11,7 +25,7 @@
 This document describes the procedures and explains their usage.
 
 [section "Version and copyright"]
-This document describes [italic Calculus], version 0.3, december 2001.
+This document describes [italic ::math::calculus], version 0.5, may 2002.
 [par]
 Usage of Calculus is free, as long as you acknowledge the
 author, Arjen Markus (e-mail: arjen.markus@wldelft.nl).
@@ -21,7 +35,7 @@
 [section "Procedures"]
 The Calculus package defines the following public procedures:
 [ulist]
-[item][italic "Integral begin end nosteps func"]
+[item][italic "integral begin end nosteps func"]
       [break]
       Determine the integral of the given function using the Simpson
       rule. The interval for the integration is [lb]begin,end[rb].
@@ -34,7 +48,7 @@
       single argument.
       [par]
 
-[item][italic "IntegralExpr begin end nosteps expression"]
+[item][italic "integralExpr begin end nosteps expression"]
       [break]
       Similar to the previous proc, this one determines the integral of
       the given [italic expression] using the Simpson rule.
@@ -48,9 +62,9 @@
       use the variable "x" as the only variable (the "integrate")
       [par]
 
-[item][italic "Integral2D xinterval yinterval func"]
+[item][italic "integral2D xinterval yinterval func"]
       [break]
-      The [italic Integral2D] procedure calculates the integral of
+      The [italic integral2D] procedure calculates the integral of
       a function of two variables over the rectangle given by the
       first two arguments, each a list of three items, the start and
       stop interval for the variable and the number of steps.
@@ -64,15 +78,15 @@
       value.
       [par]
 
-[item][italic "Integral3D xinterval yinterval zinterval func"]
+[item][italic "integral3D xinterval yinterval zinterval func"]
       [break]
-      The [italic Integral3D] procedure is the three-dimensional
-      equivalent of [italic Intergral2D]. The function taking three
+      The [italic integral3D] procedure is the three-dimensional
+      equivalent of [italic intergral2D]. The function taking three
       arguments is integrated over the block in 3D space given by the
       intervals.
       [par]
 
-[item][italic "EulerStep t tstep xvec func"]
+[item][italic "eulerStep t tstep xvec func"]
       [break]
       Set a single step in the numerical integration of a system of
       differential equations. The method used is Euler's.
@@ -89,7 +103,7 @@
       xvec and the return value of "func" must match).
       [par]
 
-[item][italic "HeunStep t tstep xvec func"]
+[item][italic "heunStep t tstep xvec func"]
       [break]
       Set a single step in the numerical integration of a system of
       differential equations. The method used is Heun's.
@@ -106,7 +120,7 @@
       xvec and the return value of "func" must match).
       [par]
 
-[item][italic "RungeKuttaStep tstep xvec func"]
+[item][italic "rungeKuttaStep tstep xvec func"]
       [break]
       Set a single step in the numerical integration of a system of
       differential equations. The method used is Runge-Kutta 4th
@@ -124,7 +138,7 @@
       xvec and the return value of "func" must match).
       [par]
 
-[item][italic "BoundaryValueSecondOrder coeff_func force_func leftbnd rightbnd nostep"]
+[item][italic "boundaryValueSecondOrder coeff_func force_func leftbnd rightbnd nostep"]
       [break]
       Solve a second order linear differential equation with boundary
       values at two sides. The equation has to be of the form:
@@ -164,7 +178,7 @@
       values of the solution.
       [par]
 
-[item][italic "SolveTriDiagonal acoeff bcoeff ccoeff dvalue"
+[item][italic "solveTriDiagonal acoeff bcoeff ccoeff dvalue"]
       [break]
       Solve a system of linear equations Ax = b with A a tridiagonal
       matrix. Returns the solution as a list.
@@ -173,13 +187,60 @@
       [italic bcoeff] - List of values on the main diagonal
       [italic ccoeff] - List of values on the upper diagonal
       [italic dvalue] - List of values on the righthand-side
+      [par]
+
+[item][italic "newtonRaphson func deriv initval"]
+      [break]
+      Determine the root of an equation given by [italic "f(x) = 0"],
+      using the Newton-Raphson method.
+      [break]
+      [italic func]    - Name of the procedure that calculates the function value
+      [italic deriv    - Name of the procedure that calculates the derivative of the function
+      [italic initval] - Initial value for the iteration
+      [par]
+
+[item][italic "newtonRaphsonParameters maxiter tolerance"]
+      [break]
+      Set new values for the two parameters that gouvern the Newton-Raphson method.
+      [break]
+      [italic maxiter]   - Maximum number of iterations
+      [italic tolerance] - Relative error in the calculation
+      [par]
 [endlist]
 
-[italic Note:]
-[break]
-The integration of the differential equations via Runge-Kutta
-is yet to be implemented.
+[italic Notes:]
 [break]
+Several of the above procedures take the [italic names] of procedures as
+arguments. To avoid problems with the [italic visibility] of these
+procedures, the fully-qualified name of these procedures is determined
+inside the calculus routines. For the user this has only one
+consequence: the named procedure must be visible in the calling
+procedure. For instance:
+
+[preserve]
+    namespace eval ::mySpace {
+       namespace export calcfunc
+       proc calcfunc { x } { return $x }
+    }
+    #
+    # Use a fully-qualified name
+    #
+    namespace eval ::myCalc {
+       proc detIntegral { begin end } {
+          return [lb]integral $begin $end 100 ::mySpace::calcfunc[rb]
+       }
+    }
+    #
+    # Import the name
+    #
+    namespace eval ::myCalc {
+       namespace import ::mySpace::calcfunc
+       proc detIntegral { begin end } {
+          return [lb]integral $begin $end 100 calcfunc[rb]
+       }
+    }
+[endpreserve]
+[par]
 Enhancements for the second-order boundary value problem:
 [ulist]
 [item]Other types of boundary conditions (zero gradient, zero flux)
@@ -193,11 +254,11 @@
 Integrate x over the interval [lb]0,100[rb] (20 steps):
 [preserve]
 proc linear_func { x } { return $x }
-puts "Integral: ::Calculus::Integral 0 100 20 linear_func"
+puts "Integral: [lb]::math::calculus::Integral 0 100 20 linear_func[rb]"
 [endpreserve]
 For simple functions, the alternative could be:
 [preserve]
-puts "Integral: ::Calculus::IntegralExpr 0 100 20 {$x}"
+puts "Integral: [lb]::math::calculus::IntegralExpr 0 100 20 {$x}[rb]"
 [endpreserve]
 Do not forget the braces!
 [par]
@@ -213,19 +274,19 @@
 Then this system can be solved with code like this:
 [preserve]
 proc dampened_oscillator { t xvec } {
-   set x  [lb]lindex $xvec 0[rb]
-   set x1 [lb]lindex $xvec 1[rb]
-   return [lb]list $x1 [lb]expr {-$x1-$x}[rb][rb]
+   set x  [lb]lindex \$xvec 0[rb]
+   set x1 [lb]lindex \$xvec 1[rb]
+   return [lb]list \$x1 [lb]expr {-\$x1-\$x}[rb][rb]
 }
 
 set xvec   { 1.0 0.0 }
 set t      0.0
 set tstep  0.1
-for { set i 0 } { $i < 20 } { incr i } {
-   set result [lb]::Calculus::EulerStep $t $tstep $xvec dampened_oscillator[rb]
-   puts "Result ($t): $result"
-   set t      [lb]expr {$t+$tstep}[rb]
-   set xvec   $result
+for { set i 0 } { \$i < 20 } { incr i } {
+   set result [lb]::math::calculus::eulerStep \$t \$tstep \$xvec dampened_oscillator[rb]
+   puts "Result (\$t): \$result"
+   set t      [lb]expr {\$t+\$tstep}[rb]
+   set xvec   \$result
 }
 [endpreserve]
 Suppose we have the boundary value problem:
@@ -239,12 +300,12 @@
 [par]
 It can be solved with the following fragment:
 [preserve]
-   proc coeffs { x } { return [list $::Diff 0.0 $::decay] }
+   proc coeffs { x } { return [lb]list \$::Diff 0.0 \$::decay[rb] }
    proc force  { x } { return 0.0 }
 
    set Diff   1.0e-2
    set decay  0.0001
    set length 100.0
-   set y [lb]::Calculus::BoundaryValueSecondOrder coeffs force {0.0 1.0} \
-      [list $length 0.0] 100
+   set y [lb]::math::calculus::boundaryValueSecondOrder coeffs force {0.0 1.0} \
+      [lb]list \$length 0.0[rb] 100[rb]
 [endpreserve]
Index: modules/math/calculus.man
===================================================================
RCS file: /cvsroot/tcllib/tcllib/modules/math/calculus.man,v
retrieving revision 1.2
diff -u -r1.2 calculus.man
--- modules/math/calculus.man	26 Mar 2002 07:40:36 -0000	1.2
+++ modules/math/calculus.man	31 May 2002 20:59:52 -0000
@@ -3,7 +3,6 @@
 [titledesc {Integration and ordinary differential equations}]
 [description]
 [para]
-
 This package implements several simple mathematical algorithms:
 
 [list_begin bullet]
@@ -14,22 +13,20 @@
 The numerical integration of a system of ordinary differential
 equations.
 
+[bullet]
+Estimating the root(s) of an equation of one variable.
+
 [list_end]
 
 [para]
-
 The package is fully implemented in Tcl. No particular attention has
 been paid to the accuracy of the calculations. Instead, well-known
 algorithms have been used in a straightforward manner.
-
 [para]
-
 This document describes the procedures and explains their usage.
 
 [section "PROCEDURES"]
-
-It defines the following public procedures:
-
+This package defines the following public procedures:
 [list_begin definitions]
 
 [call [cmd ::math::calculus::integral] [arg begin] [arg end] [arg nosteps] [arg func]]
@@ -42,6 +39,7 @@
 [nl]
 [arg func] - Function to be integrated. It should take one
 single argument.
+[para]
 
 [call [cmd ::math::calculus::integralExpr] [arg begin] [arg end] [arg nosteps] [arg expression]]
 Similar to the previous proc, this one determines the integral of
@@ -54,6 +52,7 @@
 [nl]
 [arg expression] - Expression to be integrated. It should
 use the variable "x" as the only variable (the "integrate")
+[para]
 
 [call [cmd ::math::calculus::integral2D] [arg xinterval] [arg yinterval] [arg func]]
 The [strong Integral2D] procedure calculates the integral of
@@ -68,12 +67,14 @@
 [nl]
 The function must take two arguments and return the function
 value.
+[para]
 
 [call [cmd ::math::calculus::integral3D] [arg xinterval] [arg yinterval] [arg zinterval] [arg func]]
 The [strong Integral3D] procedure is the three-dimensional
 equivalent of [strong Intergral2D]. The function taking three
 arguments is integrated over the block in 3D space given by the
 intervals.
+[para]
 
 [call [cmd ::math::calculus::eulerStep] [arg t] [arg tstep] [arg xvec] [arg func]]
 Set a single step in the numerical integration of a system of
@@ -89,6 +90,7 @@
 [arg func] - Function of t and the dependent values, returning
 a list of the derivatives of the dependent values. (The lengths of
 xvec and the return value of "func" must match).
+[para]
 
 [call [cmd ::math::calculus::heunStep] [arg t] [arg tstep] [arg xvec] [arg func]]
 Set a single step in the numerical integration of a system of
@@ -104,6 +106,7 @@
 [arg func] - Function of t and the dependent values, returning
 a list of the derivatives of the dependent values. (The lengths of
 xvec and the return value of "func" must match).
+[para]
 
 [call [cmd ::math::calculus::rungeKuttaStep] [arg tstep] [arg xvec] [arg func]]
 Set a single step in the numerical integration of a system of
@@ -120,6 +123,7 @@
 [arg func] - Function of t and the dependent values, returning
 a list of the derivatives of the dependent values. (The lengths of
 xvec and the return value of "func" must match).
+[para]
 
 [call [cmd ::math::calculus::boundaryValueSecondOrder] [arg coeff_func] [arg force_func] [arg leftbnd] [arg rightbnd] [arg nostep]]
 Solve a second order linear differential equation with boundary
@@ -164,7 +168,7 @@
 [nl]
 The procedure returns a list of x-coordinates and the approximated
 values of the solution.
-
+[para]
 
 [call [cmd ::math::calculus::solveTriDiagonal] [arg acoeff] [arg bcoeff] [arg ccoeff] [arg dvalue]]
 Solve a system of linear equations Ax = b with A a tridiagonal
@@ -177,15 +181,63 @@
 [arg ccoeff] - List of values on the upper diagonal
 [nl]
 [arg dvalue] - List of values on the righthand-side
-[nl]
-[list_end]
+[para]
 
+[call [cmd ::math::calculus::newtonRaphson] [arg func] [arg deriv] [arg initval]]
+Determine the root of an equation given by
+[example_begin]
+    func(x) = 0
+[example_end]
+using the method of Newton-Raphson. The procedure takes the following
+arguments:
+[nl]
+[arg func] - Procedure that returns the value the function at x
+[nl]
+[arg deriv] - Procedure that returns the derivative of the function at x
+[nl]
+[arg initval] - Initial value for x
 [para]
 
-[strong Note:]
+[call [cmd ::math::calculus::newtonRaphsonParameters] [arg maxiter] [arg tolerance]]
+Set the numerical parameters for the Newton-Raphson method:
+[nl]
+[arg maxiter] - Maximum number of iteration steps (defaults to 20)
+[nl]
+[arg tolerance] - Relative precision (defaults to 0.001)
+[list_end]
 [para]
-The integration of the differential equations via Runge-Kutta
-is yet to be implemented.
+
+[strong Notes:]
+[nl]
+Several of the above procedures take the [strong names] of procedures as
+arguments. To avoid problems with the [strong visibility] of these
+procedures, the fully-qualified name of these procedures is determined
+inside the calculus routines. For the user this has only one
+consequence: the named procedure must be visible in the calling
+procedure. For instance:
+[example_begin]
+    namespace eval ::mySpace {
+       namespace export calcfunc
+       proc calcfunc { x } { return $x }
+    }
+    #
+    # Use a fully-qualified name
+    #
+    namespace eval ::myCalc {
+       proc detIntegral { begin end } {
+          return [lb]integral $begin $end 100 ::mySpace::calcfunc[rb]
+       }
+    }
+    #
+    # Import the name
+    #
+    namespace eval ::myCalc {
+       namespace import ::mySpace::calcfunc
+       proc detIntegral { begin end } {
+          return [lb]integral $begin $end 100 calcfunc[rb]
+       }
+    }
+[example_end]
 [para]
 Enhancements for the second-order boundary value problem:
 [list_begin bullet]
@@ -202,16 +254,16 @@
 Integrate x over the interval [lb]0,100[rb] (20 steps):
 [example_begin]
 proc linear_func { x } { return $x }
-puts "Integral: [lb]::math::calculus::integral 0 100 20 linear_func[rb]"
+puts "Integral: [lb]::math::calculus::Integral 0 100 20 linear_func[rb]"
 [example_end]
 For simple functions, the alternative could be:
 [example_begin]
-puts "Integral: [lb]::math::calculus::integralExpr 0 100 20 {$x}[rb]"
+puts "Integral: [lb]::math::calculus::IntegralExpr 0 100 20 {$x}[rb]"
 [example_end]
 Do not forget the braces!
 [para]
 The differential equation for a dampened oscillator:
-[para]
+[nl]
 [example_begin]
 x'' + rx' + wx = 0
 [example_end]
@@ -256,9 +308,10 @@
    set decay  0.0001
    set length 100.0
 
-   set y [lb]::math::calculus::boundaryValueSecondOrder coeffs force {0.0 1.0} \
-      [lb]list $length 0.0[rb] 100[rb]
+   set y [lb]::math::calculus::boundaryValueSecondOrder \
+      coeffs force {0.0 1.0} [lb]list $length 0.0[rb] 100[rb]
 [example_end]
 
-[keywords math calculus integration "differential equations"]
+[keywords math calculus integration "differential equations" roots]
+
 [manpage_end]
Index: modules/math/calculus.tcl
===================================================================
RCS file: /cvsroot/tcllib/tcllib/modules/math/calculus.tcl,v
retrieving revision 1.2
diff -u -r1.2 calculus.tcl
--- modules/math/calculus.tcl	15 Mar 2002 21:42:44 -0000	1.2
+++ modules/math/calculus.tcl	31 May 2002 20:59:52 -0000
@@ -9,14 +9,19 @@
 # math::calculus --
 #    Namespace for the commands
 #
-
 namespace eval ::math::calculus {
-   namespace export integral integralExpr integral2D integral3D
-   namespace export eulerStep heunStep rungeKuttaStep
-   namespace export boundaryValueSecondOrder solveTriDiagonal
+   namespace export \
+          integral integralExpr integral2D integral3D \
+          eulerStep heunStep rungeKuttaStep           \
+          boundaryValueSecondOrder solveTriDiagonal   \
+          newtonRaphson newtonRaphsonParameters
+
+   variable nr_maxiter    20
+   variable nr_tolerance   0.001
+
 }
 
-# Integral --
+# integral --
 #    Integrate a function over a given interval using the Simpson rule
 #
 # Arguments:
@@ -34,11 +39,11 @@
    set hdelta   [expr {$delta/2.0}]
    set result   0.0
    set xval     $begin
-   set func_end [$func $xval]
+   set func_end [uplevel 1 [list $func $xval]]
    for { set i 1 } { $i <= $nosteps } { incr i } {
       set func_begin  $func_end
-      set func_middle [$func [expr {$xval+$hdelta}]]
-      set func_end    [$func [expr {$xval+$delta}]]
+      set func_middle [uplevel 1 [list $func [expr {$xval+$hdelta}]]]
+      set func_end    [uplevel 1 [list $func [expr {$xval+$delta}]]]
       set result      [expr  {$result+$func_begin+4.0*$func_middle+$func_end}]
 
       set xval        [expr {$begin+double($i)*$delta}]
@@ -47,7 +52,7 @@
    return [expr {$result*$delta/6.0}]
 }
 
-# IntegralExpr --
+# integralExpr --
 #    Integrate an expression with "x" as the integrate according to the
 #    Simpson rule
 #
@@ -83,7 +88,7 @@
    return [expr {$result*$delta/6.0}]
 }
 
-# Integral2D --
+# integral2D --
 #    Integrate a given fucntion of two variables over a block,
 #    using bilinear interpolation (for this moment: block function)
 #
@@ -109,7 +114,7 @@
       set y [expr {$hydelta+double($j)*$ydelta}]
       for { set i 0 } { $i < $xnumber } { incr i } {
          set x           [expr {$hxdelta+double($i)*$xdelta}]
-         set func_value  [$func $x $y]
+         set func_value  [uplevel 1 [list $func $x $y]]
          set result      [expr {$result+$func_value}]
       }
    }
@@ -117,7 +122,7 @@
    return [expr {$result*$dxdy}]
 }
 
-# Integral3D --
+# integral3D --
 #    Integrate a given fucntion of two variables over a block,
 #    using trilinear interpolation (for this moment: block function)
 #
@@ -149,7 +154,7 @@
          set y [expr {$hydelta+double($j)*$ydelta}]
          for { set i 0 } { $i < $xnumber } { incr i } {
             set x           [expr {$hxdelta+double($i)*$xdelta}]
-            set func_value  [$func $x $y $z]
+            set func_value  [uplevel 1 [list $func $x $y $z]]
             set result      [expr {$result+$func_value}]
          }
       }
@@ -158,7 +163,7 @@
    return [expr {$result*$dxdydz}]
 }
 
-# EulerStep --
+# 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.
@@ -174,7 +179,7 @@
 #
 proc ::math::calculus::eulerStep { t tstep xvec func } {
 
-   set xderiv   [$func $t $xvec]
+   set xderiv   [uplevel 1 [list $func $t $xvec]]
    set result   {}
    foreach xv $xvec dx $xderiv {
       set xnew [expr {$xv+$tstep*$dx}]
@@ -184,7 +189,7 @@
    return $result
 }
 
-# HeunStep --
+# 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.
@@ -203,13 +208,14 @@
    #
    # Predictor step
    #
-   set xpred    [eulerStep $t $tstep $xvec $func]
+   set funcq    [uplevel 1 [list ::namespace which -command $func]]
+   set xpred    [eulerStep $t $tstep $xvec $funcq]
 
    #
    # Corrector step
    #
    set tcorr    [expr {$t+$tstep}]
-   set xcorr    [eulerStep $t     $tstep $xpred $func]
+   set xcorr    [eulerStep $t     $tstep $xpred $funcq]
 
    set result   {}
    foreach xv $xvec xc $xcorr {
@@ -220,7 +226,7 @@
    return $result
 }
 
-# RungeKuttaStep --
+# 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.
@@ -236,6 +242,8 @@
 #
 proc ::math::calculus::rungeKuttaStep { t tstep xvec func } {
 
+   set funcq    [uplevel 1 [list ::namespace which -command $func]]
+
    #
    # Four steps:
    # - k1 = tstep*func(t,x0)
@@ -247,25 +255,25 @@
    set tstep2   [expr {$tstep/2.0}]
    set tstep6   [expr {$tstep/6.0}]
 
-   set xk1      [$func $t $xvec]
+   set xk1      [$funcq $t $xvec]
    set xvec2    {}
    foreach x1 $xvec xv $xk1 {
       lappend xvec2 [expr {$x1+$tstep2*$xv}]
    }
 
-   set xk2      [$func [expr {$t+$tstep2}] $xvec2]
+   set xk2      [$funcq [expr {$t+$tstep2}] $xvec2]
    set xvec3    {}
    foreach x1 $xvec xv $xk2 {
       lappend xvec3 [expr {$x1+$tstep2*$xv}]
    }
 
-   set xk3      [$func [expr {$t+$tstep2}] $xvec3]
+   set xk3      [$funcq [expr {$t+$tstep2}] $xvec3]
    set xvec4    {}
    foreach x1 $xvec xv $xk3 {
       lappend xvec4 [expr {$x1+$tstep2*$xv}]
    }
 
-   set xk4      [$func [expr {$t+$tstep}] $xvec4]
+   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}]
@@ -275,7 +283,7 @@
    return $result
 }
 
-# BoundaryValueSecondOrder --
+# boundaryValueSecondOrder --
 #    Integrate a second-order differential equation and solve for
 #    given boundary values.
 #
@@ -301,6 +309,9 @@
 proc ::math::calculus::boundaryValueSecondOrder {
    coeff_func force_func leftbnd rightbnd nostep } {
 
+   set coeffq    [uplevel 1 [list ::namespace which -command $coeff_func]]
+   set forceq    [uplevel 1 [list ::namespace which -command $force_func]]
+
    if { [llength $leftbnd] != 2 || [llength $rightbnd] != 2 } {
       error "Boundary condition(s) incorrect"
    }
@@ -322,7 +333,7 @@
    set dvalue {}
 
    set x $xleft
-   foreach {A B C} [$coeff_func $x] { break }
+   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}]
@@ -333,7 +344,7 @@
       if { [expr {abs($x)-0.5*abs($xstep)}] < 0.0 } {
          set x 0.0
       }
-      foreach {A B C} [$coeff_func $x] { break }
+      foreach {A B C} [$coeffq $x] { break }
 
       set A2 0.0
       set B2 [expr {$A/$xstep-0.5*$B+0.5*$C*$xstep}]
@@ -352,7 +363,7 @@
          set x 0.0
       }
       lappend xvec   $x
-      lappend dvalue [expr {$xstep*[$force_func $x]}]
+      lappend dvalue [expr {$xstep*[$forceq $x]}]
    }
 
    #
@@ -374,7 +385,7 @@
    return $result
 }
 
-# SolveTriDiagonal --
+# solveTriDiagonal --
 #    Solve a system of equations Ax = b where A is a tridiagonal matrix
 #
 # Arguments:
@@ -429,5 +440,63 @@
    return $yvec
 }
 
-# Now we can announce our presence.
-package provide math::calculus 0.3
+# 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
+#
+proc ::math::calculus::newtonRaphson { func deriv initval } {
+   variable nr_maxiter
+   variable nr_tolerance
+
+   set funcq  [uplevel 1 [list ::namespace which -command $func]]
+   set derivq [uplevel 1 [list ::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
+#
+proc ::math::calculus::newtonRaphsonParameters { maxiter tolerance } {
+   variable nr_maxiter
+   variable nr_tolerance
+
+   if { $maxiter > 0 } {
+      set nr_maxiter $maxiter
+   }
+   if { $tolerance > 0 } {
+      set nr_tolerance $tolerance
+   }
+}
+
+# Now we can announce our presence
+package provide math::calculus 0.5
Index: modules/math/calculus.test
===================================================================
RCS file: /cvsroot/tcllib/tcllib/modules/math/calculus.test,v
retrieving revision 1.2
diff -u -r1.2 calculus.test
--- modules/math/calculus.test	8 May 2002 16:35:25 -0000	1.2
+++ modules/math/calculus.test	31 May 2002 20:59:52 -0000
@@ -13,16 +13,13 @@
 #
 # Simple test functions - exact result predictable!
 #
-proc ::const_func { x } {
+proc const_func { x } {
    return 1
 }
-proc ::linear_func { x } {
+proc linear_func { x } {
    return $x
 }
-proc ::downward_linear { x } {
-   return [expr {100.0-$x}]
-}
-proc ::downward_linear { x } {
+proc downward_linear { x } {
    return [expr {100.0-$x}]
 }
 
@@ -46,10 +43,10 @@
 } 5000.0
 
 
-proc ::const_func2d { x y } {
+proc const_func2d { x y } {
    return 1
 }
-proc ::linear_func2d { x y } {
+proc linear_func2d { x y } {
    return $x
 }
 
@@ -64,14 +61,19 @@
 } 250000.0
 
 #
+# Note:
+# Test cases for integral3D are missing!
+#
+
+#
 # Test cases: yet to be brought into the tcltest form!
 #
 
 # xvec should one long!
-proc ::const_func { t xvec } { return 1.0 }
+proc const_func { t xvec } { return 1.0 }
 
 # xvec should be two long!
-proc ::dampened_oscillator { t xvec } {
+proc dampened_oscillator { t xvec } {
    set x  [lindex $xvec 0]
    set x1 [lindex $xvec 1]
    return [list $x1 [expr {-$x1-$x}]]
@@ -103,13 +105,49 @@
 
 #
 # Boundary value problems:
-# use simple functions
 #
-proc ::coeffs { x } { return {1.0 0.0 0.0} }
-proc ::forces { x } { return 0.0 }
+proc coeffs { x } { return {1.0 0.0 0.0} }
+proc forces { x } { return 0.0 }
 
 puts [boundaryValueSecondOrder coeffs forces {0.0 1.0} {100.0 0.0} 10]
 puts [boundaryValueSecondOrder coeffs forces {0.0 0.0} {100.0 1.0} 10]
+
+#
+# Determining the root of an equation
+# use simple functions
+#
+proc func  { x } { expr {$x*$x-1.0} }
+proc deriv { x } { expr {2.0*$x} }
+
+test "NewtonRaphson-1.0" "Result should be 1" {
+   set result [newtonRaphson func deriv 2.0]
+   if { abs($result-1.0) < 0.0001 } {
+      set answer 1
+   }
+} 1
+test "NewtonRaphson-1.1" "Result should be -1" {
+   set result [newtonRaphson func deriv -0.5]
+   if { abs($result+1.0) < 0.0001 } {
+      set answer 1
+   }
+} 1
+
+proc func2  { x } { expr {$x*exp($x)-1.0} }
+proc deriv2 { x } { expr {exp($x)+$x*exp($x)} }
+
+test "NewtonRaphson-2.1" "Result should be nearly 0.56714" {
+   set result [newtonRaphson func2 deriv2 2.0]
+   if { abs($result-0.56714) < 0.0001 } {
+      set answer 1
+   }
+} 1
+
+test "NewtonRaphson-2.1" "Result should be nearly 0.56714" {
+   set result [newtonRaphson func2 deriv2 -0.5]
+   if { abs($result-0.56714) < 0.0001 } {
+      set answer 1
+   }
+} 1
 
 cleanupTests
 }
