| math::calculus(n) | Tcl Math Library | math::calculus(n) |
math::calculus - Integration and ordinary differential equations
package require Tcl 8.4
package require math::calculus 0.7.1
::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::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
This package implements several simple mathematical algorithms:
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.
This package defines the following public procedures:
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.
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.
func(x) = 0
using the method of Newton-Raphson. The procedure takes the following
arguments: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.
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:
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]
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 SF Trackers [http://sourceforge.net/tracker/?group_id=12883]. Please also report any ideas for enhancements you may have for either package and/or documentation.
romberg
calculus, differential equations, integration, math, roots
Mathematics
Copyright (c) 2002,2003,2004 Arjen Markus
| 0.7.1 | math |