Provided by: tcllib_1.21+dfsg-1_all bug

NAME

       math::calculus - Integration and ordinary differential equations

SYNOPSIS

       package require Tcl  8.4

       package require math::calculus  0.8.2

       ::math::calculus::integral begin end nosteps func

       ::math::calculus::integralExpr begin end nosteps expression

       ::math::calculus::integral2D xinterval yinterval func

       ::math::calculus::integral2D_accurate xinterval yinterval func

       ::math::calculus::integral3D xinterval yinterval zinterval func

       ::math::calculus::integral3D_accurate xinterval yinterval zinterval func

       ::math::calculus::qk15 xstart xend func nosteps

       ::math::calculus::qk15_detailed xstart xend func nosteps

       ::math::calculus::eulerStep t tstep xvec func

       ::math::calculus::heunStep t tstep xvec func

       ::math::calculus::rungeKuttaStep t tstep xvec func

       ::math::calculus::boundaryValueSecondOrder coeff_func force_func leftbnd rightbnd nostep

       ::math::calculus::solveTriDiagonal acoeff bcoeff ccoeff dvalue

       ::math::calculus::newtonRaphson func deriv initval

       ::math::calculus::newtonRaphsonParameters maxiter tolerance

       ::math::calculus::regula_falsi f xb xe eps

_________________________________________________________________________________________________

DESCRIPTION

       This package implements several simple mathematical algorithms:

       •      The integration of a function over an interval

       •      The numerical integration of a system of ordinary differential equations.

       •      Estimating the root(s) of an equation of one variable.

       The  package  is  fully  implemented  in Tcl. No particular attention has been 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.

