<html><head>
<title>math::calculus - Tcl Math Library</title>
<style type="text/css"><!--
    HTML {
	background: 	#FFFFFF;
	color: 		black;
    }
    BODY {
	background: 	#FFFFFF;
	color:	 	black;
    }
    DIV.doctools {
	margin-left:	10%;
	margin-right:	10%;
    }
    DIV.doctools H1,DIV.doctools H2 {
	margin-left:	-5%;
    }
    H1, H2, H3, H4 {
	margin-top: 	1em;
	font-family:	sans-serif;
	font-size:	large;
	color:		#005A9C;
	background: 	transparent;
	text-align:		left;
    }
    H1.title {
	text-align: center;
    }
    UL,OL {
	margin-right: 0em;
	margin-top: 3pt;
	margin-bottom: 3pt;
    }
    UL LI {
	list-style: disc;
    }
    OL LI {
	list-style: decimal;
    }
    DT {
	padding-top: 	1ex;
    }
    UL.toc,UL.toc UL, UL.toc UL UL {
	font:		normal 12pt/14pt sans-serif;
	list-style:	none;
    }
    LI.section, LI.subsection {
	list-style: 	none;
	margin-left: 	0em;
	text-indent:	0em;
	padding: 	0em;
    }
    PRE {
	display: 	block;
	font-family:	monospace;
	white-space:	pre;
	margin:		0%;
	padding-top:	0.5ex;
	padding-bottom:	0.5ex;
	padding-left:	1ex;
	padding-right:	1ex;
	width:		100%;
    }
    PRE.example {
	color: 		black;
	background: 	#f5dcb3;
	border:		1px solid black;
    }
    UL.requirements LI, UL.syntax LI {
	list-style: 	none;
	margin-left: 	0em;
	text-indent:	0em;
	padding:	0em;
    }
    DIV.synopsis {
	color: 		black;
	background: 	#80ffff;
	border:		1px solid black;
	font-family:	serif;
	margin-top: 	1em;
	margin-bottom: 	1em;
    }
    UL.syntax {
	margin-top: 	1em;
	border-top:	1px solid black;
    }
    UL.requirements {
	margin-bottom: 	1em;
	border-bottom:	1px solid black;
    }
--></style>
</head>
<! -- Generated from file '/net/nas/data/andreask/Dev/Tcllib/tcllib/embedded/www/tcllib/files/modules/math/calculus.html' by tcllib/doctools with format 'html'
   -->
<! -- Copyright &copy; 2002,2003,2004 Arjen Markus
   -->
<! -- CVS: $Id$ math::calculus.n
   -->
