1math::calculus(n)              Tcl Math Library              math::calculus(n)
2
3
4
5______________________________________________________________________________
6

NAME

8       math::calculus - Integration and ordinary differential equations
9

SYNOPSIS

11       package require Tcl  8.4
12
13       package require math::calculus  0.8.2
14
15       ::math::calculus::integral begin end nosteps func
16
17       ::math::calculus::integralExpr begin end nosteps expression
18
19       ::math::calculus::integral2D xinterval yinterval func
20
21       ::math::calculus::integral2D_accurate xinterval yinterval func
22
23       ::math::calculus::integral3D xinterval yinterval zinterval func
24
25       ::math::calculus::integral3D_accurate   xinterval  yinterval  zinterval
26       func
27
28       ::math::calculus::qk15 xstart xend func nosteps
29
30       ::math::calculus::qk15_detailed xstart xend func nosteps
31
32       ::math::calculus::eulerStep t tstep xvec func
33
34       ::math::calculus::heunStep t tstep xvec func
35
36       ::math::calculus::rungeKuttaStep t tstep xvec func
37
38       ::math::calculus::boundaryValueSecondOrder coeff_func force_func  left‐
39       bnd rightbnd nostep
40
41       ::math::calculus::solveTriDiagonal acoeff bcoeff ccoeff dvalue
42
43       ::math::calculus::newtonRaphson func deriv initval
44
45       ::math::calculus::newtonRaphsonParameters maxiter tolerance
46
47       ::math::calculus::regula_falsi f xb xe eps
48
49______________________________________________________________________________
50

DESCRIPTION

52       This package implements several simple mathematical algorithms:
53
54       •      The integration of a function over an interval
55
56       •      The  numerical  integration of a system of ordinary differential
57              equations.
58
59       •      Estimating the root(s) of an equation of one variable.
60
61       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
65       This document describes the procedures and explains their usage.
66

PROCEDURES

