Introduction

The control mechanisms set up for NDSolve enable you to define your own numerical integration algorithms and use them as specifications for the Method option of NDSolve.

NDSolve accesses its numerical algorithms and the information it needs from them in an object-oriented manner. At each step of a numerical integration, NDSolve keeps the method in a form so that it can keep private data as needed.

The structure for method data used in NDSolve.

NDSolve does not access the data associated with an algorithm directly, so you can keep the information needed in any form that is convenient or efficient to use. The algorithm and information that might be saved in its private data are accessed only through method functions of the algorithm object. The key interaction between a method algorithm object and NDSolve is through a "Step" function.

Required "Step" function used by NDSolve.

The following is a list of important properties of the method and its "Step" function that should be defined.

Properties that should be defined for using a method algorithm object obj in NDSolve.

The table below gives specifications that can be used for step inputs and/or outputs.

short name	full name	component
"F"	"Function"	the function given the right-hand side of the ODE system or the residual for a DAE system
"T"	"Time"	current time
"H"	"TimeStep"	time step
"X"	"DependentVariables"	current values of dependent variables other than discrete
"XI"	"DependentVariablesIncrement"	such that gives the values of the dependent variables after the step; if is given as $Failed, then the step is rejected and will be recomputed (the time step returned should be smaller)
"XP"	"TemporalDerivatives"	the time derivatives of the dependent variables
"D"	"DiscreteVariables"	the discrete variables that take on a continuous range of values
"ID"	"IndexedDiscreteVariables"	the discrete variables that take on a finite set of values
"P"	"Parameters"	the parameters other than ones for which sensitivity is being computed
"SP"	"SensitivityParameters"	the parameters for which sensitivity is being computed
"SD"	"SolutionData"	the solution data list
"Obj"	"MethodObject"	a new setting for obj that will be used for the next step
"State"	"StateData"	the NDSolve`StateData object

Step input and output specifications.

The solution data is the list of values with components relating to different variable types. The solution data in the active direction can be obtained from an NDSolve`StateData object state using state@"SolutionData"["Active"]. Many of the specifications listed above correspond directly to solution data components.

For input, the specification for the function should be given with the arguments it will be given. These arguments should be solution data components. All can be used to indicate that all arguments used by NDSolve should be expected.

For output, if the method only returns some of the values in some cases (e.g. the method object does not change every step) then you can use Null for any output except the increment to indicate that it should not be used.

Classical Runge–Kutta

Here is an example of how to set up and access a simple numerical algorithm.

Generally using a fixed step size is less efficient than allowing the step size to vary with the local difficulty of the integration. Modern explicit Runge–Kutta methods (accessed in NDSolve with Method->"ExplicitRungeKutta") have a so-called embedded error estimator that makes it possible to very efficiently determine appropriate step sizes. An alternative is to use built-in step controller methods that use extrapolation. The method "DoubleStep" uses an extrapolation based on integrating a time step with a single step of size and two steps of size . The method "Extrapolation" does a more sophisticated extrapolation and modifies the degree of extrapolation automatically as the integration is performed, but is generally used with base methods of difference orders 1 and 2.

The fixed step size ended up with smaller overall error mostly because the steps are so much smaller; it required more than three times as many steps. For a problem where the local solution structure changes more significantly, the difference can be even greater.

A facility for stiffness detection is described within "DoubleStep Method for NDSolve".

For more sophisticated methods, it may be necessary or more efficient to set up some data for the method to use. When NDSolve uses a particular numerical algorithm for the first time, it calls an initialization function. You can define rules for the initialization that will set up appropriate data for your method.

Initializing a method from NDSolve.

NDSolve`InitializeMethod is protected, so to attach rules to it, you would need to unprotect it first. It is better to keep the rules associated with your method. A tidy way to do this is to make the initialization definition using TagSet as shown earlier.

As an example, suppose you want to redefine the Runge–Kutta method shown earlier so that instead of using the exact coefficients 2, 1/2, and 1/6, numerical values with the appropriate precision are used instead to make the computation slightly faster.

Saving the numerical values of the numbers gives between 5 and 10 percent speedup for a longer integration using "DoubleStep".

Crank–Nicolson

This is an example of how to set up an implicit method. To solve the system of ODEs , the scheme for a time step of size is , where and . In general, for nonlinear , the equations need to be solved with Newton iteration.

Instead of using a true Newton iteration where the Jacobian is updated at each step, the Jacobian will be computed once and its LU decomposition will be used repeatedly to solve the linear equations. As long as time steps are not too large, this will generally converge as well as the full Newton method and each iteration beyond the first is substantially cheaper.

To avoid possible problems with the Newton iteration not converging, the method will be set up to take method options that control the maximum number of iterations and the convergence tolerance.

For a linear problem, the solution is given by the first iteration. However, for a nonlinear problem, convergence may take longer.

Since the Crank–Nicolson method is usually considered regarding solving PDEs, here is an example of the method solving the wave equation. Since this is a linear equation, convergence occurs in 1 iteration so the method is quite fast.

Invoking Another Method

One of the strengths of the NDSolve framework comes from the capability to nest methods in various ways.

Using submethods.

NDSolve`InitializeSubmethod effectively calls NDSolve`InitializeMethod after separating out the method name from any suboptions that are given. If successful, it returns the method object that you can later use in NDSolve`InvokeMethod.