<body><div class="doctools">
<hr> [
  <a href="../../../../toc.html">Main Table Of Contents</a>
| <a href="../../../toc.html">Table Of Contents</a>
| <a href="../../../../index.html">Keyword Index</a>
| <a href="/tcllib">Home</a>
] <hr>
<h1 class="title">math::calculus(n) 0.7.1 tcllib &quot;Tcl Math Library&quot;</h1>
<div id="name" class="section"><h2><a name="name">Name</a></h2>
<p>math::calculus - Integration and ordinary differential equations</p>
</div>
<div id="toc" class="section"><h2><a name="toc">Table Of Contents</a></h2>
<ul class="toc">
<li class="section"><a href="#toc">Table Of Contents</a></li>
<li class="section"><a href="#synopsis">Synopsis</a></li>
<li class="section"><a href="#section1">Description</a></li>
<li class="section"><a href="#section2">PROCEDURES</a></li>
<li class="section"><a href="#section3">EXAMPLES</a></li>
<li class="section"><a href="#section4">BUGS, IDEAS, FEEDBACK</a></li>
<li class="section"><a href="#see-also">See Also</a></li>
<li class="section"><a href="#keywords">Keywords</a></li>
<li class="section"><a href="#category">Category</a></li>
<li class="section"><a href="#copyright">Copyright</a></li>
</ul>
</div>
<div id="synopsis" class="section"><h2><a name="synopsis">Synopsis</a></h2>
<div class="synopsis">
<ul class="requirements">
<li>package require <b class="pkgname">Tcl 8.4</b></li>
<li>package require <b class="pkgname">math::calculus 0.7.1</b></li>
</ul>
<ul class="syntax">
<li><a href="#1"><b class="cmd">::math::calculus::integral</b> <i class="arg">begin</i> <i class="arg">end</i> <i class="arg">nosteps</i> <i class="arg">func</i></a></li>
<li><a href="#2"><b class="cmd">::math::calculus::integralExpr</b> <i class="arg">begin</i> <i class="arg">end</i> <i class="arg">nosteps</i> <i class="arg">expression</i></a></li>
<li><a href="#3"><b class="cmd">::math::calculus::integral2D</b> <i class="arg">xinterval</i> <i class="arg">yinterval</i> <i class="arg">func</i></a></li>
<li><a href="#4"><b class="cmd">::math::calculus::integral2D_accurate</b> <i class="arg">xinterval</i> <i class="arg">yinterval</i> <i class="arg">func</i></a></li>
<li><a href="#5"><b class="cmd">::math::calculus::integral3D</b> <i class="arg">xinterval</i> <i class="arg">yinterval</i> <i class="arg">zinterval</i> <i class="arg">func</i></a></li>
<li><a href="#6"><b class="cmd">::math::calculus::integral3D_accurate</b> <i class="arg">xinterval</i> <i class="arg">yinterval</i> <i class="arg">zinterval</i> <i class="arg">func</i></a></li>
<li><a href="#7"><b class="cmd">::math::calculus::eulerStep</b> <i class="arg">t</i> <i class="arg">tstep</i> <i class="arg">xvec</i> <i class="arg">func</i></a></li>
<li><a href="#8"><b class="cmd">::math::calculus::heunStep</b> <i class="arg">t</i> <i class="arg">tstep</i> <i class="arg">xvec</i> <i class="arg">func</i></a></li>
<li><a href="#9"><b class="cmd">::math::calculus::rungeKuttaStep</b> <i class="arg">t</i> <i class="arg">tstep</i> <i class="arg">xvec</i> <i class="arg">func</i></a></li>
<li><a href="#10"><b class="cmd">::math::calculus::boundaryValueSecondOrder</b> <i class="arg">coeff_func</i> <i class="arg">force_func</i> <i class="arg">leftbnd</i> <i class="arg">rightbnd</i> <i class="arg">nostep</i></a></li>
<li><a href="#11"><b class="cmd">::math::calculus::solveTriDiagonal</b> <i class="arg">acoeff</i> <i class="arg">bcoeff</i> <i class="arg">ccoeff</i> <i class="arg">dvalue</i></a></li>
<li><a href="#12"><b class="cmd">::math::calculus::newtonRaphson</b> <i class="arg">func</i> <i class="arg">deriv</i> <i class="arg">initval</i></a></li>
<li><a href="#13"><b class="cmd">::math::calculus::newtonRaphsonParameters</b> <i class="arg">maxiter</i> <i class="arg">tolerance</i></a></li>
<li><a href="#14"><b class="cmd">::math::calculus::regula_falsi</b> <i class="arg">f</i> <i class="arg">xb</i> <i class="arg">xe</i> <i class="arg">eps</i></a></li>
</ul>
</div>
</div>
<div id="section1" class="section"><h2><a name="section1">Description</a></h2>
<p>This package implements several simple mathematical algorithms:</p>
<ul class="itemized">
<li><p>The integration of a function over an interval</p></li>
<li><p>The numerical integration of a system of ordinary differential
equations.</p></li>
<li><p>Estimating the root(s) of an equation of one variable.</p></li>
</ul>
<p>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.</p>
<p>This document describes the procedures and explains their usage.</p>
</div>
<div id="section2" class="section"><h2><a name="section2">PROCEDURES</a></h2>
<p>This package defines the following public procedures:</p>
<dl class="definitions">
<dt><a name="1"><b class="cmd">::math::calculus::integral</b> <i class="arg">begin</i> <i class="arg">end</i> <i class="arg">nosteps</i> <i class="arg">func</i></a></dt>
<dd><p>Determine the integral of the given function using the Simpson
rule. The interval for the integration is [<i class="arg">begin</i>, <i class="arg">end</i>].
The remaining arguments are:</p>
<dl class="definitions">
<dt><i class="arg">nosteps</i></dt>
<dd><p>Number of steps in which the interval is divided.</p></dd>
<dt><i class="arg">func</i></dt>
<dd><p>Function to be integrated. It should take one single argument.</p></dd>
</dl></dd>
<dt><a name="2"><b class="cmd">::math::calculus::integralExpr</b> <i class="arg">begin</i> <i class="arg">end</i> <i class="arg">nosteps</i> <i class="arg">expression</i></a></dt>
<dd><p>Similar to the previous proc, this one determines the integral of
the given <i class="arg">expression</i> using the Simpson rule.
The interval for the integration is [<i class="arg">begin</i>, <i class="arg">end</i>].
The remaining arguments are:</p>
<dl class="definitions">
<dt><i class="arg">nosteps</i></dt>
<dd><p>Number of steps in which the interval is divided.</p></dd>
<dt><i class="arg">expression</i></dt>
<dd><p>Expression to be integrated. It should
use the variable &quot;x&quot; as the only variable (the &quot;integrate&quot;)</p></dd>
</dl></dd>
<dt><a name="3"><b class="cmd">::math::calculus::integral2D</b> <i class="arg">xinterval</i> <i class="arg">yinterval</i> <i class="arg">func</i></a></dt>
<dd></dd>
<dt><a name="4"><b class="cmd">::math::calculus::integral2D_accurate</b> <i class="arg">xinterval</i> <i class="arg">yinterval</i> <i class="arg">func</i></a></dt>
<dd><p>The commands <b class="cmd">integral2D</b> and <b class="cmd">integral2D_accurate</b> 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.</p>
<p>The command <b class="cmd">integral2D</b> evaluates the function at the centre of
each rectangle, whereas the command <b class="cmd">integral2D_accurate</b> uses a
four-point quadrature formula. This results in an exact integration of
polynomials of third degree or less.</p>
<p>The function must take two arguments and return the function
value.</p></dd>
<dt><a name="5"><b class="cmd">::math::calculus::integral3D</b> <i class="arg">xinterval</i> <i class="arg">yinterval</i> <i class="arg">zinterval</i> <i class="arg">func</i></a></dt>
<dd></dd>
<dt><a name="6"><b class="cmd">::math::calculus::integral3D_accurate</b> <i class="arg">xinterval</i> <i class="arg">yinterval</i> <i class="arg">zinterval</i> <i class="arg">func</i></a></dt>
<dd><p>The commands <b class="cmd">integral3D</b> and <b class="cmd">integral3D_accurate</b> are the
three-dimensional equivalent of <b class="cmd">integral2D</b> and <b class="cmd">integral3D_accurate</b>.
The function <em>func</em> takes three arguments and is integrated over the block in
3D space given by three intervals.</p></dd>
<dt><a name="7"><b class="cmd">::math::calculus::eulerStep</b> <i class="arg">t</i> <i class="arg">tstep</i> <i class="arg">xvec</i> <i class="arg">func</i></a></dt>
<dd><p>Set a single step in the numerical integration of a system of
differential equations. The method used is Euler's.</p>
<dl class="definitions">
<dt><i class="arg">t</i></dt>
<dd><p>Value of the independent variable (typically time)
at the beginning of the step.</p></dd>
<dt><i class="arg">tstep</i></dt>
<dd><p>Step size for the independent variable.</p></dd>
<dt><i class="arg">xvec</i></dt>
<dd><p>List (vector) of dependent values</p></dd>
<dt><i class="arg">func</i></dt>
<dd><p>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 &quot;func&quot; must match).</p></dd>
</dl></dd>
<dt><a name="8"><b class="cmd">::math::calculus::heunStep</b> <i class="arg">t</i> <i class="arg">tstep</i> <i class="arg">xvec</i> <i class="arg">func</i></a></dt>
<dd><p>Set a single step in the numerical integration of a system of
differential equations. The method used is Heun's.</p>
<dl class="definitions">
<dt><i class="arg">t</i></dt>
<dd><p>Value of the independent variable (typically time)
at the beginning of the step.</p></dd>
<dt><i class="arg">tstep</i></dt>
<dd><p>Step size for the independent variable.</p></dd>
<dt><i class="arg">xvec</i></dt>
<dd><p>List (vector) of dependent values</p></dd>
<dt><i class="arg">func</i></dt>
<dd><p>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 &quot;func&quot; must match).</p></dd>
</dl></dd>
<dt><a name="9"><b class="cmd">::math::calculus::rungeKuttaStep</b> <i class="arg">t</i> <i class="arg">tstep</i> <i class="arg">xvec</i> <i class="arg">func</i></a></dt>
<dd><p>Set a single step in the numerical integration of a system of
differential equations. The method used is Runge-Kutta 4th
order.</p>
<dl class="definitions">
<dt><i class="arg">t</i></dt>
<dd><p>Value of the independent variable (typically time)
at the beginning of the step.</p></dd>
<dt><i class="arg">tstep</i></dt>
<dd><p>Step size for the independent variable.</p></dd>
<dt><i class="arg">xvec</i></dt>
<dd><p>List (vector) of dependent values</p></dd>
<dt><i class="arg">func</i></dt>
<dd><p>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 &quot;func&quot; must match).</p></dd>
</dl></dd>
<dt><a name="10"><b class="cmd">::math::calculus::boundaryValueSecondOrder</b> <i class="arg">coeff_func</i> <i class="arg">force_func</i> <i class="arg">leftbnd</i> <i class="arg">rightbnd</i> <i class="arg">nostep</i></a></dt>
<dd><p>Solve a second order linear differential equation with boundary
values at two sides. The equation has to be of the form (the
&quot;conservative&quot; form):</p>
<pre class="example">
         d      dy     d
         -- A(x)--  +  -- B(x)y + C(x)y  =  D(x)
         dx     dx     dx
