Solving ODEs in MATLAB, 4: Order, Naming Conventions

The digits in the name of a MATLAB ODE solver reflect its order and resulting accuracy. A method is said to have order p if cutting the step size in half reduces the error in one step by a factor of two to the power p+1.

Video Transcript

A Very important property of a numerical method is its order. The accuracy of the method is proportional to a power of the step size. And that power is called the order.

If h is the step size and p is the order, then the error made in one step is proportional to h to the p plus 1. And the error made in traversing an entire interval is proportional to H to the p.

So this means, if you're using a method of order p, and cut the step size in half, you can expect the overall error to be reduced by a factor of 2 to the p. The order of a numerical method is determined by analysis involving Taylor series during the derivation of the method. But we can also do an experiment to determine the order.

That's what this program does. The input is the name of an ODE solver. And then it's going to do a numerical integration of an ordinary differential equation, just involving t. So the result is the value of an integral. The integral of 1 over 1 plus t squared, from 0 to 1.

We know the exact answer is 1/2. So we integrate that differential equation twice, once with a step size of 0.1, and then with a step size of one half that. We integrate the differential equation, take the final value of y for each of those two integrations, compare those values with the exact answer, take the ratio of those two values. That shows how much the error is decreased when we cut the step size in half.

The log base 2 of that ratio is the order. It should be an integer so we can round it to the nearest integer, and return that value as the value in this function. Let's run our experiment first on ODE1.

We step size of 0.1, this method gets the integral as 0.5389, not very accurate. Cut the step size in half, it gets 0.5191. The ratio of those is two. Logarithm base 2 is 1. So ODE1 has order 1.

Now ODE2-- step size 0.1. 0.499. cut it in half, 0.4998. The ratio of those is close to 4. And so ODE2, we find with this experiment, is getting order 2.

Now let's try classical Runge-Kutta. This is why it's so popular. It's very accurate. We step size of 0.1, we get close to 1/2. Cut the step size in half, we get even closer. The ratio of these two is close to 16. To the log base 2 is 4. So ODE4 has, experimentally, order 4.

So we found that, at least with this single experiment, the ODE solvers 1, 2, and 4, have orders 1, 2, and 4. So as you probably realized, this is why we named them ODE 1, 2, and 4.

This brings us to the naming conventions in the functions in the MATLAB ODE suite. All the functions have names that are variations on the theme ODEpq. That means that the method ODEpq uses methods of order p and q So we've been getting a glimpse of that with our names, ODE 1, 2, and 4.

Here's an exercise. Modify order x to do further experiments involving the order of our ODE solvers. Change it to do other integrals. And check out the order of ODE 1, 2, and 4.

Product Focus

Other Resources

1: Euler, ODE1
ODE1 implements Euler's method. It provides an introduction to numerical methods for ODEs and to the MATLAB suite of ODE solvers. Exponential growth and compound interest are used as examples.

3: Classical Runge-Kutta, ODE4
ODE4 implements the classic Runge-Kutta method, the most widely used numerical method for ODEs over the past 100 years. Its major shortcoming is the lack of an error estimate. A simple model of the growth of a flame is an example that is used.

4: Order, Naming Conventions
The digits in the name of a MATLAB ODE solver reflect its order and resulting accuracy. A method is said to have order p if cutting the step size in half reduces the error in one step by a factor of two to the power p+1.

7: Stiffness, ODE23s, ODE15s
A problem is said to be stiff if the solution being sought varies slowly, but there are nearby solutions that vary rapidly, so the numerical method must take small steps to obtain satisfactory results. The flame model demonstrates stiffness.

8: Systems of Equations
An ODE involving higher order derivatives is rewritten as a vector system involving only first order derivatives. The classic Van der Pol nonlinear oscillator is provided as an example. The VdP equation becomes stiff as the parameter is increased.

10: Tumbling Box
Throw a rectangular box with sides of three different lengths into the air. You can get the box to tumble stably about its longest axis or its shortest axis. But if you try to make it tumble about it middle axis, you will find the motion is unstable.

11: Predator-Prey Equations
The classic Lotka-Volterra model of predator-prey competition is a nonlinear system of two equations, where one species grows exponentially and the other decays exponentially in the absence of the other. The program "predprey" studies this model.

12: Lorenz Attractor and Chaos
The Lorenz chaotic attractor was discovered by Edward Lorenz in 1963 when he was investigating a simplified model of atmospheric convection. It is a nonlinear system of three differential equations. The program "lorenzgui" studies this model.