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 (c) 2002,2003,2004 Arjen Markus