## Friday, August 12, 2011

### ODEPACK - Ordinary Differential Equation Solvers

Some information on ODEPACK.

ODEPACK - Ordinary Differential Equation Solvers: "ODEPACK
Ordinary Differential Equation Solvers"

ODEPACK is a FORTRAN77 library which implements Alan Hindmarsh's solvers for ordinary differential equations.

The library includes routines commonly referred to as

LSODE solves nonstiff or stiff systems y' = f(y,t);

LSODES is like LSODE, but in the stiff case the Jacobian matrix is assumed to be sparse, and treated with sparse routines;

LSODA automatically switches between nonstiff and stiff solvers depending on the behavior of the problem;

LSODAR is like LSODE, but includes the ability to detect and solve for solutions of a related set of algebraic equations;

LSODPK is like LSODE, but uses preconditioned Krylov space iterative methods for the linear equation solvers;

LSODKR includes the rootfinding ability of LSODA, and the Krylov solvers of LSODPK.

LSODI solves the implicit system A(y,t)*y'(t) = g(y,t).

LSOIBT is like LSODI, but assumes the matrix A is block tridiagonal.

LSODIS is like LSODI, but assumes the matrix A is sparse.

ODEPACK is a collection of Fortran solvers for the initial value problem for ordinary differential equation systems. It consists of nine solvers, namely a basic solver called LSODE and eight variants of it -- LSODES, LSODA, LSODAR, LSODPK, LSODKR, LSODI, LSOIBT, and LSODIS. The collection is suitable for both stiff and nonstiff systems. It includes solvers for systems given in explicit form, dy/dt = f(t,y), and also solvers for systems given in linearly implicit form, A(t,y) dy/dt = g(t,y). Two of the solvers use general sparse matrix solvers for the linear systems that arise. Two others use iterative (preconditioned Krylov) methods instead of direct methods for these linear systems. The most recent addition is LSODIS, which solves implicit problems with general sparse treatment of all matrices involved.

The ODEPACK solvers are written in standard Fortran 77, with a few exceptions, and with minimal machine dependencies. There are separate double and single precision versions of ODEPACK. The actual solver names are those given above with a prefix of D- or S- for the double or single precision version, respectively, i.e. DLSODE/SLSODE, etc. Each solver consists of a main driver subroutine having the same name as the solver and some number of subordinate routines. For each solver, there is also a demonstration program, which solves one or two simple problems in a somewhat self-checking manner.

Recently, the ODEPACK solvers were upgraded to improve their portability in numerous ways. Among the improvements are (a) renaming of routines and Common blocks to distinguish double and single precision versions, (b) use of generic intrinsic function names, (c) elimination of the Block Data subprogram, (d) use of a portable routine to set the unit roundoff, and (e) passing of quoted strings to the error message handler. In addition, the prologue and internal comments were reformatted, and use mixed upper/lower case. Numerous minor corrections and improvements were also made.

The above upgrade operations were applied to LSODE earlier than they were to the rest of ODEPACK, and the two upgrades were done somewhat independently. As a result, some differences will be apparent in the source files of LSODE and the other solvers -- primarily in the formatting of the comment line prologue of the main driver routine. In Subroutines DLSODE/SLSODE and their subordinate routines, the prologue was written in "SLATEC format", while for the other solvers a more relaxed style was used. The differences are entirely cosmetic, however, and do not affect performance.

Documentation on the usage of each solver is provided in the initial block of comment lines in the source file, which (in most cases) includes a simple example.

Solvers for explicit systems:

For each of the following solvers, it is assumed that the ODEs are given explicitly, so that the system can be written in the form dy/dt = f(t,y), where y is the vector of dependent variables, and t is the independent variable.

LSODE (Livermore Solver for Ordinary Differential Equations) is the basic solver of the collection. It solves stiff and nonstiff systems of the form dy/dt = f. In the stiff case, it treats the Jacobian matrix df/dy as either a dense (full) or a banded matrix, and as either user-supplied or internally approximated by difference quotients. It uses Adams methods (predictor-corrector) in the nonstiff case, and Backward Differentiation Formula (BDF) methods (the Gear methods) in the stiff case. The linear systems that arise are solved by direct methods (LU factor/solve). LSODE supersedes the older GEAR and GEARB packages, and reflects a complete redesign of the user interface and internal organization, with some algorithmic improvements.

LSODES, written jointly with A. H. Sherman, solves systems dy/dt = f and in the stiff case treats the Jacobian matrix in general sparse form. It determines the sparsity structure on its own, or optionally accepts this information from the user. It then uses parts of the Yale Sparse Matrix Package (YSMP) to solve the linear systems that arise, by a sparse (direct) LU factorization/backsolve method. LSODES supersedes, and improves upon, the older GEARS package.

LSODA, written jointly with L. R. Petzold, solves systems dy/dt = f with a dense or banded Jacobian when the problem is stiff, but it automatically selects between nonstiff (Adams) and stiff (BDF) methods. It uses the nonstiff method initially, and dynamically monitors data in order to decide which method to use.

LSODAR, also written jointly with L. R. Petzold, is a variant of LSODA with a rootfinding capability added. Thus it solves problems dy/dt = f with dense or banded Jacobian and automatic method selection, and at the same time, it finds the roots of any of a set of given functions of the form g(t,y). This is often useful for finding stop conditions, or for finding points at which a switch is to be made in the function f.

LSODPK, written jointly with Peter N. Brown, is a variant of LSODE in which the direct solvers for the linear systems have been replaced by a selection of four preconditioned Krylov (iterative) solvers. The user must supply a pair of routine to evaluate, preprocess, and solve the (left and/or right) preconditioner matrices. LSODPK also includes an option for a user-supplied linear system solver to be used without Krylov iteration.

LSODKR is a variant of LSODPK with the addition of the same rootfinding capability as in LSODAR, and also of automatic switching between functional and Newton iteration. The nonlinear iteration method-switching differs from the method-switching in LSODA and LSODAR, but provides similar savings by using the cheaper method in the non-stiff regions of the problem. LSODKR also improves on the Krylov methods in LSODPK by offering the option to save and reuse the approximate Jacobian data underlying the preconditioner.

Solvers for linearly implicit systems:

The following solvers treat systems in the linearly implicit form A(t,y) dy/dt = g(t,y), A = a square matrix, i.e. with the derivative dy/dt implicit, but linearly so. These solvers allow A to be singular, in which case the system is a differential-algebraic equation (DAE) system. In that case, the user must be very careful to supply a well-posed problem with consistent initial conditions.

LSODI, written jointly with J. F. Painter, solves linearly implicit systems in which the matrices involved (A, dg/dy, and d(A dy/dt)/dy) are all assumed to be either dense or banded. LSODI supersedes the older GEARIB solver and improves upon it in numerous ways.

LSOIBT, written jointly with C. S. Kenney, solves linearly implicit systems in which the matrices involved are all assumed to be block-tridiagonal. Linear systems are solved by the LU method.

LSODIS, written jointly with S. Balsdon, solves linearly implicit systems in which the matrices involved are all assumed to be sparse. Like LSODES, LSODIS either determines the sparsity structure or accepts it from the user, and uses parts of the Yale Sparse Matrix Package to solve the linear systems that arise, by a direct method.