</pre>
<p>Ordinarily, such an equation would be written as:</p>
<pre class="example">
             d2y        dy
         a(x)---  + b(x)-- + c(x) y  =  D(x)
             dx2        dx
</pre>
<p>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:</p>
<pre class="example">
         A(x)  =  a(x)
         B(x)  =  b(x) - a'(x)
         C(x)  =  c(x) - B'(x)  =  c(x) - b'(x) + a''(x)
</pre>
<p>Because of the differentiation, however, it is much easier to ask
the user to provide the functions A, B and C directly.</p>
<dl class="definitions">
<dt><i class="arg">coeff_func</i></dt>
<dd><p>Procedure returning the three coefficients
(A, B, C) of the equation, taking as its one argument the x-coordinate.</p></dd>
<dt><i class="arg">force_func</i></dt>
<dd><p>Procedure returning the right-hand side
(D) as a function of the x-coordinate.</p></dd>
<dt><i class="arg">leftbnd</i></dt>
<dd><p>A list of two values: the x-coordinate of the
left boundary and the value at that boundary.</p></dd>
<dt><i class="arg">rightbnd</i></dt>
<dd><p>A list of two values: the x-coordinate of the
right boundary and the value at that boundary.</p></dd>
<dt><i class="arg">nostep</i></dt>
<dd><p>Number of steps by which to discretise the
interval.
The procedure returns a list of x-coordinates and the approximated
values of the solution.</p></dd>
</dl></dd>
<dt><a name="11"><b class="cmd">::math::calculus::solveTriDiagonal</b> <i class="arg">acoeff</i> <i class="arg">bcoeff</i> <i class="arg">ccoeff</i> <i class="arg">dvalue</i></a></dt>
<dd><p>Solve a system of linear equations Ax = b with A a tridiagonal
matrix. Returns the solution as a list.</p>
<dl class="definitions">
<dt><i class="arg">acoeff</i></dt>
<dd><p>List of values on the lower diagonal</p></dd>
<dt><i class="arg">bcoeff</i></dt>
<dd><p>List of values on the main diagonal</p></dd>
<dt><i class="arg">ccoeff</i></dt>
<dd><p>List of values on the upper diagonal</p></dd>
<dt><i class="arg">dvalue</i></dt>
<dd><p>List of values on the righthand-side</p></dd>
</dl></dd>
<dt><a name="12"><b class="cmd">::math::calculus::newtonRaphson</b> <i class="arg">func</i> <i class="arg">deriv</i> <i class="arg">initval</i></a></dt>
<dd><p>Determine the root of an equation given by</p>
<pre class="example">
    func(x) = 0
