The ODE Test Suite

The ODE test suite is a collection of routines designed to test the functionality of the numerical integration package. It includes a number of integration methods and differential equations. The included methods are the most common and robust generally available. These Butcher Tables and ODEs serve as an example for deriving your own classes in addition to demonstrating the capabilities of the integrator.

The sections below describe first the methods and ODEs, and finally the theory and purpose behind each of the tests.

Methods

There are five methods included in the suite. They range from the simplest Euler method, to more advanced embedded and dense output methods, and including an implicit method (no longer supported). Table 29 lists each of these methods and shows their features; details are provided in the sections that follow.

Table 29 Integration methods used in the ODE Test Suite.

Orders with two values list the order of the method for integration, plus the order of the embedded formula.

Stages with two values list the number of stages without dense output, plus the number of stages with dense output.

Euler

The Euler method is the simplest numerical integration method, as shown in Eq. (4-6). While it is simple and fast, it suffers from stability problems and round-off errors. You should never use the Euler method for anything more than cursory examination of the equations.

Classic Runge-Kutta

This is the method which provides this type of numerical integration with its name. It is a very effective fourth-order method designed to be used in the fixed-step size mode (or occasionally with Richardson Extrapolation).

Fehlberg 4(5)

This is a more advanced fourth-order method with an embedded fifth-order method to be used for step control. It operates in the low-order mode (i.e., uses the fourth-order method to propagate the solution), and so has a very robust error estimate.

DOPRI5

This is one of the most advanced and robust RK methods available at this order. It implements both an embedded formula and dense output, and is designed to operate in the Local Extrapolation mode (i.e., uses the higher, fifth-order method to propagate the solution). While it has more stages than the Fehlberg method, it requires only a few more function evaluations thanks to its FSAL implementation.

This method should handle most of your integration needs. It is also the method employed by MATLAB in its ODE45 routine. You also must use this routine if you wish to perform event handling.

Sixth-order Implicit

This sample method is included as an example of how to create and use an implicit Butcher Table object. Implicit methods are no longer supported by this integration package (due to their inefficiency). It is recommended that use of a Linear Multi-step Method be used for stiff systems. However, LMM are not implemented in this integration package.

Differential Equations

There are five ODEs included in the suite. They range from simple linear equations, to discontinuous derivatives, to stiff systems. Table 30 lists each of these methods and shows their features; details are provided in the sections that follow.

Table 30 Differential equations used in the ODE test suite.

Event types are I for informational, R for Reset, and T for terminal.

Bouncing Ball

This system represents a ball falling due to the effects of gravity. When event finding is turned on, it looks for when the ball strikes the ground plane, at which point the velocity is reversed with a particular coefficient of restitution.

The equations of motion for this system are

(5-1)

where g₀ = 9.8 m/s. The ODE object sets a reset event to find when y₁ crosses zero from above. At this time, t_c, the velocity is reset to be

(5-2)

where k = 0.8 is the chosen coefficient of restitution.

Two additional events are defined for these equations. The second looks for when y₂ crosses zero in either direction, which locates extrema of the motion. This is an informational event only. The final event looks for when t crosses 1.5 from above (i.e., (t-1.5) crosses zero from above). This is a terminal event, and shows that not all events need to incorporate the states.

This function also supplies default values for the absolute and relative tolerances which provide good behavior for the event finder.

Discontinuous Derivative

This simple system implements a single equation which has a discontinuous derivative, namely

(5-3)

where ceil() is the ceiling operator, i.e. returns the next higher integer value.

This equation is useful in demonstrating the robustness of variable step size selection schemes. Since the method is generally very linear, the step sizes will tend to the maximum allowed. However, near the discontinuity the step size must be decreased dramatically to maintain accuracy. This will generally lead to a large number of rejected steps.

Although instructive, a better way of handling this problem would be to treat the discontinuity as a reset event. This would better maintain a consistent step size and better accuracy. This is left as an exercise to the reader.

Linear Chemical Reactor

This linear system represents a simple chemical reactor equation given by

(5-4)

This is a very well behaved linear system that is easy to analyze analytically. For example, this system has the solution

