Sundials equation solvers
| Alan C. Hindmarsh and Radu Serban (2007), Scholarpedia, 2(3):2860. | doi:10.4249/scholarpedia.2860 | revision #90765 [link to/cite this article] |
SUNDIALS (SUite of Nonlinear and DIfferential/ALgebraic equation Solvers) is a suite of advanced computational codes for solving large-scale problems that can be modeled as a system of nonlinear algebraic equations, or as initial-value problems in ordinary differential or differential-algebraic equations.
Contents |
Introduction
The SUNDIALS solvers address three problem types:
- Ordinary Differential Equation (ODE) systems (initial value problems),
- Nonlinear algebraic systems of equations, and
- Differential-algebraic equation (DAE) systems (initial value problems).
The solvers are written in ANSI C. For each problem type, there is a basic solver. For ODE and DAE systems, there are extensions that perform sensitivity analysis. All five solvers run in either a serial or a parallel machine setting, by virtue of the user's choice of the module for vector operations. The implicit methods used lead to linear systems, for which a variety of direct and Krylov iterative methods are available. The three basic solvers include support for Fortran applications; Matlab interfaces are also provided for the ODE and DAE solvers with sensitivity capabilities, as well as for the nonlinear solver. Full details (except for recent modifications) are available in the paper Hindmarsh et al. (2005).
SUNDIALS was implemented with the goal of providing robust solvers that can easily be incorporated into existing simulation codes. The primary design goals were to require minimal information from the user, allow users to easily supply their own data structures underneath the solvers, and allow for easy incorporation of user-supplied linear solvers and preconditioners.
CVODE, for ODE Systems
CVODE solves ODE initial value problems, in real \(N\)-space, written as \[ \dot{y} = f(t,y) , ~~ y(t_0) = y_0 ~~ (y \in R^N) . \] The user first selects one of two variable-order, variable-step linear multistep method families — implicit Adams methods (orders 1 to 12) or methods based on Backward Differentiation Formulas (BDFs) (orders 1 to 5). Then the user specifies either functional or Newton iteration for the treatment of the implicit nonlinear equations. For nonstiff systems, Adams method with functional iteration is sufficient. For stiff systems, characterized by at least one rapid decay mode (with time constant much smaller than the solution time scale), one must choose Newton iteration and a linear system solver that is appropriate to the problem. The user may also require that CVODE find, and stop at, the roots of a set of given functions during the integration of the ODEs.
In the Newton case, CVODE must solve linear systems of dimension \(N \times N\) of the form \((I - \gamma J) x = b \) that arise at each time step. Here \(\gamma\) is a scalar and \(J\) is the Jacobian \(\partial f / \partial y .\) The user specifies one of six algorithms to solve these systems.
The linear system solvers available are as follows. The first three are direct methods, and the last three are Krylov iterative methods.
- a dense direct solver (serial version only),
- a band direct solver (serial version only),
- a diagonal approximate Jacobian solver,
- Generalized Minimal Residual (GMRES) iteration,
- Bi-Conjugate Gradient Stable iteration (BiCGStab), and
- Transpose-Free Quasi-Minimal Residual (TFQMR) iteration.
For the dense and band solvers, the Jacobian matrix can be either supplied by the user or internally approximated by difference quotients. In the direct cases, the nonlinear iteration at each time step is Modified Newton, and the approximate Jacobian used is only updated when necessary to achieve convergence, rather than every step. The Krylov iterative methods all include scaling and (optionally) preconditioning (on the left, on the right, or on both sides). Here the nonlinear iteration is Inexact Newton iteration, and again the preconditioner is updated only when necessary.
CVODE chooses its step sizes and method orders automatically and dynamically, so as to keep estimated integration errors within given tolerances. Both relative and absolute tolerance parameters are required inputs from the user. Optionally, the user can specify that the selection of method order is augmented by an algorithm that attempts to detect when step sizes are limited by the BDF stability region boundary at order 3 or more.
KINSOL, for Nonlinear Algebraic Systems
KINSOL solves algebraic systems in real N-space, written as \[ F(u) = 0, ~~ F: R^N \rightarrow R^N, \] given an initial guess \(u_0 .\) The basic method is either a modified or an inexact Newton iteration. The linear systems that arise, \(J \delta u = -F(u)\) where \(J = \partial F / \partial u ,\) are solved with either a direct (dense or banded) solver (serial version only), or one of the Krylov iterative solvers (GMRES, BiCGStab, or TFQMR) listed above. In the Krylov case, the user can (optionally) supply a right preconditioner.
The user can specify that only full Newton corrections \(\delta u\) are to be used, or activate a line search globalization strategy. In the latter, the step vector is \(\lambda \cdot \delta u ,\) with \(\lambda \leq 1\) chosen by a set of criteria aimed at guaranteeing a decrease in the norm of \(F .\)
The user can specify a scaling on the \(u\) vector and/or a scaling on the \(F\) vector (for use in norms), and optional tolerances on the norm of the step and on the norm of \(F(u) ,\) for use in the stopping tests. There are also three choices for the stopping tests in the Krylov iterations. KINSOL allows the user to enforce inequality constraints on each component of \(u\) — either \( u^i > 0, u^i < 0, u^i \geq 0, u^i \leq 0, \) or no constraint.
IDA, for DAE Systems
IDA solves real differential-algebraic systems in \(N\)-space, in the general form \[ F(t,y,\dot{y}) = 0 , ~~ y(t_0) = y_0 , ~~ \dot{y}(t_0) = \dot{y}_0 . \] If the supplied initial conditions are not consistent with the equations, IDA has an option to correct them, for a class of problems (including semi-explicit index-one systems). The integration of the system uses variable-step variable-order BDF methods, up to order 5.
At each step, a Newton iteration leads to linear systems \(J x = b ,\) in which \(J = \partial F / \partial y + \alpha \partial F / \partial \dot{y}\) and \(\alpha\) is a scalar. These are solved by one of five methods — two direct (dense or band; serial version only) and three Krylov (GMRES, BiCGStab, or TFQMR). The direct methods allow either a user-supplied Jacobian or an internal approximation, and the Krylov methods include scaling and optional left preconditioning. The nonlinear iteration is Modified or Inexact Newton, and Jacobian information is updated only when necessary for convergence, not at every step.
As with CVODE, IDA chooses step sizes and orders dynamically to control local errors according to user tolerances (relative and absolute), and the user can have the integration stop at roots of given functions. As with KINSOL, the user may specify inequality constraints on the components of \(y .\)
Preconditioning
The Krylov linear system methods (GMRES, BiCGStab, and TFQMR) are "matrix-free" as they stand. This means that, for the solution of a system \(A x = b \) each iteration requires only the value of a matrix-vector product \(A v ,\) and that product is (on default) obtained by a difference quotient not requiring the matrix \(A\) explicitly.
To be effective, however, Krylov methods generally require preconditioning. This involves defining a matrix \(P ,\) and applying the method to the equivalent system \((P^{-1} A) x = P^{-1} b \) (left preconditioning) or \((A P^{-1})(P x) = b \) (right preconditioning). The preconditioner matrix \(P\) should capture the dominant parts of the system Jacobian, but at the same time permit a reasonably economical solution of systems \(P x = b .\) All SUNDIALS solvers allow the user to supply preconditioning. When this is done, the user supplies two separate functions — one to construct and preprocess \(P ,\) and another to solve systems \(P x = b .\) The reason for this separation into two phases is that the first operation is done much less frequently than the second.
Alternatively, the user can call on one of two preconditioner modules included in SUNDIALS. The serial versions of CVODE and CVODES offer a preconditioner that generates banded difference-quotient approximations to the Jacobian. Also, in a parallel environment, all the solvers offer a block-diagonal preconditioner with banded blocks (one per processor) generated by difference quotients.
Sensitivity Analysis
When the problem-defining function (\(f\) or \(F\)), and possibly also the initial value vector(s), involve a set of parameters \(p ,\) it is often desirable to compute the derivative with respect to \(p\) of the solution \(y\) (or of some output function of \(y\)). The SUNDIALS codes CVODES (available now) and IDAS (in development) compute such sensitivities, using either of two approaches. In forward sensitivity analysis, an ODE or DAE for \(\partial y /\partial p\) is generated and solved. In adjoint sensitivity analysis, more suitable when the number of parameters is large, the system is integrated forward and the solution saved at selected checkpoints, and another auxiliary equation (the adjoint system) is generated and integrated backwards (along with partial forward integrations between checkpoints). Then, the desired sensitivities are computed in terms of the adjoint solution using a quadrature. The codes offer various options for the way in which the auxiliary equations are generated and solved along with the original state equations.
For more details on sensitivity analysis and the CVODES solver see Serban and Hindmarsh (2005).
Software Features
SUNDIALS was designed in a modular manner with careful attention to flexibility and avoidance of duplication. For example, code modules for the various generic linear solvers (direct and Krylov), which are independent of the top-level solvers, are accessed through interface functions that are specific to the various solvers. In addition, for each top-level solver, these interfaces conform to a standard which permits a user to supply his/her own linear solver module.
Operations on \(N\)-vectors (linear sums, dot products, norms, etc.) are also isolated in an NVECTOR module. The generic NVECTOR module defines abstract vector operations, and links to an actual NVECTOR implementation. The latter can be the serial or parallel (MPI) implementation supplied with SUNDIALS, or one supplied by the user. The rest of SUNDIALS is independent of the NVECTOR implementation.
The user interface to the SUNDIALS solvers consists of various function calls that initialize modules, allocate memory, specify user functions, specify required and optional input parameters, compute and obtain solutions, get optional outputs, and free memory. When there is a default for an optional input, the user can either do nothing or call a Set function to supply a non-default value. Similar Get functions return optional outputs to the user. In the integrators CVODE, CVODES, and IDA, the user has the option of supplying values of \(t\) (output times) at which answers are to be returned, or of having answers returned after each of the time steps taken by the integrator.
For the CVODE, KINSOL, and IDA solvers, SUNDIALS includes modules that interface with Fortran applications. Thus a user's Fortran program can call routines that interface to one of these solvers, and can supply Fortran subroutines that are called by the solver through its interfaces.
Availability, Documentation, and Support
SUNDIALS is distributed open-source, under a BSD-style license. The latest release is always available from the Download page on the SUNDIALS website (https://computation.llnl.gov/casc/sundials/main.html). Older versions are also provided for users of other software packages that depend on them.
Extensive documentation is available in various formats from the Documentation page on the SUNDIALS website. This includes user guides for each individual SUNDIALS solver, documents with detailed descriptions of some of the examples distributed with the suite, preprints of articles related to SUNDIALS, as well as slides from different presentations on SUNDIALS.
Support is provided mainly through the sundials-users mailing list where subscribers can post questions and share information with other members of the SUNDIALS community. Interested users can subscribe to the mailing list from https://computation.llnl.gov/casc/sundials/support/support.html.
References
- Hindmarsh, A. C., et al., SUNDIALS: Suite of Nonlinear and Differential/Algebraic Equation Solvers, ACM Trans. Math. Softw., 31:363-396, 2005. (pdf)
- Serban, R. and Hindmarsh, A.C., CVODES: the Sensitivity-Enabled ODE Solver in SUNDIALS, Proceedings of IDETC/CIE 2005, ASME, Sept. 24-28, 2005, Long Beach, Ca, USA. (pdf)
Internal references
- Bill Gear (2007) Backward differentiation formulas. Scholarpedia, 2(8):3162.
- Eugene M. Izhikevich (2007) Equilibrium. Scholarpedia, 2(10):2014.
- Lawrence F. Shampine and Skip Thompson (2007) Initial value problems. Scholarpedia, 2(3):2861.
- Rob Schreiber (2007) MATLAB. Scholarpedia, 2(7):2929.
- Kendall E. Atkinson (2007) Numerical analysis. Scholarpedia, 2(8):3163.
- Philip Holmes and Eric T. Shea-Brown (2006) Stability. Scholarpedia, 1(10):1838.
- Lawrence F. Shampine and Skip Thompson (2007) Stiff systems. Scholarpedia, 2(3):2855.
External links
See Also
Differential-Algebraic Equations, Initial Value Problem, Numerical Analysis, Ordinary Differential Equations
