[ Tcllib Table Of Contents  Tcllib Index ]
math::calculus(n) 0.8 "Tcl Math Library"
Name
math::calculus  Integration and ordinary differential equations
Synopsis
 package require Tcl 8.4
 package require math::calculus 0.8
 ::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 leftbnd 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:

The integration of a function over an interval

The numerical integration of a system of ordinary differential equations.

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, wellknown algorithms 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 remaining arguments are:
 nosteps

Number of steps in which the interval is divided.
 func

Function to be integrated. It should take one single argument.
 ::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 fourpoint quadrature formula. This results in an exact integration 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 threedimensional 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 GaussKronrod 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 argument.
 ?nosteps?

Number of steps in which the interval is divided. Defaults to 1.
 ::math::calculus::qk15_detailed xstart xend func nosteps

Determine the integral of the given function using the GaussKronrod 15 points quadrature rule. The interval for the integration is [xstart, xend]. The procedure returns a list of four values:

The estimate of the integral over the specified interval (I).

An estimate of the absolute error in I.

The estimate of the integral of the absolute value of the function over the interval.

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 argument.
 ?nosteps?

Number of steps in which the interval is divided. Defaults 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 RungeKutta 4th order.
 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 leftbnd 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 xcoordinate.
 force_func

Procedure returning the righthand side (D) as a function of the xcoordinate.
 leftbnd

A list of two values: the xcoordinate of the left boundary and the value at that boundary.
 rightbnd

A list of two values: the xcoordinate 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 xcoordinates and the approximated 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 righthandside
 ::math::calculus::newtonRaphson func deriv initval

Determine the root of an equation given by
func(x) = 0
using the method of NewtonRaphson. The procedure takes the following 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 NewtonRaphson 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 function contained in the interval [xb,xe]. The error in this estimate is of the order of eps*abs(xexb), the actual error may be slightly larger.
The method used is the socalled 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 function at the start has a different sign than the value at the end.
In contrast to NewtonRaphson there is no need for the computation 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.0e4)
Notes:
Several of the above procedures take the names of procedures as arguments. To avoid problems with the visibility of these procedures, the fullyqualified 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:
namespace eval ::mySpace { namespace export calcfunc proc calcfunc { x } { return $x } } # # Use a fullyqualified 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 secondorder boundary value problem:

Other types of boundary conditions (zero gradient, zero flux)

Other schematisation of the firstorder term (now central differences 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 firstorder 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 decaying 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.0e2 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. Please also report any ideas for enhancements you may have for either package and/or documentation.
See Also
romberg
Category
Mathematics
Copyright
Copyright © 2002,2003,2004 Arjen Markus