math::calculus(n) 1.0 math "Math"
math::calculus - Integration and ordinary differential equations
package require Tcl 8
package require math::calculus 0.5
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, well-known
algorithms have been used in a straightforward manner.
This document describes the procedures and explains their usage.
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
-
The command integral2D 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.
The currently implemented integration is simple: the function is
evaluated at the centre of each rectangle and the content of
this block is added to the integral. In future this will be
replaced by a bilinear interpolation.
The function must take two arguments and return the function
value.
- ::math::calculus::integral3D xinterval yinterval zinterval func
-
The command Integral3D is the three-dimensional
equivalent of integral2D. The function taking three
arguments is integrated over the block in 3D space given by three
intervals.
- ::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 tstep xvec func
-
Set a single step in the numerical integration of a system of
differential equations. The method used is Runge-Kutta 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 x-coordinate.
- 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 boundary 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 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 righthand-side
- ::math::calculus::newtonRaphson func deriv initval
-
Determine the root of an equation given by
using the method of Newton-Raphson. 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 Newton-Raphson method:
- maxiter
-
Maximum number of iteration steps (defaults to 20)
- tolerance
-
Relative precision (defaults to 0.001)
Notes:
Several of the above procedures take the names of procedures as
arguments. To avoid problems with the 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:
| |
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:
-
Other types of boundary conditions (zero gradient, zero flux)
-
Other schematisation of the first-order term (now central
differences are used, but upstream differences might be useful too).
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:
can be split into a system of first-order equations:
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.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]
|
calculus, differential equations, integration, math, roots