NDSolve

Switch off messages for variables with similar names and load packages with some test problems and utility functions.

In[4]:=

Euler's method

One of the first and simplest methods for solving initial value problems was proposed by Euler:

Euler's method is not very accurate.

Local accuracy is measured by how high terms are matched with the Taylor expansion of the solution. Euler's method is first-order accurate, so that errors occur one order higher starting at powers of .

Euler's method is implemented in NDSolve as ExplicitEuler.

In[7]:=

Out[7]=

Generalizing Euler's method

The idea of Runge-Kutta methods is to take successive (weighted) Euler steps to approximate a Taylor series. In this way function evaluations (and not derivatives) are used.

For example, consider the one-step formulation of the midpoint method:

The midpoint method can be shown to have a local error of so it is second-order accurate.

The midpoint method is implemented in NDSolve as ExplicitMidpoint.

In[8]:=

Out[8]=

Runge-Kutta methods

Extending the approach in (-1), repeated function evaluation can be used to obtain higher-and-higher-order methods.

Denote the Runge-Kutta method for the approximate solution to an initial value problem at by:

where s is the number of stages.

It is generally assumed that the row-sum conditions hold:

These effectively determine the points in time at which the function is sampled and are a particularly useful device in the derivation of high-order Runge-Kutta methods.

The coefficients of the method are free parameters that are chosen to satisfy a Taylor series expansion through some order in the time step h. In practice other conditions such as stability can also constrain the coefficients.

Explicit Runge-Kutta methods are a special case where the matrix is strictly lower triangular:

It has become customary to denote the method coefficients , , and using a Butcher table, which has the following form for explicit Runge-Kutta methods:

The row-sum conditions can be visualized as summing across the rows of the table.

Notice that a consequence of explicitness is , so that the function is sampled at the beginning of the current integration step.

Example

The Butcher table for the explicit midpoint method (-1) is given by:

FSAL schemes

A particularly interesting special class of explicit Runge-Kutta methods, used in most modern codes, are those for which the coefficient have a special structure known as First Same As Last (FSAL):

For consistent FSAL schemes the Butcher table (-1) has the form:

The advantage of FSAL methods is that the function value at the end of one integration step is the same as the first function value at the next integration step.

The function values at the beginning and end of each integration step are required anyway when constructing the InterpolatingFunction that is used for dense output in NDSolve.

Embedded pairs and local error estimation

An efficient means of obtaining local error estimates for adaptive step size control is to consider two methods of different orders and that share the same coefficient matrix (and hence function values).

These give two solutions:

A commonly used notation is , typically with or .

In most modern codes, including the default choice in NDSolve, the solution is advanced with the more accurate formula so that , which is known as local extrapolation.

The vector of coefficients gives an error estimator avoiding subtractive cancellation of in floating-point arithmetic when forming the difference between (-1) and (-1).

The quantity gives a scalar measure of the error that can be used for step size selection.

Step control

The classical, Integral or I step size controller uses the formula:

where .

The error estimate is therefore used to determine the next step size to use from the current step size.

The notation is explained within and .

Overview

Explicit Runge Kutta pairs of orders 2(1) through 9(8) have been implemented.

Formula pairs have the following properties:

First Same As Last strategy.

Local extrapolation mode, that is the higher-order formula is used to propagate the solution.

Stiffness detection capability (see [HW96], [P83], [R87], [S84]).

Proportional-Integral step size controller for stiff and quasi-stiff systems [G91].

Optimal formula pairs of orders 2(1), 3(2), and 4(3) subject to the above requirements have been derived using Mathematica, and are described in [SS03b].

The 5(4) pair selected is due to Bogacki and Shampine [BS89b, S94] and the 6(5), 7(6), 8(7), and 9(8) pairs are due to Verner.

For the selection of higher-order pairs, issues such as local truncation error ratio and stability region compatibility should be considered (see [S94]) . Various tools have been written to assess these qualitative features.

Methods are interchangeable so that, for example, it is possible to substitute the 5(4) method of Bogacki and Shampine with a method of Dormand and Prince.

Example

Define the Brusselator ODE problem, which models a chemical reaction.

In[9]:=

Out[9]=

This solves the system using an explicit Runge-Kutta method.

In[10]:=

Out[10]=

Extract the interpolating functions from the solution.

In[11]:=

Plot the solution components.

In[12]:=

Method comparison

Sometimes you may be interested to find out what methods are being used in NDSolve.

Here you can see the coefficients of the default 2(1) embedded pair.

In[13]:=

Out[13]=

Another issue is that you may want to compare some of the different methods to see how they perform for a specific problem.

Utilities

We will make use of a utility function CompareMethods for comparing various methods. Some useful NDSolve features of this function for comparing methods are:

The option EvaluationMonitor, which is used to count the number of function evaluations.

The option StepMonitor, which is used to count the number of accepted and rejected integration steps.

This displays the results of the method comparison using a GridBox.

In[14]:=

Reference solution

A number of examples for comparing numerical methods in the literature rely on the fact that a closed-form solution is available, which is obviously quite limiting.

