The Design of the NDSolve Framework
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)
- Separation of method initialization phase and run-time computation
- Hierarchical and re-entrant numerical methods
- Uniform treatment of rounding errors (see [HLW02], [SS03] and the references therein)
- Vectorized framework based on a generalization of the BLAS model [LAPACK99] using optimized in-place arithmetic.
- Tensor framework that allows families of methods to share one implementation
- Type and precision dynamic for all methods
- Plug-in capabilities that 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 toward 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 re-entrant 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.
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 the problem, the relative and absolute local error tolerances, and the initial step size estimate.
Here is a list of options appropriate for
ExplicitRungeKutta.
| 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 "Norms in NDSolve".
This system dynamically switches type from real to complex during the numerical integration, automatically recompiling as needed.
| Out[2]= |  |
Some basic methods
| | |
| 1 | Explicit Euler | yn+1=yn+hnf (tn, yn) |
| 2 | Explicit Midpoint |  |
| 1 | Backward or Implicit Euler (1-stage RadauIIA) | yn+1=yn+hnf (tn+1, yn+1) |
| 2 | Implicit Midpoint (1-stage Gauss) |  |
| 2 | Trapezoidal (2-stage Lobatto IIIA) |  |
| 1 | Linearly Implicit Euler | (I-hnJ) (yn+1-yn)=hnf (tn, yn) |
| 2 | Linearly Implicit Midpoint |  |
Some of the one-step methods that have been implemented.
Here
yn=yn+1-yn+1/2,
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 Gauss scheme of the
ImplicitRungeKutta method.
