math::calculus(3tcl) Tcl Math Library math::calculus(3tcl)
______________________________________________________________________________
NAME
math::calculus - Integration and ordinary differential equations
SYNOPSIS
package require Tcl 8.4
package require math::calculus 0.8.2
::math::calculus::integral begin end nosteps func
::math::calculus::integralExpr begin end nosteps expression
::math::calculus::integral2D xinterval yinterval func
::math::calculus::integral2D_accurate xinterval yinterval func
::math::calculus::integral3D xinterval yinterval zinterval func
::math::calculus::integral3D_accurate xinterval yinterval zinterval
func
::math::calculus::qk15 xstart xend func nosteps
::math::calculus::qk15_detailed xstart xend func nosteps
::math::calculus::eulerStep t tstep xvec func
::math::calculus::heunStep t tstep xvec func
::math::calculus::rungeKuttaStep t tstep xvec func
::math::calculus::boundaryValueSecondOrder coeff_func force_func left-
bnd rightbnd nostep
::math::calculus::solveTriDiagonal acoeff bcoeff ccoeff dvalue
::math::calculus::newtonRaphson func deriv initval
::math::calculus::newtonRaphsonParameters maxiter tolerance
::math::calculus::regula_falsi f xb xe eps
______________________________________________________________________________
DESCRIPTION
This package implements several simple mathematical algorithms:
o The integration of a function over an interval
o The numerical integration of a system of ordinary differential
equations.
o Estimating the root(s) of an equation of one variable.
The package is fully implemented in Tcl. No particular attention has
been paid to the accuracy of the calculations. Instead, well-known al-
gorithms have been used in a straightforward manner.
This document describes the procedures and explains their usage.
PROCEDURES
This package defines the following public procedures:
::math::calculus::integral begin end nosteps func
Determine the integral of the given function using the Simpson
rule. The interval for the integration is [begin, end]. The re-
maining arguments are:
nosteps
Number of steps in which the interval is divided.
func Function to be integrated. It should take one single ar-
gument.
::math::calculus::integralExpr begin end nosteps expression
Similar to the previous proc, this one determines the integral
of the given expression using the Simpson rule. The interval
for the integration is [begin, end]. The remaining arguments
are:
nosteps
Number of steps in which the interval is divided.
expression
Expression to be integrated. It should use the variable
"x" as the only variable (the "integrate")
::math::calculus::integral2D xinterval yinterval func
::math::calculus::integral2D_accurate xinterval yinterval func
The commands integral2D and integral2D_accurate calculate 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.
The command integral2D evaluates the function at the centre of
each rectangle, whereas the command integral2D_accurate uses a
four-point quadrature formula. This results in an exact integra-
tion of polynomials of third degree or less.
The function must take two arguments and return the function
value.
::math::calculus::integral3D xinterval yinterval zinterval func
::math::calculus::integral3D_accurate xinterval yinterval zinterval
func
The commands integral3D and integral3D_accurate are the three-
dimensional equivalent of integral2D and integral3D_accurate.
The function func takes three arguments and is integrated over
the block in 3D space given by three intervals.
::math::calculus::qk15 xstart xend func nosteps
Determine the integral of the given function using the Gauss-
Kronrod 15 points quadrature rule. The returned value is the
estimate of the integral over the interval [xstart, xend]. The
remaining arguments are:
func Function to be integrated. It should take one single ar-
gument.
?nosteps?
Number of steps in which the interval is divided. De-
faults to 1.
::math::calculus::qk15_detailed xstart xend func nosteps
Determine the integral of the given function using the Gauss-
Kronrod 15 points quadrature rule. The interval for the inte-
gration is [xstart, xend]. The procedure returns a list of four
values:
o The estimate of the integral over the specified interval
(I).
o An estimate of the absolute error in I.
o The estimate of the integral of the absolute value of the
function over the interval.
o The estimate of the integral of the absolute value of the
function minus its mean over the interval.
The remaining arguments are:
func Function to be integrated. It should take one single ar-
gument.
?nosteps?
Number of steps in which the interval is divided. De-
faults to 1.
::math::calculus::eulerStep t tstep xvec func
Set a single step in the numerical integration of a system of
differential equations. The method used is Euler's.
t Value of the independent variable (typically time) at the
beginning of the step.
tstep Step size for the independent variable.
xvec List (vector) of dependent values
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).
::math::calculus::heunStep t tstep xvec func
Set a single step in the numerical integration of a system of
differential equations. The method used is Heun's.
t Value of the independent variable (typically time) at the
beginning of the step.
tstep Step size for the independent variable.
xvec List (vector) of dependent values
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).
::math::calculus::rungeKuttaStep t tstep xvec func
Set a single step in the numerical integration of a system of
differential equations. The method used is Runge-Kutta 4th or-
der.
t Value of the independent variable (typically time) at the
beginning of the step.
tstep Step size for the independent variable.
xvec List (vector) of dependent values
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).
::math::calculus::boundaryValueSecondOrder coeff_func force_func left-
bnd rightbnd nostep
Solve a second order linear differential equation with boundary
values at two sides. The equation has to be of the form (the
"conservative" form):
d dy d
-- A(x)-- + -- B(x)y + C(x)y = D(x)
dx dx dx
Ordinarily, such an equation would be written as:
d2y dy
a(x)--- + b(x)-- + c(x) y = D(x)
dx2 dx
The first form is easier to discretise (by integrating over a
finite volume) than the second form. The relation between the
two forms is fairly straightforward:
A(x) = a(x)
B(x) = b(x) - a'(x)
C(x) = c(x) - B'(x) = c(x) - b'(x) + a''(x)
Because of the differentiation, however, it is much easier to
ask the user to provide the functions A, B and C directly.
coeff_func
Procedure returning the three coefficients (A, B, C) of
the equation, taking as its one argument the x-coordi-
nate.
force_func
Procedure returning the right-hand side (D) as a function
of the x-coordinate.
leftbnd
A list of two values: the x-coordinate of the left bound-
ary and the value at that boundary.
rightbnd
A list of two values: the x-coordinate of the right
boundary and the value at that boundary.
nostep Number of steps by which to discretise the interval. The
procedure returns a list of x-coordinates and the approx-
imated values of the solution.
::math::calculus::solveTriDiagonal acoeff bcoeff ccoeff dvalue
Solve a system of linear equations Ax = b with A a tridiagonal
matrix. Returns the solution as a list.
acoeff List of values on the lower diagonal
bcoeff List of values on the main diagonal
ccoeff List of values on the upper diagonal
dvalue List of values on the righthand-side
::math::calculus::newtonRaphson func deriv initval
Determine the root of an equation given by
func(x) = 0
using the method of Newton-Raphson. The procedure takes the fol-
lowing arguments:
func Procedure that returns the value the function at x
deriv Procedure that returns the derivative of the function at
x
initval
Initial value for x
::math::calculus::newtonRaphsonParameters maxiter tolerance
Set the numerical parameters for the Newton-Raphson method:
maxiter
Maximum number of iteration steps (defaults to 20)
tolerance
Relative precision (defaults to 0.001)
::math::calculus::regula_falsi f xb xe eps
Return an estimate of the zero or one of the zeros of the func-
tion contained in the interval [xb,xe]. The error in this esti-
mate is of the order of eps*abs(xe-xb), the actual error may be
slightly larger.
The method used is the so-called regula falsi or false position
method. It is a straightforward implementation. The method is
robust, but requires that the interval brackets a zero or at
least an uneven number of zeros, so that the value of the func-
tion at the start has a different sign than the value at the
end.
In contrast to Newton-Raphson there is no need for the computa-
tion of the function's derivative.
command f
Name of the command that evaluates the function for which
the zero is to be returned
float xb
Start of the interval in which the zero is supposed to
lie
float xe
End of the interval
float eps
Relative allowed error (defaults to 1.0e-4)
Notes:
Several of the above procedures take the names of procedures as argu-
ments. To avoid problems with the visibility of these procedures, the
fully-qualified name of these procedures is determined inside the cal-
culus routines. For the user this has only one consequence: the named
procedure must be visible in the calling procedure. For instance:
namespace eval ::mySpace {
namespace export calcfunc
proc calcfunc { x } { return $x }
}
#
# Use a fully-qualified name
#
namespace eval ::myCalc {
proc detIntegral { begin end } {
return [integral $begin $end 100 ::mySpace::calcfunc]
}
}
#
# Import the name
#
namespace eval ::myCalc {
namespace import ::mySpace::calcfunc
proc detIntegral { begin end } {
return [integral $begin $end 100 calcfunc]
}
}
Enhancements for the second-order boundary value problem:
o Other types of boundary conditions (zero gradient, zero flux)
o Other schematisation of the first-order term (now central dif-
ferences are used, but upstream differences might be useful
too).
EXAMPLES
Let us take a few simple examples:
Integrate x over the interval [0,100] (20 steps):
proc linear_func { x } { return $x }
puts "Integral: [::math::calculus::integral 0 100 20 linear_func]"
For simple functions, the alternative could be:
puts "Integral: [::math::calculus::integralExpr 0 100 20 {$x}]"
Do not forget the braces!
The differential equation for a dampened oscillator:
x'' + rx' + wx = 0
can be split into a system of first-order equations:
x' = y
y' = -ry - wx
Then this system can be solved with code like this:
proc dampened_oscillator { t xvec } {
set x [lindex $xvec 0]
set x1 [lindex $xvec 1]
return [list $x1 [expr {-$x1-$x}]]
}
set xvec { 1.0 0.0 }
set t 0.0
set tstep 0.1
for { set i 0 } { $i < 20 } { incr i } {
set result [::math::calculus::eulerStep $t $tstep $xvec dampened_oscillator]
puts "Result ($t): $result"
set t [expr {$t+$tstep}]
set xvec $result
}
Suppose we have the boundary value problem:
Dy'' + ky = 0
x = 0: y = 1
x = L: y = 0
This boundary value problem could originate from the diffusion of a de-
caying substance.
It can be solved with the following fragment:
proc coeffs { x } { return [list $::Diff 0.0 $::decay] }
proc force { x } { return 0.0 }
set Diff 1.0e-2
set decay 0.0001
set length 100.0
set y [::math::calculus::boundaryValueSecondOrder \
coeffs force {0.0 1.0} [list $length 0.0] 100]
BUGS, IDEAS, FEEDBACK
This document, and the package it describes, will undoubtedly contain
bugs and other problems. Please report such in the category math ::
calculus of the Tcllib Trackers [http://core.tcl.tk/tcllib/reportlist].
Please also report any ideas for enhancements you may have for either
package and/or documentation.
When proposing code changes, please provide unified diffs, i.e the out-
put of diff -u.
Note further that attachments are strongly preferred over inlined
patches. Attachments can be made by going to the Edit form of the
ticket immediately after its creation, and then using the left-most
button in the secondary navigation bar.
SEE ALSO
romberg
KEYWORDS
calculus, differential equations, integration, math, roots
CATEGORY
Mathematics
COPYRIGHT
Copyright (c) 2002,2003,2004 Arjen Markus
tcllib 0.8.2 math::calculus(3tcl)