Solving ODEs in MATLAB, 5: Estimating Error, ODE23

ODE23 compares methods of order two and three to automatically choose the step size and maintain a specified accuracy. It is the simplest MATLAB solver that has modern features such as automatic error estimate and continuous interpolant. ODE23 is suitable for coarse accuracy requirements such as computer graphics.

Video Transcript

Software that implements modern numerical methods has two features that aren't present in codes like ODE4 and classical Runge-Kutta. The methods in the software can estimate error and provide automatic step size control. You don't specify the step size h. You specify an accuracy you want. And the methods estimate the errors as they go along and adjust the step size accordingly. And they provide a fully accurate continuous interpolant. They don't just provide the solution at the discrete set of points. They provide a function that defines the solution everywhere in the interval. And so you can plot it, find zeroes of the function, provide a facility called event handling, and so on.

Larry Shampine is an authority on the numerical solution of ordinary differential equations. He is the principal author of this textbook about solving ODEs with MATLAB. He's a, now, emeritus professor at the Southern Methodist University in Dallas. And he's been a long time consultant to the MathWorks about the development of our ODE Suite. Shampine and his student, Przemyslaw Bogacki, published this method in 1989. And it's the basis for ODE23, the first of the methods we will use out of the MATLAB ODE Suite.

The basic method is order three. And the error estimate is based on the difference between the order three method and then the underlying order two method. There are four slopes involved.

The first one is the value of the function at the start of the interval. But that's based on something called FSAL, first same as last, where that slope is most likely left over from the previous step. If the previous step was successful, this function value is the same as the last function value from the previous step.

That slope is used to step into the middle of the interval, function is evaluated there. That slope is used to step 3/4 of the way across the interval and a third slope obtained there. Then these three values are used to take the step. yn plus 1 is a linear combination of these three function values. Then the function is evaluated to get a fourth slope at the end of the interval. And then, these four slopes are used to estimate the error.

The error estimate here is the difference between yn plus 1 and another estimate of the solution that's obtained from a second order method that we don't actually evaluate. We just need the difference between that method and yn plus 1 to estimate the error.

This estimated error is compared with a user-supplied tolerance. If the estimated error is less than a tolerance, then the step is successful. And this fourth slope, s4, becomes the s1 of the next step.

If the answer is bigger than the tolerance, then the error could be the basis for adjusting the step size. In either case, the error estimate is the basis for adjusting the step size for the next step. This is the Bogacki-Shampine Order 3(2) Method that's the basis for ODE 23.

Let's look at some very simple uses of ODE23 just to get started. I'm going to take the differential equation y prime is equal to y. So I'm going to compute e to the t. And just call ODE23 on the interval from 0 to 1, with initial value 1. No output arguments. If I call it ODE23, it just plots the solution. Here it is. It just produces a plot. It picks a step size, goes from 0 to 1, and here it gets the final value of e-- 2.7 something.

If I do supply output arguments. I say t comma y equals ODE23, it comes back with values of t and y. ODE23 picks the values of t it wants. This is a trivial problem. It ends up picking a step size of 0.1. After it gets started, it chooses an initial step size of .08 for whatever error tolerances. And the final value of y is 2.718, which is the value of e.

So these are the two simple uses of ODE23. If you don't supply any output arguments, it draws a graph. If you do supply output arguments, t and y, it comes back with the values of t and y choosing the values of t to meet the error. The default error tolerances is 10 to the minus 3. So this value is going to be accurate to three digits. And sure enough that's what we got.

Now let's try something a little more challenging to see the automatic error-controlled step size choice in action. Set a equal to a quarter. And then set y0 equal to 15.9. If I would set it to 16, which is 1 over a squared, I'd run into a singularity. Now the differential equation is y prime is equal to 2 (a minus t) times y squared. I'm going to integrate this with the ODE23 on the interval from 0 to 1 starting at y0, and saving the results in t and y, and then plotting them. So here's my plot command, and there is the solution.

So there is a near singularity at a. It nearly blows up. And then it settles back down. So the points are bunched together as you go up to the singularity and come back down, but then get farther apart as the solution settles down. And the ODE solver is able to take bigger steps.

To see what steps were actually taken, let's compute the difference of t, and then plot that. So here are the step sizes that were taken. And we see that a small step size was taken near the almost singularity at that 0.25. And then as we get towards the end of the interval, a larger step size is taken. And then, finally, the step size just to reach the end of the interval is taken as the last step. So that's the automatic step size choice of ODE23.

BS23 has a nice natural interpolant that goes along with it that's actually been known for over 100 years. It's called Hermite Cubic Interpolation. We know that two points determine a straight line. Well, two points and two slopes determine a cubic.

On each interval we have the values of y and yn plus 1. We also have two slopes, namely this. We have the derivatives at the end points, yn prime and yn plus 1 prime, that's the values of the differential equation at those points. So those four values determine a cubic that goes through those two points and has those two slopes.

This cubic allows the software to evaluate the solution at any point in the interval without additional cost as defined by addition evaluations of the function f. This can be used to draw graphs of the solution, nice smooth graphs of the solution, find zeroes of the solution, do event handling, and so on. Another feature provided by ODE23.

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.

This website uses cookies to improve your user experience, personalize content and ads, and analyze website traffic. By continuing to use this website, you consent to our use of cookies. Please see our Privacy Policy to learn more about cookies and how to change your settings.