NAME
math::calculus - Integration and ordinary differential equations
SYNOPSIS
package require TTccll 88 package require mmaatthh::::ccaallccuulluuss 00..55..11 ::::mmaatthh::::ccaallccuulluuss::::iinntteeggrraall begin end nosteps func ::::mmaatthh::::ccaallccuulluuss::::iinntteeggrraallEExxpprr begin end nosteps expression ::::mmaatthh::::ccaallccuulluuss::::iinntteeggrraall22DD xinterval yinterval func ::::mmaatthh::::ccaallccuulluuss::::iinntteeggrraall33DD xinterval yinterval zinterval func ::::mmaatthh::::ccaallccuulluuss::::eeuulleerrSStteepp t tstep xvec func ::::mmaatthh::::ccaallccuulluuss::::hheeuunnSStteepp t tstep xvec func ::::mmaatthh::::ccaallccuulluuss::::rruunnggeeKKuuttttaaSStteepp t tstep xvec func::::mmaatthh::::ccaallccuulluuss::::bboouunnddaarryyVVaalluueeSSeeccoonnddOOrrddeerr coefffunc forcefunc left-
bnd rightbnd nostep ::::mmaatthh::::ccaallccuulluuss::::ssoollvveeTTrriiDDiiaaggoonnaall acoeff bcoeff ccoeff dvalue ::::mmaatthh::::ccaallccuulluuss::::nneewwttoonnRRaapphhssoonn func deriv initval ::::mmaatthh::::ccaallccuulluuss::::nneewwttoonnRRaapphhssoonnPPaarraammeetteerrss maxiter toleranceDESCRIPTION
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 hasbeen 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. PPRROOCCEEDDUURREESS This package defines the following public procedures: ::::mmaatthh::::ccaallccuulluuss::::iinntteeggrraall 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. ::::mmaatthh::::ccaallccuulluuss::::iinntteeggrraallEExxpprr 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") ::::mmaatthh::::ccaallccuulluuss::::iinntteeggrraall22DD xinterval yinterval func The command iinntteeggrraall22DD calculates the integral of a function oftwo variables over the rectangle given by the first two argu-
ments, 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. ::::mmaatthh::::ccaallccuulluuss::::iinntteeggrraall33DD xinterval yinterval zinterval funcThe command IInntteeggrraall33DD is the three-dimensional equivalent of
iinntteeggrraall22DD. The function taking three arguments is integrated over the block in 3D space given by three intervals. ::::mmaatthh::::ccaallccuulluuss::::eeuulleerrSStteepp 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). ::::mmaatthh::::ccaallccuulluuss::::hheeuunnSStteepp 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). ::::mmaatthh::::ccaallccuulluuss::::rruunnggeeKKuuttttaaSStteepp t tstep xvec func Set a single step in the numerical integration of a system ofdifferential 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).::::mmaatthh::::ccaallccuulluuss::::bboouunnddaarryyVVaalluueeSSeeccoonnddOOrrddeerr coefffunc forcefunc 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 dya(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. coefffunc Procedure returning the three coefficients (A, B, C) ofthe equation, taking as its one argument the x-coordi-
nate. forcefuncProcedure returning the right-hand side (D) as a function
of the x-coordinate.
leftbndA list of two values: the x-coordinate of the left bound-
ary and the value at that boundary. rightbndA 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. Theprocedure returns a list of x-coordinates and the approx-
imated values of the solution. ::::mmaatthh::::ccaallccuulluuss::::ssoollvveeTTrriiDDiiaaggoonnaall 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 diagonaldvalue List of values on the righthand-side
::::mmaatthh::::ccaallccuulluuss::::nneewwttoonnRRaapphhssoonn func deriv initval Determine the root of an equation given by func(x) = 0using 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 ::::mmaatthh::::ccaallccuulluuss::::nneewwttoonnRRaapphhssoonnPPaarraammeetteerrss maxiter toleranceSet 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 argu-
ments. To avoid problems with the visibility of these procedures, thefully-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 calcfuncproc 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). EEXXAAMMPPLLEESS Let us take a few simple examples: Integrate x over the interval [0,100] (20 steps):proc linearfunc { x } { return $x }
puts "Integral: [::math::calculus::integral 0 100 20 linearfunc]"
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 = 0can be split into a system of first-order equations:
x' = yy' = -ry - wx
Then this system can be solved with code like this: proc dampenedoscillator { 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.1for { set i 0 } { $i < 20 } { incr i } {
set result [::math::calculus::eulerStep $t $tstep $xvec dampenedoscillator]
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.0set y [::math::calculus::boundaryValueSecondOrder \
coeffs force {0.0 1.0} [list $length 0.0] 100]
KKEEYYWWOORRDDSScalculus, differential equations, integration, math, roots
math 0.5.1 math::calculus(n)