</pre>
<p>using the method of Newton-Raphson. The procedure takes the following
arguments:</p>
<dl class="definitions">
<dt><i class="arg">func</i></dt>
<dd><p>Procedure that returns the value the function at x</p></dd>
<dt><i class="arg">deriv</i></dt>
<dd><p>Procedure that returns the derivative of the function at x</p></dd>
<dt><i class="arg">initval</i></dt>
<dd><p>Initial value for x</p></dd>
</dl></dd>
<dt><a name="13"><b class="cmd">::math::calculus::newtonRaphsonParameters</b> <i class="arg">maxiter</i> <i class="arg">tolerance</i></a></dt>
<dd><p>Set the numerical parameters for the Newton-Raphson method:</p>
<dl class="definitions">
<dt><i class="arg">maxiter</i></dt>
<dd><p>Maximum number of iteration steps (defaults to 20)</p></dd>
<dt><i class="arg">tolerance</i></dt>
<dd><p>Relative precision (defaults to 0.001)</p></dd>
</dl></dd>
<dt><a name="14"><b class="cmd">::math::calculus::regula_falsi</b> <i class="arg">f</i> <i class="arg">xb</i> <i class="arg">xe</i> <i class="arg">eps</i></a></dt>
<dd><p>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.</p>
<p>The method used is the so-called <em>regula falsi</em> or
<em>false position</em> 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.</p>
<p>In contrast to Newton-Raphson there is no need for the computation of
the function's derivative.</p>
<dl class="arguments">
<dt>command <i class="arg">f</i></dt>
<dd><p>Name of the command that evaluates the function for
which the zero is to be returned</p></dd>
<dt>float <i class="arg">xb</i></dt>
<dd><p>Start of the interval in which the zero is supposed
to lie</p></dd>
<dt>float <i class="arg">xe</i></dt>
<dd><p>End of the interval</p></dd>
<dt>float <i class="arg">eps</i></dt>
<dd><p>Relative allowed error (defaults to 1.0e-4)</p></dd>
</dl></dd>
</dl>
<p><em>Notes:</em></p>
<p>Several of the above procedures take the <em>names</em> of procedures as
arguments. To avoid problems with the <em>visibility</em> 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:</p>
<pre class="example">
    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]
       }
    }
