Design

Features

Supporting a large number of numerical integration methods for differential equations is a lot of work.

In order to cut down on maintenance and duplication of code, common components are shared between methods.

This approach also allows code optimization to be carried out in just a few central routines.

The principal features of the NDSolve framework are:

Uniform design and interface.

Code reuse (common code base).

Objection orientation (method property specification and communication).

Data hiding.

Separation of method initialization phase and run-time computation.

Numerical methods are hierarchical and reentrant.

Uniform treatment of rounding errors (see [HLW02], [SS02] and the references therein).

Vectorized framework based on a generalization of the BLAS model [LAPACK99] using optimized in-place arithmetic.

Many schemes such as ExplicitRungeKutta and ImplicitRungeKutta use a tensor framework that allows families of methods to share one implementation.

All methods are type and precision dynamic.

Plug-in capabilities allow user extensibility and prototyping.

Specialized data structures.

Common time stepping

A common time stepping mechanism is used for all one-step methods. The routine handles a number of different criteria including:

step sizes in a numerical integration do not become too small in value, which may happen in solving stiff systems

step sizes do not change sign unexpectedly, which may be a consequence of user programming error

step sizes are not increased after a step rejection

step sizes are not decreased drastically towards the end of an integration

specified (or detected) singularities are handled by restarting the integration

divergence of iterations in implicit methods (e.g. using fixed large step sizes)

unrecoverable integration errors (e.g. numerical exceptions)

rounding error feedback (compensated summation) is particularly advantageous for high order methods or methods that conserve specific quantities during the numerical integration.

Data encapsulation

Each method has its own data object that contains information that is needed for the invocation of the method. This includes, but is not limited to, coefficients, workspaces, step size control parameters, step size acceptance/rejection information, and Jacobian matrices. This is a generalization of the ideas used in codes like LSODA ([H83], [P83]).

Method hierarchy

Methods are reentrant and hierarchical, meaning that one method can call another. This is a generalization of the ideas used in the Generic ODE Solving System, Godess (see [O95], [O98] and the references therein), which is implemented in C++.

Initial design

The original method framework design allowed a number of methods to be invoked in the solver.

First revision

This was later extended to allow one method to call another in a sequential fashion, with an arbitrary number of levels of nesting.

The construction of compound integration methods is particularly useful in geometric numerical integration.

Second revision

A more general tree invocation process was required to implement composition methods:

This is an example of a method composed with its adjoint.

Current state

The tree invocation process was extended to allow for a subfield to be solved by each method, instead of the entire vector field.

This example turns up in the ABC flow (see Composition and Splitting).

User extensibility

Built-in methods can be used as building blocks for the efficient construction of special purpose (compound) integrators. User defined methods can also be added.

Method classes

Methods such as ExplicitRungeKutta include a number of schemes of different orders. Moreover, alternative coefficient choices can be specified by the user. This is a generalization of the ideas found in RKSUITE [BGS93].

Automatic selection and user controllability

The framework provides automatic step size selection and method order selection. Methods are user configurable via method options.

For example a user can select the class of ExplicitRungeKutta methods and the code will automatically attempt to ascertain the "optimal" order according to problem, the relative and absolute local error tolerances, and the initial step size estimate.

Here is a list of options appropriate for ExplicitRungeKutta.

In[1]:=

Out[1]=

Shared features

These features are not necessarily restricted to NDSolve since they can also be used for other types of numerical methods.

Function evaluation is performed using a NumericalFunction that dynamically changes type as needed, such as when IEEE floating point overflow or underflow occurs. It also calls Mathematica's compiler Compile for efficiency when appropriate.

Jacobian evaluation uses symbolic differentiation or finite difference approximations, including automatic or user specifiable sparsity detection.

Dense linear algebra is based on LAPACK and sparse linear algebra uses special purpose packages such as UMFPACK.

Common subexpressions in the numerical evaluation of the function representing a differential system are detected and collected to avoid repeated work.

Other supporting functionality that has been implemented is described in Norm and ScaledVectorNorm.

This system dynamically switches type from real to complex during the numerical integration, automatically recompiling as needed.

In[1]:=

Out[1]=

Some basic methods

The following table gives some of the one-step methods that have been implemented:

Here , I denotes the identity matrix and J denotes the Jacobian matrix

Although the implicit midpoint method has not been implemented as a separate method, it is available through the one stage Guass scheme of the ImplicitRungeKutta method.

Acknowledgements

A number of people have helped in shaping the design and implementation of the methods and framework in NDSolve. We would like to thank them for sharing their experience and their enthusiasm for the project.

Thanks to Giulia Spaletta who collaborated in the development of a number of methods.

Thanks to Alan Hindmarsh for allowing us to use the codes LSODA, CVODE and IDA.

Thanks also to the following people for fruitful discussions: John Butcher, Jeff Cash, Ernst Hairer, Robert McLachlan, Per Christian Moan, Reinout Quispel.

Aspects of the new framework have been presented at various symposia and workshops over the years. We would particularly like to thank the groups at Bergen and Trondheim, Norway, and Bari, Italy.

Thanks also to Andrew Hunt of Wolfram Research for a number of stylistic improvements to the documentation.

Robert Knapp and Mark Sofroniou.