This package defines the following public procedures:
::math::calculus::integralbeginendnostepsfunc
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::integralExprbeginendnostepsexpression
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::integral2Dxintervalyintervalfunc::math::calculus::integral2D_accuratexintervalyintervalfunc
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::integral3Dxintervalyintervalzintervalfunc::math::calculus::integral3D_accuratexintervalyintervalzintervalfunc
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::qk15xstartxendfuncnosteps
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_detailedxstartxendfuncnosteps
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::eulerStepttstepxvecfunc
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::heunStepttstepxvecfunc
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::rungeKuttaStepttstepxvecfunc
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::boundaryValueSecondOrdercoeff_funcforce_funcleftbndrightbndnostep
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::solveTriDiagonalacoeffbcoeffccoeffdvalue
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::newtonRaphsonfuncderivinitval
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::newtonRaphsonParametersmaxitertolerance
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_falsifxbxeeps
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 regulafalsi or falseposition 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)
::math::calculus::root_bisectionfxbxeeps
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 bisection. For properties: see the regulafalsi procedure.
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-7)
::math::calculus::root_secantfxbxeeps
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 secant method. For properties: see the regulafalsi procedure.
Note that this method is not guaranteed to produce a zero, but it is fast.
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-7)
::math::calculus::root_brentfxbxeeps
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 Brent method. For properties: see the regulafalsi procedure.
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-7)
::math::calculus::root_chandrupatlafxbxeeps
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 Chandrupatla method. It is a variant of the Brent method and is
reputed to be faster. See https://www.embeddedrelated.com/showarticle/855.php for a discussion.
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-7)
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).
• The various root finding methods differ in robustness, ease of use and convergence rates. While
the Newton-Raphson method is quite fast (quadratic convergence), it is not guaranteed to produce
an answer. The methods that bracket the root, require two starting points, but except for the
secant method are guaranteed to deliver a good estimate.
Install Now
PNG Icons