</pre>
<p>Enhancements for the second-order boundary value problem:</p>
<ul class="itemized">
<li><p>Other types of boundary conditions (zero gradient, zero flux)</p></li>
<li><p>Other schematisation of the first-order term (now central
differences are used, but upstream differences might be useful too).</p></li>
</ul>
</div>
<div id="section3" class="section"><h2><a name="section3">EXAMPLES</a></h2>
<p>Let us take a few simple examples:</p>
<p>Integrate x over the interval [0,100] (20 steps):</p>
<pre class="example">
proc linear_func { x } { return $x }
puts &quot;Integral: [::math::calculus::integral 0 100 20 linear_func]&quot;
</pre>
<p>For simple functions, the alternative could be:</p>
<pre class="example">
puts &quot;Integral: [::math::calculus::integralExpr 0 100 20 {$x}]&quot;
</pre>
<p>Do not forget the braces!</p>
<p>The differential equation for a dampened oscillator:</p>
<pre class="example">
x'' + rx' + wx = 0
</pre>
<p>can be split into a system of first-order equations:</p>
<pre class="example">
x' = y
y' = -ry - wx
</pre>
<p>Then this system can be solved with code like this:</p>
<pre class="example">
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 &lt; 20 } { incr i } {
   set result [::math::calculus::eulerStep $t $tstep $xvec dampened_oscillator]
   puts &quot;Result ($t): $result&quot;
   set t      [expr {$t+$tstep}]
   set xvec   $result
}
</pre>
<p>Suppose we have the boundary value problem:</p>
<pre class="example">
    Dy'' + ky = 0
    x = 0: y = 1
    x = L: y = 0
</pre>
<p>This boundary value problem could originate from the diffusion of a
decaying substance.</p>
<p>It can be solved with the following fragment:</p>
<pre class="example">
   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]
</pre>
</div>
<div id="section4" class="section"><h2><a name="section4">BUGS, IDEAS, FEEDBACK</a></h2>
<p>This document, and the package it describes, will undoubtedly contain
bugs and other problems.
Please report such in the category <em>math :: calculus</em> of the
<a href="http://sourceforge.net/tracker/?group_id=12883">Tcllib SF Trackers</a>.
Please also report any ideas for enhancements you may have for either
package and/or documentation.</p>
</div>
<div id="see-also" class="section"><h2><a name="see-also">See Also</a></h2>
<p>romberg</p>
</div>
<div id="keywords" class="section"><h2><a name="keywords">Keywords</a></h2>
<p><a href="../../../../index.html#key722">calculus</a>, <a href="../../../../index.html#key723">differential equations</a>, <a href="../../../../index.html#key721">integration</a>, <a href="../../../../index.html#key67">math</a>, <a href="../../../../index.html#key724">roots</a></p>
</div>
<div id="category" class="section"><h2><a name="category">Category</a></h2>
<p>Mathematics</p>
</div>
<div id="copyright" class="section"><h2><a name="copyright">Copyright</a></h2>
<p>Copyright &copy; 2002,2003,2004 Arjen Markus</p>
</div>
</div></body></html>