PROCEDURES

       This package defines the following public procedures:

       ::math::calculus::integral 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.

       ::math::calculus::integralExpr 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")

       ::math::calculus::integral2D xinterval yinterval func

       ::math::calculus::integral2D_accurate xinterval yinterval func
              The  commands  integral2D  and  integral2D_accurate  calculate  the  integral  of a
              function of two variables over the rectangle given by the first two arguments, each
              a  list of three items, the start and stop interval for the variable and the number
              of steps.

              The command integral2D evaluates the function at  the  centre  of  each  rectangle,
              whereas  the command integral2D_accurate uses a four-point quadrature formula. This
              results in an exact integration of polynomials of third degree or less.

              The function must take two arguments and return the function value.

       ::math::calculus::integral3D xinterval yinterval zinterval func

       ::math::calculus::integral3D_accurate xinterval yinterval zinterval func
              The  commands  integral3D  and  integral3D_accurate   are   the   three-dimensional
              equivalent  of  integral2D  and integral3D_accurate.  The function func takes three
              arguments and is integrated over the block in 3D space given by three intervals.

       ::math::calculus::qk15 xstart xend func nosteps
              Determine the integral of the given function  using  the  Gauss-Kronrod  15  points
              quadrature  rule.   The  returned  value  is  the estimate of the integral over the
              interval [xstart, xend].  The remaining arguments are:

              func   Function to be integrated. It should take one single argument.

              ?nosteps?
                     Number of steps in which the interval is divided. Defaults to 1.

       ::math::calculus::qk15_detailed xstart xend func nosteps
              Determine the integral of the given function  using  the  Gauss-Kronrod  15  points
              quadrature  rule.   The  interval  for  the  integration  is  [xstart,  xend].  The
              procedure returns a list of four values:

              •      The estimate of the integral over the specified interval (I).

              •      An estimate of the absolute error in I.

              •      The estimate of the integral of the absolute value of the function over  the
                     interval.

              •      The estimate of the integral of the absolute value of the function minus its
                     mean over the interval.

              The remaining arguments are:

              func   Function to be integrated. It should take one single argument.

              ?nosteps?
                     Number of steps in which the interval is divided. Defaults to 1.

       ::math::calculus::eulerStep 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).

       ::math::calculus::heunStep 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).

       ::math::calculus::rungeKuttaStep t tstep xvec func
              Set a single step  in  the  numerical  integration  of  a  system  of  differential
              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).

       ::math::calculus::boundaryValueSecondOrder coeff_func force_func leftbnd 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        dy
                       a(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.

              coeff_func
                     Procedure returning the three coefficients (A, B, C) of the equation, taking
                     as its one argument the x-coordinate.

              force_func
                     Procedure returning the  right-hand  side  (D)  as  a  function  of  the  x-
                     coordinate.

              leftbnd
                     A list of two values: the x-coordinate of the left boundary and the value at
                     that boundary.

              rightbnd
                     A 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.  The procedure returns
                     a list of x-coordinates and the approximated values of the solution.

       ::math::calculus::solveTriDiagonal 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 diagonal

              dvalue List of values on the righthand-side

       ::math::calculus::newtonRaphson func deriv initval
              Determine the root of an equation given by

                  func(x) = 0

              using the method of Newton-Raphson. The procedure takes the following 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

       ::math::calculus::newtonRaphsonParameters maxiter tolerance
              Set the numerical parameters for the Newton-Raphson method:

              maxiter
                     Maximum number of iteration steps (defaults to 20)

              tolerance
                     Relative precision (defaults to 0.001)

       ::math::calculus::regula_falsi f xb xe eps
              Return an estimate of the zero or one of the zeros of the function contained in the
              interval [xb,xe]. The error in this estimate is of the order of eps*abs(xe-xb), the
              actual error may be slightly larger.

              The  method  used  is  the so-called regula falsi or false position method. It is a
              straightforward implementation.  The  method  is  robust,  but  requires  that  the
              interval  brackets  a zero or at least an uneven number of zeros, so that the value
              of the function at the start has a different sign than the value at the end.

              In contrast to  Newton-Raphson  there  is  no  need  for  the  computation  of  the
              function's derivative.

              command f
                     Name  of the command that evaluates the function for which the zero is to be
                     returned

              float xb
                     Start of the interval in which the zero is supposed to lie

              float xe
                     End of the interval

              float eps
                     Relative allowed error (defaults to 1.0e-4)

       Notes:

       Several of the above procedures take the  names  of  procedures  as  arguments.  To  avoid
       problems  with  the  visibility  of  these  procedures,  the fully-qualified name of these
       procedures is determined inside the calculus 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 calcfunc
                     proc 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:

       •      Other types of boundary conditions (zero gradient, zero flux)

       •      Other schematisation of the first-order term (now central differences are used, but
              upstream differences might be useful too).

EXAMPLES

       Let us take a few simple examples:

       Integrate x over the interval [0,100] (20 steps):

              proc linear_func { x } { return $x }
              puts "Integral: [::math::calculus::integral 0 100 20 linear_func]"

       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 = 0

       can be split into a system of first-order equations:

              x' = y
              y' = -ry - wx

       Then this system can be solved with code like this:

              proc dampened_oscillator { 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.1
              for { set i 0 } { $i < 20 } { incr i } {
                 set result [::math::calculus::eulerStep $t $tstep $xvec dampened_oscillator]
                 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.0

                 set y [::math::calculus::boundaryValueSecondOrder \
                    coeffs force {0.0 1.0} [list $length 0.0] 100]

BUGS, IDEAS, FEEDBACK

       This document, and the package it describes,  will  undoubtedly  contain  bugs  and  other
       problems.   Please  report  such  in  the category math :: calculus of the Tcllib Trackers
       [http://core.tcl.tk/tcllib/reportlist].  Please also report any ideas for enhancements you
       may have for either package and/or documentation.

       When proposing code changes, please provide unified diffs, i.e the output of diff -u.

       Note further that attachments are strongly preferred over inlined patches. Attachments can
       be made by going to the Edit form of the ticket immediately after its creation,  and  then
       using the left-most button in the secondary navigation bar.

SEE ALSO

       romberg

KEYWORDS

       calculus, differential equations, integration, math, roots

CATEGORY

       Mathematics

COPYRIGHT

       Copyright (c) 2002,2003,2004 Arjen Markus