NDSolve

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 which might be saved in its private data are accessed only through method functions of the algorithm object.

Required method functions for algorithms used from NDSolve.

These method functions must be defined for the algorithm work with NDSolve. The Step 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 rhs[t + h, y + y] at the step end, so to avoid recomputation, you can add that to the list.

Interpretation of Step method output.

Classical Runge-Kutta

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

This defines a method function to take a single step towards integrating an ODE using the classical fourth order RungeKutta method. Since the method is so simple, it is not necessary to save any private data.

In[15]:=

This defines a method function so that NDSolve can obtain the proper difference order to use for the method. The ___ template is used because the difference order for the method is always 4.

In[16]:=

This defines a method function for the step mode so that NDSolve will know how to control time steps. This algorithm method does not

In[17]:=

This integrates the simple harmonic oscillator equation with fixed step size.

In[18]:=

Out[18]=

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 h and two steps of size h/2. 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.

This integrates the simple harmonic oscillator using the classical fourth-order Runge-Kutta method with steps controlled by using the DoubleStep method.

In[19]:=

Out[19]=

This makes a plot comparing the error in the computed solutions at the step ends. The error for the DoubleStep method is shown in blue.

In[20]:=

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.

In computing the solution above, the option StiffnessTest->False was used for the DoubleStep adaptive step method. For simple equations like the harmonic oscillator that are known not to be stiff, this setting is fine. However, if you plan to work with equations where the properties of the solution are not known, it is better to leave the stiffness test active, since otherwise you may get solutions corrupted from instability. The stiffness test requires the LinearStabilityBoundary property be defined for the method.

One way to gain some understanding of the LinearStabilityBoundary property is to use the OrderStar package, which can plot the linear stability region for methods.

This loads the OrderStar package.

In[22]:=

The linear stability region for a method is the region in the complex z = h plane where the method with step h applied to the test equation x' = x has modulus less than |x|, so that errors are not amplified.

This computes symbolically a step for the CRK4 method as defined above.

In[23]:=

Out[23]=

This computes the factor multiplying x, converting h->z. Note that x had to be added to the result of the method because the "Step" function was defined to only return the increment.

In[24]:=

Out[24]=

The absolute stability region for the method is shown in white.

In[25]:=

You can see from the plot that the linear stability boundary is at about -2.8. You could use this value, but for higher precision solutions, it is good to have an exact value, if possible. Mathematica's symbolic capabilities make this quite feasible.

This computes the factor exactly using the symbolic command Reduce.

In[26]:=

Out[26]=

This sets the "LinearStabilityBoundary" property

In[27]:=

Out[27]=

Now it is possible to compute the solution without the stiffness test disabled. The solution is the same as the one computed above because this equation is not stiff.

In[28]:=

Out[28]=

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.

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

As an example, suppose we want to redefine the Runge-Kutta method shown above 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.

This defines a method function to take a single step towards integrating an ODE using the classical fourth-order RungeKutta method using saved numerical values for the required coefficients.

In[13]:=

This defines a rule that that initialize the algorithm object with the data to be used later.

In[11]:=

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

EventLocator

A method for event location is described within .

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 k and Adams-Moulton corrector of order k + 1, 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.

In[16]:=

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.

This defines a function that computes and saves the values of the constant step size Adams-Moulton coefficients.

In[17]:=

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 NDSolveStateData object. Some other data such as precision is saved for optimization and convenience.

This defines a rule for initializing the AdamsBM method from NDSolve.

In[19]:=

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 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, we set the default of the option to be .

This sets the default for the MaxDifferenceOrder option of the AdamsBM method.

In[20]:=

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

This defines the step mode for the AdamsBM method to always be Automatic. This means that it cannot be called from step controller methods that request fixed step sizes of possibly varying sizes.

In[21]:=

This defines the difference order for the AdamsBM method. This varies with the number of past values saved.

In[22]:=

Finally, 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.

This defines the Step method function for AdamsBM that returns step data according to the templates described above.

In[23]:=

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 j k-2, k-1, and k and determines if the order should be lowered or not.

In[24]:=

This defines a function that determines the best order to use after a successful step.

In[25]:=

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

In[27]:=

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

This solves the harmonic oscillator equation with the Adams method defined earlier.

In[28]:=

Out[28]=

This shows the error of the computed solution. It is apparent that the error is kept within reasonable bounds. Note that after the first few points, the step size has been increased.

In[38]:=

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.

In[4]:=

A lot of time is required for coefficient computation.

In[7]:=

Out[7]=

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, we need to figure out the coefficient update.

In[10]:=

Out[10]=