http://www.opencontent.org/openpub -->

tclode::odesolv(n) 1.0 "Tcl ODE Solver"

NAME

tclode::odesolv - Tcl ODE Solver - Initial Value Problems

SYNOPSIS

package require Tcl 8.4
package require tclode 1.0

::tclode::odesolv ?-option value?... ?--? varName derivativeExpr ?varName derivativeExpr?...
::tclode::odesolv ?-option value?... ?--? system
$solver run t
$solver continue t
$solver status

DESCRIPTION

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.

PROCEDURES

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.

PARAMETERS

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.

OPTIONS

-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 SOLVER 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 DICTIONARY

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.

ISTATE VALUES

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.

EXAMPLES

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

Copyright © 2007 Kevin B. Kenny <kennykb@acm.org>. Redistribution permitted under the terms of the Open Publication License http://www.opencontent.org/openpub