Attachment "calc.patch" to
ticket [553773ffff]
added by
dgp
2002-06-01 04:01:51.
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 ([email protected]) 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: [email protected]).
@@ -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
}