In NDSolve it is possible to get very accurate approximations using arbitrary-precision adaptive step size, adaptive order methods based on Extrapolation.

The following reference solution is computed with a method that switches between a pair of Extrapolation methods, depending on whether the problem appears to be stiff.

In[15]:=

Automatic order selection

When you select DifferenceOrder -> Automatic, the code will automatically attempt to choose the optimal order method for the integration.

Two algorithms have been implemented for this purpose and are described within .

Example

Here is an example that compares built-in methods of various orders, together with the method that is selected automatically.

In[17]:=

Compute the number of integration steps, function evaluations, and the endpoint global error.

In[19]:=

Display the results in a table.

In[20]:=

Out[21]//DisplayForm=

The default method has order seven, which is close to the optimal order of eight in this example; seven function evaluations are needed for the estimates to determine the order.

The selection algorithms are heuristic in that the optimal order may change through the integration but, as the above example illustrates, a reasonable default choice is usually made.

Coefficient plug-in

The implementation of ExplicitRungeKutta provides a default method pair at each order.

Sometimes, however, it is convenient to use a different method, for example:

To replicate the results of someone else.

To use a special-purpose method that works well for a specific problem.

To experiment with a new method

The classical Runge-Kutta method

This shows how to define the coefficients of the classical explicit Runge-Kutta method of order four, approximated to precision p.

In[22]:=

The method has no embedded error estimate and hence there is no specification of the coefficient error vector. This means that the method is invoked with fixed step sizes.

Here is an example of the calling syntax.

In[26]:=

Out[26]=

ode23

This defines the coefficients for a 3(2) FSAL explicit Runge-Kutta pair.

The third-order formula is due to Ralston, and the embedded method was derived by Bogacki and Shampine [BS89a].

In[27]:=

The method is used in the Texas Instruments TI-85 pocket calculator, Matlab and RKSUITE [S94].

Unfortunately it does not allow for the form of stiffness detection that has been chosen.

A method of Fehlberg

This defines the coefficients for a 4(5) explicit Runge-Kutta pair of Fehlberg that was popular in the 1960s [F69].

The fourth order formula is used to propagate the solution and the fifth-order formula is used only for the purpose of error estimation.

In[32]:=

In contrast to the classical Runge-Kutta method of order four, the coefficients include an additional entry that is used for error estimation.

The Fehlberg method is not a FSAL scheme since the coefficient matrix is not of the form (-1); it is a six-stage scheme, but it requires six function evaluations per step because of the function evaluation that is required at the end of the step to construct the InterpolatingFunction.

A DOrmand-PRInce method

Here is how to define a 5(4) pair of Dormand and Prince [DP80]. This is currently the method used by ode45 in Matlab.

In[37]:=

The Dormand-Prince method is a FSAL scheme since the coefficient matrix is of the form (-1); it is a seven-stage scheme, but effectively uses only six function evaluations.

Here is how to call the coefficients of Dormand and Prince in place of the built-in choice. Since the structure of the coefficients includes an error vector, the implementation is able to ascertain that adaptive step sizes can be computed.

In[42]:=

Out[42]=

Method comparison

Here we solve a system using several explicit Runge-Kutta pairs.

For the Fehlberg 4(5) pair, the option EmbeddedDifferenceOrder is used to specify the order of the embedded method.

In[43]:=

The Dormand and Prince 5(4) pair is defined as follows.

In[44]:=

The 5(4) pair of Bogacki and Shampine is the default order-five method.

In[45]:=

Put the methods and some descriptive names together in a list.

In[46]:=

Compute the number of integration steps, function evaluations, and the endpoint global error.

In[48]:=

Display the results in a table.

In[49]:=

Out[50]//DisplayForm=

The default method was the least expensive and provided the most accurate solution.

Method plug-in

This shows how to implement the classical explicit Runge-Kutta method of order four using the method plug-in environment.

This definition is optional since the method in fact has no data. However, any expression can be stored inside the data object. For example, the coefficients could be approximated here to avoid coercion from rational to floating-point numbers at each integration step.

In[51]:=

The actual method implementation is written using a stepping procedure.

In[52]:=

Notice that the implementation closely resembles the description that you might find in a textbook. There are no memory allocation/deallocation statements or type declarations, for example. In fact the implementation works for machine real numbers, machine complex numbers, and even using arbitrary precision software arithmetic.

Here is an example of the calling syntax. For simplicity the method only uses fixed step sizes so we need to specify what step sizes to take.

In[53]:=

Out[53]=

Many of the methods that have been built into NDSolve were first prototyped using top-level code before being implemented in the kernel for efficiency.

Stiffness

Stiffness is a combination of problem, initial data, numerical method, and error tolerances.

Stiffness can arise, for example, in the translation of diffusion terms by divided differences into a large system of ODEs.

In order to understand more about the nature of stiffnes it is useful to study how methods behave when applied to a simple problem.

Linear stability

Consider applying a Runge-Kutta method to a linear scalar equation known as Dahlquist's equation:

The result is a rational function of polynomials where (see for example [L87]).

