Mathematica Tutorial

ExplicitRungeKutta Method for NDSolve

Introduction

This loads packages containing some test problems and utility functions.

In[3]:=

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 h².

Euler's method is implemented in NDSolve as ExplicitEuler.

In[5]:=

Out[5]=

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 O (h³) so it is second-order accurate.

The midpoint method is implemented in NDSolve as ExplicitMidpoint.

In[6]:=

Out[6]=

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 t_n+1=t_n+h, 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 A is strictly lower triangular:

It has become customary to denote the method coefficients c=[c_i]^T, b=[b_i]^T, and A=[a_{i, j}] 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 c₁=0, so that the function is sampled at the beginning of the current integration step.

Example

The Butcher table for the explicit midpoint method (2) 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 coefficients have a special structure known as First Same As Last (FSAL):

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

The advantage of FSAL methods is that the function value k_s at the end of one integration step is the same as the first function value k₁ 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 p and

that share the same coefficient matrix (and hence function values).

These give two solutions:

A commonly used notation is

, typically with

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 y_n in floating-point arithmetic when forming the difference between (4) and (5).

The quantity

err_n

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 Tol/

err_n

is explained within "Norms in NDSolve".

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 already stated requirements have been derived using Mathematica, and are described in [SS04].

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[7]:=

Out[7]=

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

In[8]:=

Out[8]=

Extract the interpolating functions from the solution.

In[9]:=

Plot the solution components.

In[10]:=

Out[10]=

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[11]:=

Out[11]=

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

Utilities

You 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[12]:=

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[13]:=

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 "SymplecticPartitionedRungeKutta Method for NDSolve".

Example 1

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

This selects the order of the methods to choose between and makes a list of method options to pass to NDSolve.

In[15]:=

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

In[17]:=

Display the results in a table.

In[18]:=

Out[19]//DisplayForm=

The default method has order nine, which is close to the optimal order of eight in this example. One function evaluation is needed during the initialization phase to determine the order.

Example 2

A limitation of the previous example is that it did not take into account the accuracy of the solution obtained by each method, so that it did not give a fair reflection of the cost.

Rather than taking a single tolerance to compare methods, it is preferable to use a range of tolerances.

The following example compares various ExplicitRungeKutta methods of different orders using a variety of tolerances.

This selects the order of the methods to choose between and makes a list of method options to pass to NDSolve.

In[20]:=

The data comparing accuracy and work is computed using CompareMethods for a range of tolerances.

In[22]:=

The work-error comparison data for the various methods is displayed in the following logarithmic plot, where the global error is displayed on the vertical axis and the number of function evaluations on the horizontal axis. The default order selected (displayed in red) can be seen to increase as the tolerances are increased.

In[23]:=

Out[23]=

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

Ideally, a selection of different problems should be used for benchmarking.

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[24]:=

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[25]:=

Out[25]=

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

This defines function for computing the coefficients to a desired precision.

In[26]:=

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.

This defines the function for computing the coefficients to a desired precision.

In[31]:=

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 (6); 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.

This defines function for computing the coefficients to a desired precision.

In[36]:=

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

Here is how the coefficients of Dormand and Princeton be used 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[41]:=

Out[41]=

Method comparison

Here you 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[42]:=

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

In[43]:=

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

In[44]:=

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

In[45]:=

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

In[47]:=

Display the results in a table.

In[48]:=

Out[49]//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[50]:=

The actual method implementation is written using a stepping procedure.

In[51]:=

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 you need to specify what step sizes to take.

In[52]:=

Out[52]=

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 stiffness 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 R (z) where z=h

(see for example [L87]).

This utility function finds the linear stability function R (z) 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[53]:=

Out[53]=

This function finds the linear stability function R (z) 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[54]:=

You can now visualize the absolute stability region

R (z)

=1.

In[55]:=

Out[55]=

Depending on the magnitude of

in (8), if you choose the step size h such that

R (h

)

<1, then errors in successive steps will be damped and the method is said to be absolutely stable.

R (h

)

>1, 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 "StiffnessTest Method Option for NDSolve".

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

with g_i and k_i defined in (9).

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

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 s-stage explicit Runge-Kutta has a form suitable for (10) if you require that c_s-1=c_s=1.

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

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

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

Examples

Select a stiff system modeling a chemical reaction.

In[56]:=

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[57]:=

Out[57]=

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

In[58]:=

The reason is that the stiffness detection device needs to know where the border of the linear stability boundary is. Therefore, you need to determine the point at which the stability region crosses the negative real axis.

You 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[59]:=

Out[59]=

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[60]:=

The coefficients can now be used with stiffness detection enabled for the default setting StiffnessTest->Automatic.

In[61]:=

Out[61]=

The Fehlberg 4(5) method does not have the correct coefficient structure (16) required for stiffness detection, since c_s=1/2≠1.

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 (17) when considering mildly stiff problems.

This system models a chemical reaction.

In[62]:=

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

In[63]:=

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

In[64]:=

Out[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-order 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 (18) 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 (19) is a special case of (20) and is used if a step is rejected:

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

The scaled error estimate in (21) is taken to be

err_n-1

err_n

for 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 you use generic control parameters suggested by Gustafsson:

This specifies the step control parameters.

In[66]:=

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

In[67]:=

Out[68]=

Nonstiff problem

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

Select and solve a non stiff system using the I step controller.

In[69]:=

In[70]:=

Out[71]=

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

In[72]:=

Out[73]=

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

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

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

Fine-tuning

Instead of using (26) 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->{s₁, s₂} specifies the safety factors to use in the step size estimate so that (27) becomes:

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

The option StepSizeRatioBounds->{sr_min, sr_max} specifies bounds on the next step size to take such that:

Option summary

option name	default value
Coefficients	EmbeddedExplicitRungeKuttaCoefficients	specify the coefficients of the explicit Runge-Kutta method
DifferenceOrder	Automatic	specify the order of local accuracy
EmbeddedDifferenceOrder	Automatic	specify the order of the embedded method in a pair of explicit Runge-Kutta methods
StepSizeControlParameters	Automatic	specify the PI step control parameters (see (28))
StepSizeRatioBounds	{,4}	specify the bounds on a relative change in the new step size (see (29))
StepSizeSafetyFactors	Automatic	specify the safety factors to use in the step size estimate (see (30))
StiffnessTest	Automatic	specify whether to use the stiffness detection capability

Options of the method ExplicitRungeKutta.

The default setting of Automatic for the option DifferenceOrder selects the default coefficient order based on the problem, initial values and local error tolerances, balanced against the work of the method for each coefficient set.

The default setting of Automatic for the option EmbeddedDifferenceOrder specifies that the default order of the embedded method is one lower than the method order. This depends on the value of the DifferenceOrder option.

The default setting of Automatic for the option StepSizeControlParameters uses the values {1, 0} if stiffness detection is active and {3/10, 2/5} otherwise.

The default setting of Automatic for the option StepSizeSafetyFactors uses the values {17/20, 9/10} if the I step controller (31) is used and {9/10, 9/10} if the PI step controller (32) is used. The step controller used depends on the values of the options StepSizeControlParameters and StiffnessTest.

The default setting of Automatic for the option StiffnessTest will activate the stiffness test if a LinearStabilityBoundary method property is specified and if the coefficients have the form (33).