1math::calculus(n) Tcl Math Library math::calculus(n)
2______________________________________________________________________________
6
8 math::calculus - Integration and ordinary differential equations
9
11 package require Tcl 8.4
12 package require math::calculus 0.7
14 ::math::calculus::integral begin end nosteps func
16 ::math::calculus::integralExpr begin end nosteps expression
18 ::math::calculus::integral2D xinterval yinterval func
20 ::math::calculus::integral2D_accurate xinterval yinterval func
22 ::math::calculus::integral3D xinterval yinterval zinterval func
24 ::math::calculus::integral3D_accurate xinterval yinterval zinterval
26 func
27 ::math::calculus::eulerStep t tstep xvec func
29 ::math::calculus::heunStep t tstep xvec func
31 ::math::calculus::rungeKuttaStep t tstep xvec func
33 ::math::calculus::boundaryValueSecondOrder coeff_func force_func left‐
35 bnd rightbnd nostep
36 ::math::calculus::solveTriDiagonal acoeff bcoeff ccoeff dvalue
38 ::math::calculus::newtonRaphson func deriv initval
40 ::math::calculus::newtonRaphsonParameters maxiter tolerance
42 ::math::calculus::regula_falsi f xb xe eps
44_________________________________________________________________
46
48 This package implements several simple mathematical algorithms:
49 · The integration of a function over an interval
51 · The numerical integration of a system of ordinary differential
53 equations.
54 · Estimating the root(s) of an equation of one variable.
56 The package is fully implemented in Tcl. No particular attention has
58 been paid to the accuracy of the calculations. Instead, well-known
59 algorithms have been used in a straightforward manner.
60 This document describes the procedures and explains their usage.
62
64 This package defines the following public procedures:
65 ::math::calculus::integral begin end nosteps func
67 Determine the integral of the given function using the Simpson
68 rule. The interval for the integration is [begin, end]. The
69 remaining arguments are:
70 nosteps
72 Number of steps in which the interval is divided.
73 func Function to be integrated. It should take one single
75 argument.
76 ::math::calculus::integralExpr begin end nosteps expression
79 Similar to the previous proc, this one determines the integral
80 of the given expression using the Simpson rule. The interval
81 for the integration is [begin, end]. The remaining arguments
82 are:
83 nosteps
85 Number of steps in which the interval is divided.
86 expression
88 Expression to be integrated. It should use the variable
89 "x" as the only variable (the "integrate")
90 ::math::calculus::integral2D xinterval yinterval func
93 ::math::calculus::integral2D_accurate xinterval yinterval func
95 The commands integral2D and integral2D_accurate calculate the
96 integral of a function of two variables over the rectangle given
97 by the first two arguments, each a list of three items, the
98 start and stop interval for the variable and the number of
99 steps.
100 The command integral2D evaluates the function at the centre of
102 each rectangle, whereas the command integral2D_accurate uses a
103 four-point quadrature formula. This results in an exact integra‐
104 tion of polynomials of third degree or less.
105 The function must take two arguments and return the function
107 value.
108 ::math::calculus::integral3D xinterval yinterval zinterval func
110 ::math::calculus::integral3D_accurate xinterval yinterval zinterval
112 func
113 The commands integral3D and integral3D_accurate are the three-
114 dimensional equivalent of integral2D and integral3D_accurate.
115 The function func takes three arguments and is integrated over
116 the block in 3D space given by three intervals.
117 ::math::calculus::eulerStep t tstep xvec func
119 Set a single step in the numerical integration of a system of
120 differential equations. The method used is Euler's.
121 t Value of the independent variable (typically time) at the
123 beginning of the step.
124 tstep Step size for the independent variable.
126 xvec List (vector) of dependent values
128 func Function of t and the dependent values, returning a list
130 of the derivatives of the dependent values. (The lengths
131 of xvec and the return value of "func" must match).
132 ::math::calculus::heunStep t tstep xvec func
135 Set a single step in the numerical integration of a system of
136 differential equations. The method used is Heun's.
137 t Value of the independent variable (typically time) at the
139 beginning of the step.
140 tstep Step size for the independent variable.
142 xvec List (vector) of dependent values
144 func Function of t and the dependent values, returning a list
146 of the derivatives of the dependent values. (The lengths
147 of xvec and the return value of "func" must match).
148 ::math::calculus::rungeKuttaStep t tstep xvec func
151 Set a single step in the numerical integration of a system of
152 differential equations. The method used is Runge-Kutta 4th
153 order.
154 t Value of the independent variable (typically time) at the
156 beginning of the step.
157 tstep Step size for the independent variable.
159 xvec List (vector) of dependent values
161 func Function of t and the dependent values, returning a list
163 of the derivatives of the dependent values. (The lengths
164 of xvec and the return value of "func" must match).
165 ::math::calculus::boundaryValueSecondOrder coeff_func force_func left‐
168 bnd rightbnd nostep
169 Solve a second order linear differential equation with boundary
170 values at two sides. The equation has to be of the form (the
171 "conservative" form):
172 d dy d
174 -- A(x)-- + -- B(x)y + C(x)y = D(x)
175 dx dx dx
176 Ordinarily, such an equation would be written as:
178 d2y dy
180 a(x)--- + b(x)-- + c(x) y = D(x)
181 dx2 dx
182 The first form is easier to discretise (by integrating over a
184 finite volume) than the second form. The relation between the
185 two forms is fairly straightforward:
186 A(x) = a(x)
188 B(x) = b(x) - a'(x)
189 C(x) = c(x) - B'(x) = c(x) - b'(x) + a''(x)
190 Because of the differentiation, however, it is much easier to
192 ask the user to provide the functions A, B and C directly.
193 coeff_func
195 Procedure returning the three coefficients (A, B, C) of
196 the equation, taking as its one argument the x-coordi‐
197 nate.
198 force_func
200 Procedure returning the right-hand side (D) as a function
201 of the x-coordinate.
202 leftbnd
204 A list of two values: the x-coordinate of the left bound‐
205 ary and the value at that boundary.
206 rightbnd
208 A list of two values: the x-coordinate of the right
209 boundary and the value at that boundary.
210 nostep Number of steps by which to discretise the interval. The
212 procedure returns a list of x-coordinates and the approx‐
213 imated values of the solution.
214 ::math::calculus::solveTriDiagonal acoeff bcoeff ccoeff dvalue
217 Solve a system of linear equations Ax = b with A a tridiagonal
218 matrix. Returns the solution as a list.
219 acoeff List of values on the lower diagonal
221 bcoeff List of values on the main diagonal
223 ccoeff List of values on the upper diagonal
225 dvalue List of values on the righthand-side
227 ::math::calculus::newtonRaphson func deriv initval
230 Determine the root of an equation given by
231 func(x) = 0
233 using the method of Newton-Raphson. The procedure takes the fol‐
235 lowing arguments:
236 func Procedure that returns the value the function at x
238 deriv Procedure that returns the derivative of the function at
240 x
241 initval
243 Initial value for x
244 ::math::calculus::newtonRaphsonParameters maxiter tolerance
247 Set the numerical parameters for the Newton-Raphson method:
248 maxiter
250 Maximum number of iteration steps (defaults to 20)
251 tolerance
253 Relative precision (defaults to 0.001)
254 ::math::calculus::regula_falsi f xb xe eps
256 Return an estimate of the zero or one of the zeros of the func‐
257 tion contained in the interval [xb,xe]. The error in this esti‐
258 mate is of the order of eps*abs(xe-xb), the actual error may be
259 slightly larger.
260 The method used is the so-called regula falsi or false position
262 method. It is a straightforward implementation. The method is
263 robust, but requires that the interval brackets a zero or at
264 least an uneven number of zeros, so that the value of the func‐
265 tion at the start has a different sign than the value at the
266 end.
267 In contrast to Newton-Raphson there is no need for the computa‐
269 tion of the function's derivative.
270 f command Name of the command that evaluates the function for
272 which the zero is to be returned
273 xb float Start of the interval in which the zero is supposed
275 to lie
276 xe float End of the interval
278 eps float Relative allowed error (defaults to 1.0e-4)
280 Notes:
282 Several of the above procedures take the names of procedures as argu‐
284 ments. To avoid problems with the visibility of these procedures, the
285 fully-qualified name of these procedures is determined inside the cal‐
286 culus routines. For the user this has only one consequence: the named
287 procedure must be visible in the calling procedure. For instance:
288 namespace eval ::mySpace {
290 namespace export calcfunc
291 proc calcfunc { x } { return $x }
292 }
293 #
294 # Use a fully-qualified name
295 #
296 namespace eval ::myCalc {
297 proc detIntegral { begin end } {
298 return [integral $begin $end 100 ::mySpace::calcfunc]
299 }
300 }
301 #
302 # Import the name
303 #
304 namespace eval ::myCalc {
305 namespace import ::mySpace::calcfunc
306 proc detIntegral { begin end } {
307 return [integral $begin $end 100 calcfunc]
308 }
309 }
310 Enhancements for the second-order boundary value problem:
313 · Other types of boundary conditions (zero gradient, zero flux)
315 · Other schematisation of the first-order term (now central dif‐
317 ferences are used, but upstream differences might be useful
318 too).
319
321 Let us take a few simple examples:
322 Integrate x over the interval [0,100] (20 steps):
324 proc linear_func { x } { return $x }
326 puts "Integral: [::math::calculus::integral 0 100 20 linear_func]"
327 For simple functions, the alternative could be:
329 puts "Integral: [::math::calculus::integralExpr 0 100 20 {$x}]"
331 Do not forget the braces!
333 The differential equation for a dampened oscillator:
335 x'' + rx' + wx = 0
337 can be split into a system of first-order equations:
340 x' = y
342 y' = -ry - wx
343 Then this system can be solved with code like this:
346 proc dampened_oscillator { t xvec } {
348 set x [lindex $xvec 0]
349 set x1 [lindex $xvec 1]
350 return [list $x1 [expr {-$x1-$x}]]
351 }
352 set xvec { 1.0 0.0 }
354 set t 0.0
355 set tstep 0.1
356 for { set i 0 } { $i < 20 } { incr i } {
357 set result [::math::calculus::eulerStep $t $tstep $xvec dampened_oscillator]
358 puts "Result ($t): $result"
359 set t [expr {$t+$tstep}]
360 set xvec $result
361 }
362 Suppose we have the boundary value problem:
365 Dy'' + ky = 0
367 x = 0: y = 1
368 x = L: y = 0
369 This boundary value problem could originate from the diffusion of a
372 decaying substance.
373 It can be solved with the following fragment:
375 proc coeffs { x } { return [list $::Diff 0.0 $::decay] }
377 proc force { x } { return 0.0 }
378 set Diff 1.0e-2
380 set decay 0.0001
381 set length 100.0
382 set y [::math::calculus::boundaryValueSecondOrder \
384 coeffs force {0.0 1.0} [list $length 0.0] 100]
385
388 romberg
389
391 calculus, differential equations, integration, math, roots
392
394 Copyright (c) 2002,2003,2004 Arjen Markus
395math 0.7 math::calculus(n)