68       This package defines the following public procedures:
69
70       ::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
75              nosteps
76                     Number of steps in which the interval is divided.
77
78              func   Function to be integrated. It should take one single  ar‐
79                     gument.
80
81
82       ::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
88              nosteps
89                     Number of steps in which the interval is divided.
90
91              expression
92                     Expression to be integrated. It should use  the  variable
93                     "x" as the only variable (the "integrate")
94
95
96       ::math::calculus::integral2D xinterval yinterval func
97
98       ::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
105              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
110              The function must take two arguments  and  return  the  function
111              value.
112
113       ::math::calculus::integral3D xinterval yinterval zinterval func
114
115       ::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
122       ::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
128              func   Function to be integrated. It should take one single  ar‐
129                     gument.
130
131              ?nosteps?
132                     Number  of  steps  in  which the interval is divided. De‐
133                     faults to 1.
134
135
136       ::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
142              •      The  estimate of the integral over the specified interval
143                     (I).
144
145              •      An estimate of the absolute error in I.
146
147              •      The estimate of the integral of the absolute value of the
148                     function over the interval.
149
150              •      The estimate of the integral of the absolute value of the
151                     function minus its mean over the interval.
152
153              The remaining arguments are:
154
155              func   Function to be integrated. It should take one single  ar‐
156                     gument.
157
158              ?nosteps?
159                     Number  of  steps  in  which the interval is divided. De‐
160                     faults to 1.
161
162
163       ::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
167              t      Value of the independent variable (typically time) at the
168                     beginning of the step.
169
170              tstep  Step size for the independent variable.
171
172              xvec   List (vector) of dependent values
173
174              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
178
179       ::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
183              t      Value of the independent variable (typically time) at the
184                     beginning of the step.
185
186              tstep  Step size for the independent variable.
187
188              xvec   List (vector) of dependent values
189
190              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
194
195       ::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
200              t      Value of the independent variable (typically time) at the
201                     beginning of the step.
202
203              tstep  Step size for the independent variable.
204
205              xvec   List (vector) of dependent values
206
207              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
211
212       ::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
218
219                       d      dy     d
220                       -- A(x)--  +  -- B(x)y + C(x)y  =  D(x)
221                       dx     dx     dx
222
223
224              Ordinarily, such an equation would be written as:
225
226
227                           d2y        dy
228                       a(x)---  + b(x)-- + c(x) y  =  D(x)
229                           dx2        dx
230
231
232              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
236
237                       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
241
242              Because of the differentiation, however, it is  much  easier  to
243              ask the user to provide the functions A, B and C directly.
244
245              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
250              force_func
251                     Procedure returning the right-hand side (D) as a function
252                     of the x-coordinate.
253
254              leftbnd
255                     A list of two values: the x-coordinate of the left bound‐
256                     ary and the value at that boundary.
257
258              rightbnd
259                     A  list  of  two  values:  the  x-coordinate of the right
260                     boundary and the value at that boundary.
261
262              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
266
267       ::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
271              acoeff List of values on the lower diagonal
272
273              bcoeff List of values on the main diagonal
274
275              ccoeff List of values on the upper diagonal
276
277              dvalue List of values on the righthand-side
278
279
280       ::math::calculus::newtonRaphson func deriv initval
281              Determine the root of an equation given by
282
283
284                  func(x) = 0
285
286
287              using the method of Newton-Raphson. The procedure takes the fol‐
288              lowing arguments:
289
290              func   Procedure that returns the value the function at x
291
292              deriv  Procedure that returns the derivative of the function  at
293                     x
294
295              initval
296                     Initial value for x
297
298
299       ::math::calculus::newtonRaphsonParameters maxiter tolerance
300              Set the numerical parameters for the Newton-Raphson method:
301
302              maxiter
303                     Maximum number of iteration steps (defaults to 20)
304
305              tolerance
306                     Relative precision (defaults to 0.001)
307
308       ::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
314              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
321              In  contrast to Newton-Raphson there is no need for the computa‐
322              tion of the function's derivative.
323
324              command f
325                     Name of the command that evaluates the function for which
326                     the zero is to be returned
327
328              float xb
329                     Start  of  the  interval in which the zero is supposed to
330                     lie
331
332              float xe
333                     End of the interval
334
335              float eps
336                     Relative allowed error (defaults to 1.0e-4)
337
338       Notes:
339
340       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
346
347                  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
369
370       Enhancements for the second-order boundary value problem:
371
372       •      Other types of boundary conditions (zero gradient, zero flux)
373
374       •      Other schematisation of the first-order term (now  central  dif‐
375              ferences  are  used,  but  upstream  differences might be useful
376              too).
377

EXAMPLES

379       Let us take a few simple examples:
380
381       Integrate x over the interval [0,100] (20 steps):
382
383
384              proc linear_func { x } { return $x }
385              puts "Integral: [::math::calculus::integral 0 100 20 linear_func]"
386
387       For simple functions, the alternative could be:
388
389
390              puts "Integral: [::math::calculus::integralExpr 0 100 20 {$x}]"
391
392       Do not forget the braces!
393
394       The differential equation for a dampened oscillator:
395
396              x'' + rx' + wx = 0
397
398
399       can be split into a system of first-order equations:
400
401              x' = y
402              y' = -ry - wx
403
404
405       Then this system can be solved with code like this:
406
407              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
413              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
423
424       Suppose we have the boundary value problem:
425
426                  Dy'' + ky = 0
427                  x = 0: y = 1
428                  x = L: y = 0
429
430
431       This boundary value problem could originate from the diffusion of a de‐
432       caying substance.
433
434       It can be solved with the following fragment:
435
436                 proc coeffs { x } { return [list $::Diff 0.0 $::decay] }
437                 proc force  { x } { return 0.0 }
438
439                 set Diff   1.0e-2
440                 set decay  0.0001
441                 set length 100.0
442
443                 set y [::math::calculus::boundaryValueSecondOrder \
444                    coeffs force {0.0 1.0} [list $length 0.0] 100]
445
446

BUGS, IDEAS, FEEDBACK

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
454       When proposing code changes, please provide unified diffs, i.e the out‐
455       put of diff -u.
456
457       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

SEE ALSO

463       romberg
464

KEYWORDS

466       calculus, differential equations, integration, math, roots
467

CATEGORY

469       Mathematics
470
472       Copyright (c) 2002,2003,2004 Arjen Markus
473
474
475
476
477tcllib                               0.8.2                   math::calculus(n)
Impressum