<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 © 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 "Tcl Math Library"</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 "x" as the only variable (the "integrate")</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 "func" 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 "func" 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 "func" 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
"conservative" 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 "Integral: [::math::calculus::integral 0 100 20 linear_func]"
</pre>
<p>For simple functions, the alternative could be:</p>
<pre class="example">
puts "Integral: [::math::calculus::integralExpr 0 100 20 {$x}]"
</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 < 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
}
</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 © 2002,2003,2004 Arjen Markus</p>
</div>
</div></body></html>