This utility function finds the linear stability function for Runge-Kutta methods. The form depends on the coefficients and is a polynomial if the Runge-Kutta method is explicit.

Here is the stability function for the fifth-order scheme in the Dormand-Prince 5(4) pair.

In[54]:=

Out[54]=

This function finds the linear stability function for Runge-Kutta methods. The form depends on the coefficients and is a polynomial if the Runge-Kutta method is explicit.

The following package is useful for visualizing linear stability regions for numerical methods for differential equations.

In[55]:=

We can now visualize the absolute stability region .

In[56]:=

Depending on the magnitude of in (-1), if we choose the step size such that , then errors in successive steps will be damped and the method is said to be absolutely stable.

If , then step size selection will be restricted by stability and not by local accuracy.

Stiffness detection

The device for stiffness detection that is used with the option StiffnessTest is described within .

Recast in terms of explicit Runge-Kutta methods the condition for stiffness detection can be formulated as:

with and defined in (-1).

The difference - can be shown to correspond to a number of applications of the power method applied to .

The difference is therefore a good approximation of the eigenvector corresponding to the leading eigenvalue.

The product gives an estimate that can be compared to the stability boundary in order to detect stiffness.

An -stage explicit Runge-Kutta has a form suitable for (-1) if we require that .

The default embedded pairs used in ExplicitRungeKutta all have the form (-1).

An important point is that (-1) is very cheap and convenient; it uses already-available information from the integration and requires no additional function evaluations.

Another advantage of (-1) is that it is straightforward to make use of consistent FSAL methods (-1).

Examples

Select a stiff system modeling a chemical reaction.

In[57]:=

This applies a built-in explicit Runge-Kutta method to the stiff system.

By default stiffness detection is enabled, since it only has a small impact on the running time.

In[58]:=

Out[58]=

The coefficients of the Dormand-Prince 5(4) pair are of the form (-1). However, using these coefficients gives the following message.

In[59]:=

The reason is that the stiffness detection device needs to know where the border of the linear stability boundary is.

Therefore, we need to determine the point at which the stability region crosses the negative real axis.

We can set up an equation in terms of the linear stability function and solve it exactly to find the point where the contour crosses the negative real axis.

In[60]:=

Out[60]=

In general, there may be more than one point of intersection and it may be necessary to choose the appropriate solution.

The following definition sets the value of the linear stability boundary.

In[61]:=

The coefficients can now be used with stiffness detection enabled.

In[62]:=

Out[62]=

The Fehlberg 4(5) method does not have the correct coefficient structure (-1) required for stiffness detection, since .

The default value StiffnessTest -> Automatic checks to see if the method coefficients provide a stiffness detection capability and if they do then stiffness detection is enabled.

Step control revisited

There are some reasons to look at alternatives to the standard Integral step controller (-1) when considering mildly stiff problems.

In[63]:=

This defines an explicit Runge-Kutta method based on the Dormand-Prince coefficients that does not use stiffness detection.

In[64]:=

This solves the system and plots the step sizes that are taken using the utility function StepDataPlot.

In[65]:=

Solving a stiff or mildly stiff problem with the standard step size controller leads to large oscillations, sometimes leading to a number of undesirable step size rejections.

The study of this issue is known as step-control stability.

It can be studied by matching the linear stability regions for the high- and low-sorder methods in an embedded pair.

One approach to addressing the oscillation is to derive special methods, but this compromises the local accuracy.

PI step control

An appealing alternative to Integral step control (-1) is Proportional-Integral or PI step control.

In this case the step size is selected using the local error in two successive integration steps according to the formula:

This has the effect of damping and hence gives a smoother step size sequence.

Note that Integral step control (-1) is a special case of (-1) and is used if a step is rejected:

The option StepSizeControlParameters -> {,} can be used to specify the values of and .

The option StepSizeStartingError can be used to specify the initial scaled error estimate in (-1), which is needed in the first integration step.

Examples

Stiff problem

This defines a method similar to IERK that uses the option StepSizeControlParameters to specify a PI controller.

Here we use generic control parameters suggested by Gustaffson:

In[67]:=

Solving the system again, it can be observed that the step size sequence is now much smoother.

In[68]:=

Nonstiff problem

In general the I step controller (-1) is able to take larger steps for a nonstiff problem than the PI step controller (-1) as the following example illustrates.

In[70]:=

In[71]:=

Using the PI step controller the step sizes are slightly smaller.

In[73]:=

For this reason, the default setting for StepSizeControlParameters is Automatic , which is interpreted as:

Use the I step controller (-1) if StiffnessTest -> False.

Use the PI step controller (-1) if StiffnessTest -> True.

Fine-tuning

Instead of using (-1) directly, it is common practice to use safety factors to ensure that the error is acceptable at the next step with high probability, thereby preventing unwanted step rejections.

The option StepSizeSafetyFactors -> {

} specifies the safety factors to use in the step size estimate so that (-1) becomes:

Here is an absolute factor and typically scales with the order of the method.

The option StepSizeRatioBounds -> {

srmin

srmax

} specifies bounds on the next step size to take such that:

Option summary

Options of the method ExplicitRungeKutta.