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.8.2
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::qk15 xstart xend func nosteps
29 ::math::calculus::qk15_detailed xstart xend func nosteps
31 ::math::calculus::eulerStep t tstep xvec func
33 ::math::calculus::heunStep t tstep xvec func
35 ::math::calculus::rungeKuttaStep t tstep xvec func
37 ::math::calculus::boundaryValueSecondOrder coeff_func force_func left‐
39 bnd rightbnd nostep
40 ::math::calculus::solveTriDiagonal acoeff bcoeff ccoeff dvalue
42 ::math::calculus::newtonRaphson func deriv initval
44 ::math::calculus::newtonRaphsonParameters maxiter tolerance
46 ::math::calculus::regula_falsi f xb xe eps
48______________________________________________________________________________
50
52 This package implements several simple mathematical algorithms:
53 • The integration of a function over an interval
55 • The numerical integration of a system of ordinary differential
57 equations.
58 • Estimating the root(s) of an equation of one variable.
60 The package is fully implemented in Tcl. No particular attention has
62 been paid to the accuracy of the calculations. Instead, well-known al‐
63 gorithms have been used in a straightforward manner.
64 This document describes the procedures and explains their usage.
66
68 This package defines the following public procedures:
69 ::math::calculus::integral begin end nosteps func
71 Determine the integral of the given function using the Simpson
72 rule. The interval for the integration is [begin, end]. The re‐
73 maining arguments are:
74 nosteps
76 Number of steps in which the interval is divided.
77 func Function to be integrated. It should take one single ar‐
79 gument.
80 ::math::calculus::integralExpr begin end nosteps expression
83 Similar to the previous proc, this one determines the integral
84 of the given expression using the Simpson rule. The interval
85 for the integration is [begin, end]. The remaining arguments
86 are:
87 nosteps
89 Number of steps in which the interval is divided.
90 expression
92 Expression to be integrated. It should use the variable
93 "x" as the only variable (the "integrate")
94 ::math::calculus::integral2D xinterval yinterval func
97 ::math::calculus::integral2D_accurate xinterval yinterval func
99 The commands integral2D and integral2D_accurate calculate the
100 integral of a function of two variables over the rectangle given
101 by the first two arguments, each a list of three items, the
102 start and stop interval for the variable and the number of
103 steps.
104 The command integral2D evaluates the function at the centre of
106 each rectangle, whereas the command integral2D_accurate uses a
107 four-point quadrature formula. This results in an exact integra‐
108 tion of polynomials of third degree or less.
109 The function must take two arguments and return the function
111 value.
112 ::math::calculus::integral3D xinterval yinterval zinterval func
114 ::math::calculus::integral3D_accurate xinterval yinterval zinterval
116 func
117 The commands integral3D and integral3D_accurate are the three-
118 dimensional equivalent of integral2D and integral3D_accurate.
119 The function func takes three arguments and is integrated over
120 the block in 3D space given by three intervals.
121 ::math::calculus::qk15 xstart xend func nosteps
123 Determine the integral of the given function using the Gauss-
124 Kronrod 15 points quadrature rule. The returned value is the
125 estimate of the integral over the interval [xstart, xend]. The
126 remaining arguments are:
127 func Function to be integrated. It should take one single ar‐
129 gument.
130 ?nosteps?
132 Number of steps in which the interval is divided. De‐
133 faults to 1.
134 ::math::calculus::qk15_detailed xstart xend func nosteps
137 Determine the integral of the given function using the Gauss-
138 Kronrod 15 points quadrature rule. The interval for the inte‐
139 gration is [xstart, xend]. The procedure returns a list of four
140 values:
141 • The estimate of the integral over the specified interval
143 (I).
144 • An estimate of the absolute error in I.
146 • The estimate of the integral of the absolute value of the
148 function over the interval.
149 • The estimate of the integral of the absolute value of the
151 function minus its mean over the interval.
152 The remaining arguments are:
154 func Function to be integrated. It should take one single ar‐
156 gument.
157 ?nosteps?
159 Number of steps in which the interval is divided. De‐
160 faults to 1.
161 ::math::calculus::eulerStep t tstep xvec func
164 Set a single step in the numerical integration of a system of
165 differential equations. The method used is Euler's.
166 t Value of the independent variable (typically time) at the
168 beginning of the step.
169 tstep Step size for the independent variable.
171 xvec List (vector) of dependent values
173 func Function of t and the dependent values, returning a list
175 of the derivatives of the dependent values. (The lengths
176 of xvec and the return value of "func" must match).
177 ::math::calculus::heunStep t tstep xvec func
180 Set a single step in the numerical integration of a system of
181 differential equations. The method used is Heun's.
182 t Value of the independent variable (typically time) at the
184 beginning of the step.
185 tstep Step size for the independent variable.
187 xvec List (vector) of dependent values
189 func Function of t and the dependent values, returning a list
191 of the derivatives of the dependent values. (The lengths
192 of xvec and the return value of "func" must match).
193 ::math::calculus::rungeKuttaStep t tstep xvec func
196 Set a single step in the numerical integration of a system of
197 differential equations. The method used is Runge-Kutta 4th or‐
198 der.
199 t Value of the independent variable (typically time) at the
201 beginning of the step.
202 tstep Step size for the independent variable.
204 xvec List (vector) of dependent values
206 func Function of t and the dependent values, returning a list
208 of the derivatives of the dependent values. (The lengths
209 of xvec and the return value of "func" must match).
210 ::math::calculus::boundaryValueSecondOrder coeff_func force_func left‐
213 bnd rightbnd nostep
214 Solve a second order linear differential equation with boundary
215 values at two sides. The equation has to be of the form (the
216 "conservative" form):
217 d dy d
220 -- A(x)-- + -- B(x)y + C(x)y = D(x)
221 dx dx dx
222 Ordinarily, such an equation would be written as:
225 d2y dy
228 a(x)--- + b(x)-- + c(x) y = D(x)
229 dx2 dx
230 The first form is easier to discretise (by integrating over a
233 finite volume) than the second form. The relation between the
234 two forms is fairly straightforward:
235 A(x) = a(x)
238 B(x) = b(x) - a'(x)
239 C(x) = c(x) - B'(x) = c(x) - b'(x) + a''(x)
240 Because of the differentiation, however, it is much easier to
243 ask the user to provide the functions A, B and C directly.
244 coeff_func
246 Procedure returning the three coefficients (A, B, C) of
247 the equation, taking as its one argument the x-coordi‐
248 nate.
249 force_func
251 Procedure returning the right-hand side (D) as a function
252 of the x-coordinate.
253 leftbnd
255 A list of two values: the x-coordinate of the left bound‐
256 ary and the value at that boundary.
257 rightbnd
259 A list of two values: the x-coordinate of the right
260 boundary and the value at that boundary.
261 nostep Number of steps by which to discretise the interval. The
263 procedure returns a list of x-coordinates and the approx‐
264 imated values of the solution.
265 ::math::calculus::solveTriDiagonal acoeff bcoeff ccoeff dvalue
268 Solve a system of linear equations Ax = b with A a tridiagonal
269 matrix. Returns the solution as a list.
270 acoeff List of values on the lower diagonal
272 bcoeff List of values on the main diagonal
274 ccoeff List of values on the upper diagonal
276 dvalue List of values on the righthand-side
278 ::math::calculus::newtonRaphson func deriv initval
281 Determine the root of an equation given by
282 func(x) = 0
285 using the method of Newton-Raphson. The procedure takes the fol‐
288 lowing arguments:
289 func Procedure that returns the value the function at x
291 deriv Procedure that returns the derivative of the function at
293 x
294 initval
296 Initial value for x
297 ::math::calculus::newtonRaphsonParameters maxiter tolerance
300 Set the numerical parameters for the Newton-Raphson method:
301 maxiter
303 Maximum number of iteration steps (defaults to 20)
304 tolerance
306 Relative precision (defaults to 0.001)
307 ::math::calculus::regula_falsi f xb xe eps
309 Return an estimate of the zero or one of the zeros of the func‐
310 tion contained in the interval [xb,xe]. The error in this esti‐
311 mate is of the order of eps*abs(xe-xb), the actual error may be
312 slightly larger.
313 The method used is the so-called regula falsi or false position
315 method. It is a straightforward implementation. The method is
316 robust, but requires that the interval brackets a zero or at
317 least an uneven number of zeros, so that the value of the func‐
318 tion at the start has a different sign than the value at the
319 end.
320 In contrast to Newton-Raphson there is no need for the computa‐
322 tion of the function's derivative.
323 command f
325 Name of the command that evaluates the function for which
326 the zero is to be returned
327 float xb
329 Start of the interval in which the zero is supposed to
330 lie
331 float xe
333 End of the interval
334 float eps
336 Relative allowed error (defaults to 1.0e-4)
337 Notes:
339 Several of the above procedures take the names of procedures as argu‐
341 ments. To avoid problems with the visibility of these procedures, the
342 fully-qualified name of these procedures is determined inside the cal‐
343 culus routines. For the user this has only one consequence: the named
344 procedure must be visible in the calling procedure. For instance:
345 namespace eval ::mySpace {
348 namespace export calcfunc
349 proc calcfunc { x } { return $x }
350 }
351 #
352 # Use a fully-qualified name
353 #
354 namespace eval ::myCalc {
355 proc detIntegral { begin end } {
356 return [integral $begin $end 100 ::mySpace::calcfunc]
357 }
358 }
359 #
360 # Import the name
361 #
362 namespace eval ::myCalc {
363 namespace import ::mySpace::calcfunc
364 proc detIntegral { begin end } {
365 return [integral $begin $end 100 calcfunc]
366 }
367 }
368 Enhancements for the second-order boundary value problem:
371 • Other types of boundary conditions (zero gradient, zero flux)
373 • Other schematisation of the first-order term (now central dif‐
375 ferences are used, but upstream differences might be useful
376 too).
377
379 Let us take a few simple examples:
380 Integrate x over the interval [0,100] (20 steps):
382 proc linear_func { x } { return $x }
385 puts "Integral: [::math::calculus::integral 0 100 20 linear_func]"
386 For simple functions, the alternative could be:
388 puts "Integral: [::math::calculus::integralExpr 0 100 20 {$x}]"
391 Do not forget the braces!
393 The differential equation for a dampened oscillator:
395 x'' + rx' + wx = 0
397 can be split into a system of first-order equations:
400 x' = y
402 y' = -ry - wx
403 Then this system can be solved with code like this:
406 proc dampened_oscillator { t xvec } {
408 set x [lindex $xvec 0]
409 set x1 [lindex $xvec 1]
410 return [list $x1 [expr {-$x1-$x}]]
411 }
412 set xvec { 1.0 0.0 }
414 set t 0.0
415 set tstep 0.1
416 for { set i 0 } { $i < 20 } { incr i } {
417 set result [::math::calculus::eulerStep $t $tstep $xvec dampened_oscillator]
418 puts "Result ($t): $result"
419 set t [expr {$t+$tstep}]
420 set xvec $result
421 }
422 Suppose we have the boundary value problem:
425 Dy'' + ky = 0
427 x = 0: y = 1
428 x = L: y = 0
429 This boundary value problem could originate from the diffusion of a de‐
432 caying substance.
433 It can be solved with the following fragment:
435 proc coeffs { x } { return [list $::Diff 0.0 $::decay] }
437 proc force { x } { return 0.0 }
438 set Diff 1.0e-2
440 set decay 0.0001
441 set length 100.0
442 set y [::math::calculus::boundaryValueSecondOrder \
444 coeffs force {0.0 1.0} [list $length 0.0] 100]
445
448 This document, and the package it describes, will undoubtedly contain
449 bugs and other problems. Please report such in the category math ::
450 calculus of the Tcllib Trackers [http://core.tcl.tk/tcllib/reportlist].
451 Please also report any ideas for enhancements you may have for either
452 package and/or documentation.
453 When proposing code changes, please provide unified diffs, i.e the out‐
455 put of diff -u.
456 Note further that attachments are strongly preferred over inlined
458 patches. Attachments can be made by going to the Edit form of the
459 ticket immediately after its creation, and then using the left-most
460 button in the secondary navigation bar.
461
463 romberg
464
466 calculus, differential equations, integration, math, roots
467
469 Mathematics
470
472 Copyright (c) 2002,2003,2004 Arjen Markus
473tcllib 0.8.2 math::calculus(n)