(5-5)

for the initial conditions .

This function supplies a default set of initial conditions and requested output times.

Logistic Curve

This is single-variable logistic curve is another which is relatively simple to analyze, even though it is non-linear. The equation is given by

(5-6)

which has the solution

(5-7)

for the initial condition y₀ = 1.

This function supplies a default set of initial conditions and requested output times.

Stiff ODE

This is a stiff non-linear equation based on the Van der Pol equation

(5-8)

This function supplies a default set of initial conditions and requested output times.

Test Descriptions

This section describes each of the tests performed in the ODE test suite. Seven tests are performed, each with a different objective, although some tests do overlap. Table 31 gives a brief summary of the major items examined by each test. The sections below then provide some theory and detail behind the tests, and the major objective of each.

Table 31 ODE Test Suite test matrix.

Function Test

This test verifies the operation of the integration routines. The linear chemical reactor equation is integrated forward and backward in time, using fixed and variable steps, embedded formula and Richardson Extrapolation. Initial conditions are set using each of the accessor functions.

The results of these tests are shown in a matrix with the pass or fail results.

Accuracy Test

This test checks the accuracy of each of the methods by integrating the linear chemical reactor and logistic equations with many different fixed step sizes. By comparing the ratio of errors between the different step sizes, an estimate of the order of the method can be determined. The theory of this is briefly described here.

Recall from Eq. (4-9) that the local truncation error (the difference between the estimated solution and the true solution given the same initial conditions) is given by

(5-9)

This is the error for a single step, but what about for the method? Each step accumulates this error, so the total error is

(5-10)

where N is the number of steps taken. Thus, the total error is of the proportional to the step size to the power of the order of the method. The error condition given in Eq. (5-10) results in

(5-11)

where M is a constant independent of the step size, h.

Consider now performing the same integration with step sizes that are half the values of the previous integration. The error of this new calculation is

(5-12)

Given Eqs. (5-11) and (5-12), it should be clear that the order of the method, p, can be estimated from

(5-13)

The accuracy test integrates the equations with each of the integrators using step sizes of , for i = 2 to 10. It monitors the size of the error after each step, taking the largest value. Then, the estimate of the order is computed from successive runs according to Eq. (5-13). Note that this is merely an estimate. For some equations, the order is overestimated because the equation has super-linear convergence. Also, when the absolute error becomes very small, round-off errors accumulate, and the estimate again becomes meaningless.

Step Size Test

The step size test measures the accuracy of the step size selection method for a couple of different equations. The Fehlberg 4(5) method is applied to the linear chemical reactor and the discontinuous derivative ODEs. The steps are chosen according to

(5-14)

for various values of r.

At each step, the ratio of the error to the error estimate, R, is computed. The results are tabulated and the percentage of steps that satisfy , for various a, are displayed. Ideally, most of the steps should be less than or equal to 1, with many steps much smaller than that. Richardson Extrapolation should be much more accurate than the embedded formula, although its computational expense generally outweighs this benefit. Also, the chemical reactor, being linear, will be much better behaved than the discontinuous derivative equation. This indicates that there would be a large number of rejected steps dealing with this ODE.

Event Handling Test

This test checks the validity of the event handling capabilities of the integrator and its components. The bouncing ball ODE is integrated using the DOPRI5 method (the only one of the methods presented here that can do event searches). The computed values are determined and compared to the exact analytical values. The two values (for time, position, and velocity) should all match to within the five decimal place of accuracy displayed. The results are written to the MATLAB script eventtst.m for examination.

Output Test

This test checks the output and dense output capabilities of the integrator. The linear chemical reactor equation is integrated using the DOPRI5 method. Variable step size control is chosen, and dense output is used. The computed and exact solutions are written out to the MATLAB script integtst.m for examination.

Error State Test

This function exercises all of the various illegal input and integration errors to verify that the integrator can correctly identify and correct them. It also verifies that integration can be accomplished even if there is an illegal input, and that it is halted if there is an error state. The exception handling mechanism is verified. The results of each of these error checks is displayed.

Miscellaneous Tests

No additional tests are currently performed.