| 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.
Ordinarily, such an equation would be written as:
d dy d
-- A(x)-- + -- B(x)y + C(x)y = D(x)
dx dx 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:
d2y dy
a(x)--- + b(x)-- + c(x) y = D(x)
dx2 dx
Because of the differentiation, however, it is much easier to ask the user to provide the functions A, B and C directly.
A(x) = a(x)
B(x) = b(x) - a'(x)
C(x) = c(x) - B'(x) = c(x) - b'(x) + a''(x)
using the method of Newton-Raphson. The procedure takes the following arguments:
func(x) = 0
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 |