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.

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.

These method functions must be defined for the algorithm to work with NDSolve. The method function should always return a list, but the length of the list depends on whether the step was successful or not. Also, some methods may need to compute the function value at the step end, so to avoid recomputation, you can add that to the list.

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) 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 uses an extrapolation based on integrating a time step with a single step of size and two steps of size . The method 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.

As a system symbol, 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 .

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.

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 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 Mathematica, 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 method functions, it makes sense to define the simple and method functions.

Finally, here is the definition of the 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.

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

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.