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 (RefLink[NDSolve, paclet:ref/NDSolve]) framework are:

Common Time Stepping

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

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.

(RefLink[NDSolve, paclet:ref/NDSolve]) ⟶ "ExplicitRungeKutta"

(RefLink[NDSolve, paclet:ref/NDSolve]) ⟶ "ImplicitRungeKutta"

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.

(RefLink[NDSolve, paclet:ref/NDSolve]) ⟶ "Extrapolation" ⟶ "ExplicitMidpoint"

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

(RefLink[NDSolve, paclet:ref/NDSolve]) ⟶ "Projection" ⟶ "ExplicitRungeKutta"

Second Revision

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

(RefLink[NDSolve, paclet:ref/NDSolve]) ({{⟶, "Composition" ({{↗, "ExplicitEuler"}, {⋮, ⋮}, {⟶, ({{"ImplicitEuler", }})}, {⋮, ⋮}, {↘, ({{"ExplicitEuler", }})}})}})

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 subsection of "Composition and Splitting Methods for NDSolve".

(RefLink[NDSolve, paclet:ref/NDSolve]) ({{⟶, "Splitting"  f=(Subscript[f, 1])+(Subscript[f, 2])  ({{↗, "LocallyExact"  (Subscript[f, 1])}, {⟶, ({{"ImplicitMidpoint"  (Subscript[f, 2]), }})}, {↘, ({{"LocallyExact"  (Subscript[f, 1]), }})}})}})

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".

Options[NDSolve`ExplicitRungeKutta]

MethodMonitor

In order to illustrate the low-level behavior of some methods, such as stiffness switching or order variation that occurs at run time, a new "MethodMonitor" has been added.

This fits between the relatively coarse resolution of "StepMonitor" and the fine resolution of "EvaluationMonitor" .

This feature is not officially documented and the functionality may change in future versions.

Shared Features

These features are not necessarily restricted to (RefLink[NDSolve, paclet:ref/NDSolve]) since they can also be used for other types of numerical methods.

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

y[1/2]/.NDSolve[{y'[t]==Sqrt[y[t]]-1,y[0]==1/10},y,{t,0,1},Method->"ExplicitRungeKutta"]

Some Basic Methods

Some of the one-step methods that have been implemented.

Here Δ (Subscript[y, n])==(Subscript[y, n+1])-(Subscript[y, n+1/2]), I denotes the identity matrix, and J denotes the Jacobian matrix ((∂f)/(∂y))((Subscript[t, n]),(Subscript[y, n])).

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.