http://www.opencontent.org/openpub
-->
tclode::odesolv(n) 1.0 "Tcl ODE Solver"
tclode::odesolv - Tcl ODE Solver - Initial Value Problems
package require Tcl 8.4
package require tclode 1.0
The odesolv command in the tclode package performs numerical
integration for initial value problems (IVP's). It is a version of the
Livermore Solver for Ordinary Differential Equations with Automatic
Rootfinding (LSODAR), available from http://www.netlib.org/odepack, but
adapted to calling from Tcl rather than Fortran.
The tclode package implements a single command:
- ::tclode::odesolv ?-option value?... ?--? varName derivativeExpr ?varName derivativeExpr?...
-
Sets up to integrate an Initial Value Problem. The return value from
::tclode::odesolv is the name of a solver, an object that can be
used to find values of the variables at particular times (and perform other
tasks, according to the options presented). See THE SOLVER OBJECT
below for details of how the solver is used.
- ::tclode::odesolv ?-option value?... ?--? system
-
An alternative form to the one above. In this form, system is a list comprising variable names and derivative expressions in alternation.
- varName
-
Name of a variable in the system of ODE's to integrate. All the variable names
must be distinct.
- derivativeExpr
-
Expression that gives the derivative of a model variable. The
expression can take any form that is acceptable to the expr
command. The model variables may be included in the expression using
\$-substitution, and in fact, any other variables in the caller's scope
may also be included.
- -atol value
-
If value is a floating-point number, requests that the
integration machinery return results that are accurate to within an
absolute tolerance of value. Otherwise value must be a
list of floating point numbers, which are absolute tolerances to apply
to the model variables, in the order in which they appear in the
tclode::odesolv command.
- -findroots exprList
-
The exprList argument must be a list of expressions for
which roots are to be found. This option, if present, requests that the
solver stop whenever the value of one of the expressions changes sign.
Any expression acceptable to the expr command can appear as an
element of exprList.
- -h0 value
-
The value argument must be a floating-point number. This option
requests that the solver use the given value as the initial step size
when commencing integration.
- -hmax value
-
The value argument must be a floating-point number. This option
requests that the solver avoid taking integration steps larger than the
given value.
- -hmin value
-
The value argument must be a floating-point number. This options
requests that the solver avoid taking integration steps smaller than the
given value.
- -indvar varName
-
This option requests that the variable designated by varName be the
independent, or time, variable in the initial value problem. If this option
is not supplied, the independent variable is named, 't'.
- -ixpr level
-
This option asks the solver to return verbose data when the
integration method is changed. If the level argument is 0, then
no extra data is returned. If it is 1 or greater, then the return value
from calls to the solver may include messages indicating the choice of
integration method and the value of the independent variable when the
method is changed.
- -jacobian matrix
-
This option supplies a symbolic Jacobian matrix to the solver. The
Jacobian is the matrix of second derivatives of the equations to be solved.
In a system having n model variables, the Jacobian must be an
n-element list. Each list element is a row of the Jacobian matrix,
and is itself an n-element list of Tcl expressions. The ith
element of the jth row gives the second derivative of the jth
variable with respect to the ith variable and the independent variable.
If, for instance, the system has two model variables, 'x' and 'y', then
the Jacobian should comprise two lists of two expressions each. In order,
four expressions
should be: d/dx(dx/dt), d/dy(dx/dt), d/dx(dy/dt) and d/dy(dy/dt).
If this option is not supplied, the solver will first attempt to use the
symdiff package, if available, to differentiate the model derivatives
with respect to the model variables. If that fails, the solver will resort
to approximating the derivatives numerically.
- -mxhnil value
-
This option controls the maximum number of messages that will be returned from
a single call to the solver. The default is 10. The value must be a
non-negative integer.
- -mxordn value
-
This option controls the maximum order to be allowed when using
Adams-Moulton integration for non-stiff systems of equations.
The value must be an integer in the range [1..12].
The default is 12.
- -mxords value
-
This option controls the maximum order to be allowed when using
Gear backward-differencing formulas for stiff systems of equations.
The value must be an integer in the range [1..5].
The default is 5.
- -mxstep value
-
This option controls the maximum number of integration steps to be
taken in one call to the solver. The value must be a positive
integer. The default is 500.
- -roots varName
-
This option asks the solver to return indications of what roots
were found when the -findroots option was supplied. The varName
argument is the name of a variable in the calling scope. After
each call to the solver, this variable will be set to a list of
indices in the list of expressions presented to -findroots
corresponding to expressions whose values changed sign on the
current integration step.
- -rtol value
-
If value is a floating-point number, requests that the
integration machinery return results that are accurate to within a
relative tolerance of value. Otherwise value must be a
list of floating point numbers, which are relative tolerances to apply
to the model variables, in the order in which they appear in the
tclode::odesolv command.
The return value from tclode::odesolv is itself a Tcl command,
that can be used to manipulate the given solver. The solver command
accepts subcommands:
- $solver run t
-
Requests that the given system of equations be run until the independent
variable has advanced to t. The run may be terminated prematurely
if the -findroots option is specified and a root is found,
if the maximum number of steps given by the -mxstep option is
exceeded, or if an error occurs. The return value from the run command
is a human-readable set of messages describing what happened in the current
run.
- $solver continue t
-
Requests that the integration continue from where it left off on the
last call to run or continue.
The integration is to proceed until the independent
variable has advanced to t. The run may be terminated prematurely
if the -findroots option is specified and a root is found,
if the maximum number of steps given by the -mxstep option is
exceeded, or if an error occurs. The return value from the run command
is a human-readable set of messages describing what happened in the current
run.
- $solver status
-
Requests machine-readable information about what happened on the last call
to run or continue. The information is returned as a dictionary.
See THE STATUS DICTIONARY below for the meanings of the
values.
Prior to starting a run of the solver, the program must initialize all
model variables to their initial values, and initialize the
independent variable to its starting point. The program should not change
the values of the model variables or independent variable before a call
to continue, but may do so before a call to run.
When the run or continue commands return successfully, the
independent variable is set to the point at which integration stopped,
and the model variables are left set to the values that they took at that
point.
The status command when given to a solver returns a dictionary
that gives the results of the last integration in macine-readable form.
The keys of the dictionary are as follows:
- istate
-
Value is an integer that describes the reason for success or failure of
the integration. See ISTATE VALUES below for a table of possible values.
- hu
-
Value is a floating-point number giving the step size used on the last
successful integration step.
- hcur
-
Value is a floating-point number giving the step size that will be attempted
on the next integration step, if the continue command is used.
- tcur
-
Value is a floating-point number giving the value of the independent variable
that the solver has actually reached. This value will usually be farther than
the "current" point, which is often found by interpolation over the mesh
points.
- tolsf
-
If the integration has failed because too much precision was requested in
the -atol or -rtol options, this value will be a scale
factor. If the tolerances presented in -atol or -rtol are
uniformly increased by a factor of 'tolsf', then the integration is likely
to succeed.
- tsw
-
Value is a floating-point number giving the value of the independent variable
at the last time that the solver switched between integration methods.
- nge
-
Value is an integer giving the number of times that the expressions
for which roots are to be found have been evaluated.
- nst
-
Value is an integer giving the number of integration steps that have been
taken so far.
- nfe
-
Value is an integer giving the number of times that the derivative expressions
have been evaluated so far.
- nje
-
Value is an integer giving the number of times that the Jacobian matrix
has been evaluated so far.
- nqu
-
Value is an integer giving the order of the last successful integration step.
- nqcur
-
Value is an integer giving the order that will be tried in the next
integration step if the continue command is used.
- imxer
-
Value is an integer giving the index in the set of model variables
corresponding to the component of largest magnitude in the weighted
local error vector. This can give information if the integration fails
to converge.
- mused
-
Value is 1 if the last successful integration step used a non-stiff
(Adams-Moulton) solver, or 2 if it used a stiff (Gear backward-difference
formula) solver.
- mcur
-
Value is 1 if the next integration step will attempt to use a non-stiff
(Adams-Moulton) solver, or 2 if it will use a stiff (Gear backward-difference
formula) solver.
In the dictionary returned from the status command on a solver,
the 'istate' value is an integer that gives the reason the last integration
step halted. Positive values indicate success; negative values indicate
some sort of trouble. The possible values, and there meanings, are:
- 1
-
Nothing was done. The current value of the independent variable is already
at or beyond the target point.
- 2
-
Integration was performed successfully and stopped at the target value of the
independent variable.
- 3
-
Integration was performed successfully and stopped because a root was
found for one of the expressions given in the -findroots option.
- -1
-
Integration stopped because an excessive number of steps were requested.
It was otherwise successful up to the current value of the independent
variable. It may be continued using the continue command to the solver.
- -2
-
Integration stopped because too much accuracy was requested. It may
be continued if the -rtol or -atol values are reset.
The 'tolsf' member of the status dictionary gives a suggestion regarding
what scale factor might be appropriate.
- -3
-
Integration stopped because the Fortran code detected an illegal combination
of parameters. The return value from the run or continue
command should contain an error message with further information about
what is wrong.
- -4
-
Integration stopped because there were repeated error test failures on
one attempted step. It was successful as far as the current value of the
independent variable. The most likely cause of this status code is a
singularity in the system of equations being integrated.
- -5
-
Integration stopped because of repeated convergence test failures.
The commonest cause of failure to converge is a Jacobian matrix that is
either inaccurate or outright incorrect.
- -6
-
The -rtol option was specified for some model variable, and the
-atol option was not specified. The model variable has vanished
at the given value of the independent variable, so the error control can
not be evaluated. The integration was successful as far as the current
value of the independent variable.
- -7
-
If this state occurs, it indicates an error in the internal memory management
of the tclode package itself. Please report a bug to the maintainers.
The tclode::solver command is a complicated command, and
it is difficult to give a full spectrum of the use cases. A few simple ones
follow.
| |
set solver [tclode::odesolv -atol 1.0e-6 -rtol 1.0e-6 -findroots {{$s - 0.5} {$c - sqrt(0.5)}} -roots which -indvar theta -- s {$c} c {-$s}
set c 1.0
set s 0.0
set theta 0.0
$solver run 3.14159
rename $solver {}
|
This example sets up the trivial set of differential equations
ds/dt=c; dc/dt=-s. The initial conditions are c=1 and s=0.
The solver is asked to integrate until either the independent
variable theta reaches 3.14159, or the value of c reaches
sqrt(0.5), or the value of s reaches 0.5.
At the conclusion of this block of code, the value of 'which' is 0, indicating
that '$s-0.5' (the expression at index 0 in the ?-findroots? argument)
reached zero. The value of 's' is very nearly 0.5, that of 'c' is very close
to sqrt(3)/2, and that of 'theta' is very close to 'pi/6'.
| |
set solver [tclode::odesolv -rtol 1.0e-6 -atol 1.0e-6 -ixpr 1 -jacobian {{ 998. 1998.}
{-999. -1999.}} u {998. * $u + 1998. * $v} v {-999. * $u - 1999. * $v}]
set t 0.0
set u 1.0
set v 0.0
for {set tout 0.05} {$tout < 1.000} {set tout [expr {$tout + 0.05}]} {
set msg [$solver continue $tout]
if {$msg ne {}} {
puts $msg
}
puts [list $t $u $v]
}
|
This example sets up a classical "stiff" system of equations: du/dt=998u+1998v
and dv/dt=-999u-1999v, with initial values of u=1 and v=0. It runs the
solver through closely spaced values of t, printing a table of results.
It also prints the message when the solver detects stiffness and switches
from Adams-Moulton integration to Gear backward differencing.
Copyright © 2007 Kevin B. Kenny <kennykb@acm.org>. Redistribution permitted under the terms of the Open Publication License http://www.opencontent.org/openpub