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.1
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 command f
272 Name of the command that evaluates the function for which
273 the zero is to be returned
274 float xb
276 Start of the interval in which the zero is supposed to
277 lie
278 float xe
280 End of the interval
281 float eps
283 Relative allowed error (defaults to 1.0e-4)
284 Notes:
286 Several of the above procedures take the names of procedures as argu‐
288 ments. To avoid problems with the visibility of these procedures, the
289 fully-qualified name of these procedures is determined inside the cal‐
290 culus routines. For the user this has only one consequence: the named
291 procedure must be visible in the calling procedure. For instance:
292 namespace eval ::mySpace {
294 namespace export calcfunc
295 proc calcfunc { x } { return $x }
296 }
297 #
298 # Use a fully-qualified name
299 #
300 namespace eval ::myCalc {
301 proc detIntegral { begin end } {
302 return [integral $begin $end 100 ::mySpace::calcfunc]
303 }
304 }
305 #
306 # Import the name
307 #
308 namespace eval ::myCalc {
309 namespace import ::mySpace::calcfunc
310 proc detIntegral { begin end } {
311 return [integral $begin $end 100 calcfunc]
312 }
313 }
314 Enhancements for the second-order boundary value problem:
317 · Other types of boundary conditions (zero gradient, zero flux)
319 · Other schematisation of the first-order term (now central dif‐
321 ferences are used, but upstream differences might be useful
322 too).
323
325 Let us take a few simple examples:
326 Integrate x over the interval [0,100] (20 steps):
328 proc linear_func { x } { return $x }
330 puts "Integral: [::math::calculus::integral 0 100 20 linear_func]"
331 For simple functions, the alternative could be:
333 puts "Integral: [::math::calculus::integralExpr 0 100 20 {$x}]"
335 Do not forget the braces!
337 The differential equation for a dampened oscillator:
339 x'' + rx' + wx = 0
341 can be split into a system of first-order equations:
344 x' = y
346 y' = -ry - wx
347 Then this system can be solved with code like this:
350 proc dampened_oscillator { t xvec } {
352 set x [lindex $xvec 0]
353 set x1 [lindex $xvec 1]
354 return [list $x1 [expr {-$x1-$x}]]
355 }
356 set xvec { 1.0 0.0 }
358 set t 0.0
359 set tstep 0.1
360 for { set i 0 } { $i < 20 } { incr i } {
361 set result [::math::calculus::eulerStep $t $tstep $xvec dampened_oscillator]
362 puts "Result ($t): $result"
363 set t [expr {$t+$tstep}]
364 set xvec $result
365 }
366 Suppose we have the boundary value problem:
369 Dy'' + ky = 0
371 x = 0: y = 1
372 x = L: y = 0
373 This boundary value problem could originate from the diffusion of a
376 decaying substance.
377 It can be solved with the following fragment:
379 proc coeffs { x } { return [list $::Diff 0.0 $::decay] }
381 proc force { x } { return 0.0 }
382 set Diff 1.0e-2
384 set decay 0.0001
385 set length 100.0
386 set y [::math::calculus::boundaryValueSecondOrder \
388 coeffs force {0.0 1.0} [list $length 0.0] 100]
389
392 This document, and the package it describes, will undoubtedly contain
393 bugs and other problems. Please report such in the category math ::
394 calculus of the Tcllib SF Trackers [http://source‐
395 forge.net/tracker/?group_id=12883]. Please also report any ideas for
396 enhancements you may have for either package and/or documentation.
397
399 romberg
400
402 calculus, differential equations, integration, math, roots
403
405 Copyright (c) 2002,2003,2004 Arjen Markus
406math 0.7.1 math::calculus(n)