NDSolve`InvokeMethod effectively constructs a "Step" function for the submethod. Since the arguments for the function are specified by the "StepInput" specification for the submethod, the function f given to NDSolve`InvokeMethod should have all the arguments used by NDSolve and can be specified in the caller "StepInput" specification by "Function"[All]. Note that for built-in methods, there may not be an explicit "Step" function and NDSolve`InvokeMethod handles these using efficient internal code. The result returned by NDSolve`InvokeMethod is always constructed from the submethod's actual output to match a "StepOutput" of {"TimeStep","SolutionData","MethodObject"}.

Here is an example that simply calls a submethod to provide simple monitoring of the steps and difference order.

For a more complicated integration the method can be used with Reap to collect data.

Adams Methods

In terms of the NDSolve framework, it is not really any more difficult to write an algorithm that controls steps automatically. However, the requirements for estimating error and determining an appropriate step size usually make this much more difficult from both the mathematical and programming standpoints. The following example is a partial adaptation of the Fortran DEABM code of Shampine and Watts to fit into the NDSolve framework. The algorithm adaptively chooses both step size and order based on criteria described in [SG75].

The first stage is to define the coefficients. The integration method uses variable step-size coefficients. Given a sequence of step sizes , where is the current step to take, the coefficients for the method with Adams–Bashforth predictor of order and Adams–Moulton corrector of order , such that

where the are the divided differences.

This defines a function that computes the coefficients and , along with , that are used in error estimation. The formulas are from [HNW93] and use essentially the same notation:

hlist is the list of step sizes from past steps. The constant-coefficient Adams coefficients can be computed once, and much more easily. Since the constant step size Adams–Moulton coefficients are used in error prediction for changing the method order, it makes sense to define them once with rules that save the values.

The next stage is to set up a data structure that will keep the necessary information between steps and define how that data should be initialized. The key information that needs to be saved is the list of past step sizes, hlist, and the divided differences, . Since the method does the error estimation, it needs to get the correct norm to use from the NDSolve`StateData object. Some other data such as precision is saved for optimization and convenience.

hlist is initialized to {} since at initialization time there have been no steps. is initialized to the derivative of the solution vector at the initial condition since the 0 ^th divided difference is just the function value. Note that is a matrix. The third element, which is initialized to a zero vector, is used for determining the best order for the next step. It is effectively an additional divided difference. The use of the other parts of the data is clarified in the definition of the stepping function.

The initialization is also set up to get the value of an option that can be used to limit the maximum order of the method to use. In the code DEABM, this is limited to 12, because this is a practical limit for machine-precision calculations. However, in the Wolfram System, computations can be done in higher precision where higher-order methods may be of significant advantage, so there is no good reason for an absolute limit of this sort. Thus, you set the default of the option to be .

Before proceeding to the more complicated "Step" method functions, it makes sense to define the simple "StepMode" and "DifferenceOrder" method functions.

Finally, here is the definition of the "Step" function. The actual process of taking a step is only a few lines. The rest of the code handles the automatic step size and order selection following very closely the DEABM code of Shampine and Watts.

There are a few deviations from DEABM in the code here. The most significant is that coefficients are recomputed at each step, whereas DEABM computes only those that need updating. This modification was made to keep the code simpler, but does incur a clear performance loss, particularly for small to moderately sized systems. A second significant modification is that much of the code for limiting rejected steps is left to NDSolve, so there are no checks in this code to see if the step size is too small or the tolerances are too large. The stiffness detection heuristic has also been left out. The order and step-size determination code has been modularized into separate functions.

This defines a function that constructs error estimates for , , and and determines if the order should be lowered or not:

This defines a function that determines the best step size to use after a successful step of size :

Once these definitions are entered, you can access the method in NDSolve by simply using Method->AdamsBM.

Where this method has the potential to outperform some of the built-in methods is with high-precision computations with strict tolerances. This is because the built-in methods are adapted from codes with the restriction to order 12.

A lot of time is required for coefficient computation.

This is not using as high an order as might be expected.

In any case, about half the time is spent generating coefficients, so to make it better, you need to figure out the coefficient update.