MATH2071: LAB #5: Implicit ODE Methods

TABLE OF CONTENTS

• Introduction
• Stiff Systems
• The Backward Euler Method
• The Trapezoid Method
• The Starting Point Problem
• BDF Methods
• Direction Field Plots
• Summary
• ASSIGNMENT

Introduction

The explicit methods that we discussed last time are well suited to handling a large set of ODE's. The Adams-Bashforth methods are very inexpensive, using old data to predict the next value. To get a higher order method, we simply use more data. The Runge-Kutta methods have a complicated form, and are expensive, because they evaluate the derivative function at lots of extra points, but these methods can be a bit more accurate, since they are sampling the derivative function right near the place where the next point is to be computed. And high-order versions of Runge-Kutta methods are available to give very accurate results.

This two sets of methods perform poorly, however, for a class of "stiff" problems, that can come up in chemistry and biology. We will examine two families of implicit methods suitable for such problems. We will find that the implementation of an implicit method has at least two extra complications we didn't see with the explicit method: getting a starting value, and solving (or not solving!) the nonlinear equation.

Stiff Systems

I've warned you that there are problems which defeat the explicit Runge-Kutta and Adams-Bashforth methods. In fact, for such problems, the higher order methods perform even more poorly than the low order methods. These problems are called "stiff" ODE's. We will only look at some very simple examples.

Consider the differential equation

y' = - 5 * y
y(0) = 1

FUNCTION M FILE: Make a copy of model_ode.m, called stiff_ode.m, which has the form

function yprime = stiff_ode ( x, y )
yprime = - 5 * y;

COMPUTATION: Compute the numerical solution of the stiff ODE, from x = 0.0 to x = 2.0. Record the value of the solution at x = 2.0. Try the Euler method with several stepsizes. Try several methods at the same stepsize. Plot your results and compare them visually to the exact solution.

        Stepsize   Euler           AB2             RK3    

        0.2        __________      __________      __________ 
        0.1        __________      
        0.05       __________      
        0.025      __________      
      

The Backward Euler Method

The Backward Euler method is a variation of Euler's method. Before we say anything more about it, let's take a long hard look at the algorithm:

x_k+1 = x_k + dx
y_k+1 = y_k + dx * f(x_k+1,y_k+1)

If you sneeze at the wrong time, you might think there is no difference between this method and Euler's method. But look carefully. This is not a "recipe", the way most formulas are. The equation defining y_k+1 is implicit, and has to be handled in a more careful way than we are used to. So our instinctive reaction might be, well who needs this kind of complication? It happens that there are ODE's for which this approach is the best known way to produce a solution, where Runge-Kutta and Adams-Bashforth methods would produce only gibberish.

Another thing to notice about this method is that it's exactly Euler's method "time-reversed". In other words, if we started at (x₀,y₀), with a given differential equation, then we could use Euler's method to go "right" to (x_n,y_n), or we could start at the left endpoint and use the backward Euler method to go back to the starting point.

But to use the backward Euler method in a real situation, we have to deal with the implicit nature of the formula. In fact, what most algorithms do is "fudge" it. They use a reasonably good guess for what the answer is, make one small adjustment that they hope is an improvement, and declare victory. This process is so common that it has a name, and these methods are often called predictor-corrector methods.

So now let's look at a very simple implementation of the backward Euler method. We are going to treat the formula as a machine. We dump an estimated value of y_k+1 into the right hand side of the formula, and out comes a "corrected" estimate. Where do we get our initial estimate? We use the forward Euler method, as a "predictor". Now generally this is not going to produce a quantity that actually satisfies the backward Euler formula. If we put the corrected value in the right hand side, we'd get yet another value out, and so on. The hope is, that at least for small enough steps, this is a convergent process and that one corrector step might be "good enough".

This may sound like a pretty rickety idea, but it's a start. This technique of handling an implicit ODE approximation formula by a prediction step and a single correction step is called PECE. (The E's represent evaluations of the derivative function).

Here's what our PECE version of the backward Euler method looks like:

x_k+1 = x_k + dx
y_p = y_k + dx * f(x_k,y_k) (predict)
y_k+1 = y_k + dx * f(x_k+1,y_p) (correct)

SOLVER M FILE: make a copy of your euler.m file, naming it be_pece.m. Modify it to have the form:

function [ x, y ] = be_pece ( f, x_range, y_initial, nstep )

COMPUTATION: Compute the numerical solution of the model ODE, from x = 0.0 to x = 2.0, using the backward Euler PECE method, with stepsizes of 0.2, 0.1, 0.05 and 0.025. For each case, record the value of the numerical solution at x = 2.0. (The exact value is exp(2.0) = 7.389056, and the Euler method values are included for comparison.)

        Stepsize   Euler           BE_PECE    

        0.2        6.191736      __________ 
        0.1        6.727500      __________ 
        0.05       7.039989      __________ 
        0.025      7.209566      __________ 
      

COMPUTATION: Repeat the previous computation, but now do it for stiff_ode, whose exact solution at 2 is exp(-10.0) or about 0.000045:

        Stepsize   BE_PECE    

        0.2        __________ 
        0.1        __________ 
        0.05       __________ 
        0.025      __________ 
      

COMPUTATION: From the web page, copy the files: stiff_ode_partial.m and be_newton.m. In order to solve the problem with this method, you need to use a slightly more complicated command, such as:

[x,y] = be_newton ( 'stiff_ode', 'stiff_ode_partial', [0,2], 1, 10 );

        Stepsize   BE_NEWTON 

        0.2        __________ 
        0.1        __________ 
        0.05       __________ 
        0.025      __________ 
      

The Trapezoid Method

The trapezoid method can be derived from the trapezoid rule for integration. It has a simple form:

x_k+1 = x_k + dx
y_k+1 = y_k + dx * ( f(x_k,y_k) + f(x_k+1,y_k+1) ) / 2

To properly use the trapezoidal formula, we must be prepared to solve a nonlinear equation. But it is very common to use a watered-down version of the formula, which is not guaranteed to have the good stability properties, but which is easy to solve. Instead of using Newton's method, we can use the method of successive substitution or functional iteration. In such a method, we have a nonlinear equation which can be written as

y = f(y)

With these modificiations, our version of the trapezoidal method can be termed a "PECECE" method, and can be written:

x_k+1 = x_k + dx
y_k+1 = y_k + dx * f(x_k,y_k) (predict)
y_k+1 = y_k + dx * ( f(x_k,y_k) + f(x_k+1,y_k+1) ) / 2 (correction #1)
y_k+1 = y_k + dx * ( f(x_k,y_k) + f(x_k+1,y_k+1) ) / 2 (correction #1)

SOLVER M FILE: make a copy of your euler.m file, naming it trapezoid.m. Modify it to have the form:

function [ x, y ] = trapezoid ( f, x_range, y_initial, nstep )

FUNCTION M FILE: make a copy of your model_ode.m file, naming it bump_ode.m. It should evaluate the right hand side of the ODE

y' = - y²

COMPUTATION: Test your method by computing the numerical solution of bump_ode.m, using initial condition y(0)=1, to x = 5.0, using the trapezoid PECECE method, with stepsize 0.25. Record the value of the numerical solution at regularly spaced intervals. The book carried out this experiment, and we list the values it gives.

        X       Exact Y        Book Y     TRAPEZOID_PECECE

        0.0       1.000        1.000         1.0
        1.0       0.500        0.496         __________ 
        2.0       0.333        0.331         __________ 
        3.0       0.250        0.249         __________ 
        4.0       0.200        0.199         __________ 
        5.0       0.167        0.166         __________ 
      

The Starting Point Problem

If we actually try to solve the nonlinear equation defining an implicit formula, then it doesn't matter too much what initial value we use, as long as we're close enough for the nonlinear equation solver to converge. So I could use an Euler prediction with be_newton, or even just the old value of the solution, y_k.

However, in a PECE method (or PECECE, or whatever) method, where we only take a few corrector steps, a poor starting value may mean the corrector can't get us a good result. So it becomes more important to try to get a reasonably good predictor method. What is common practice, when using an Adams-Moulton corrector, is to use the Adams-Bashforth method of the same order as the predictor. Thus, we might be better off in the trapezoid rule using the AB2 method for the predictor...but that's not what the book does. This kind of predictor/corrector method is called an Adams-Bashforth-Moulton method or ABM.

COMPUTATION: Try out a simple PECE version of an Adams-Bashforth-Moulton method of order 3. Copy the code abm3_pece.m from the web page. Test the method on the bump_ode.m problem, using initial condition y(0)=1, to x = 5.0, with stepsize 0.25.

        X       Exact Y        ABM3__PECE

        0.0       1.000        1.0
        1.0       0.500        __________ 
        2.0       0.333        __________ 
        3.0       0.250        __________ 
        4.0       0.200        __________
        5.0       0.167        __________ 
      

BDF Methods

BDF methods (backward differentiation formulas) are similar to Adams-Moulton methods, in that they are implicit, and intended to be stable in the face of stiff problems. While the backward Euler method qualifies as the first member of both of these families, the two families differ in that the Adams-Moulton methods look back at values of the derivative y' at old points, while the BDF methods look back at values of y itself. (This makes a BDF method easier to use with a variable step size, which we haven't even started to think about yet!)

The second member of the BDF family, "BDF2", can be summarized as

x_k+1 = x_k + dx
y_k+1 = y_k + 1/3 * y_k - 1/3 * y_k-1 + dx * (2/3) * f(x_k+1,y_k+1)

BDF2 is implicit, and so there are the same problems of getting a starting point, and deciding how well to solve the nonlinear equation at each step, as with the backward Euler method. A new problem occurs, however, because we can't use BDF2 until we have two old solution values. (This is the same problem we had with the explicit multistep methods, like AB2).

When we had this startup problem with the Adams-Bashforth methods, we simply used a Runge-Kutta method. But with a stiff problem, even one RK2 step may get us into trouble! So for now, let's decide only to use implicit methods. Then we're stuck with using a first step of BDF1 before we can switch to our desired method of BDF2. (If stability is really a problem, and we want to use a BDF method of order N, but we don't want to lose accuracy, our first N-1 steps could be very small BDF1 steps, assuming that our BDF-N method can handle a variable stepsize).

COMPUTATION: copy bdf2_newton.m from the web page. Try your code out on the stiff ODE, going from x = 0 to 2, with initial condition y = 1. Compare your results with the results from be_newton.m. You should expect better accuracy.

Direction Field Plots

One way of looking at a differential equation is to plot its direction field. At any point (x,y), we can plot a little arrow in the direction that the differential equation is telling us to go. This is something like the wind direction. What is that direction? Think of how the Euler method works: if you take a step of dx in the x direction, then you must move (dy/dx)*dx in the y direction. In other words, the direction vector is some multiple of (dx,(dy/dx)*dx).

So to make a direction field plot of, say, our stiff ODE, we might want to pick the points (x,y) in the unit square that are spaced dx = 1/10 apart. The MATLAB meshgrid command will help us do that. Then we need to evaluate the direction vector (u,v). The u vector will be all 1's (use the MATLAB ones function), and the v vector comes from our right hand side function. We probably want to use a scale factor of about dx to keep the arrows from getting muddled. Finally, we use the special MATLAB command quiver to display the vector plot. (I don't think "quiver" is available in OCTAVE!). Here's how to plot the direction field for the [0,2] by [-1,1] range of (x,y) values:

dx = 0.1;
[xp,yp] = meshgrid ( 0:dx:2, -1:dx:1 );
u = dx * ones ( size ( xp ) );
v = dx * stiff_ode ( xp, yp );
quiver ( xp, yp, u, v )

u = 2 * u;
v = 2 * v;
quiver ( xp, yp, u, v )

axis equal

GRAPHICS: Use the above commands to make a plot of the direction field of the stiff ODE. Compute the solution of the ODE from 0 to 2, using 10 steps, and the bdf2_newton code. Issue the hold on command, and then plot your solution on top of the direction field. Can you see why this is a hard problem to solve?

If the differential equation is higher order, then it's more difficult to figure out a sensible way of plotting a direction field. Consider, for instance, how you would try to plot the direction field for the predator-prey problem! However, that's a nice problem for phase plane plots!

Summary

It can be hard to keep track of the different ODE methods and their properties. Here is a summary of the four important families of numerical ODE solvers we have discussed in the labs and in the class.

• Runge-Kutta methods are explicit, and require y_k and y'_k, and values of y and y' at auxilliary points between y_k and y_k+1. (RK1 = Euler's method, RK2 = Euler Halfstep, RK3, ...)
• Adams-Bashforth methods are explicit, and require y_k, and values of y' at previously computed points k, k-1, ... (AB1 = Euler's method, AB2, ...)
• Adams-Moulton methods are implicit, and require y_k, an estimate of y_k+1, and values of y' at points k+1, k, ... (AM1 = Backward Euler Method, AM2 = trapezoid method, AM3, ...)
• BDF methods are implicit, and require an estimate of y_k+1, the value of y'_k+1, and values of y at points k, k-1, .... (BDF1 = Backward Euler Method, BDF2, ...)

The midpoint method doesn't really fit into any of these categories. It's an explicit method, but it's not an Adams-Bashforth method, because it uses an older solution value, y_k-1.

Assignment

Compute the numerical solution of the model ODE, from x = 0.0 to x = 2.0, using the trapezoid PECECE method, with the given stepsizes. Record the value of the numerical solution at x = 2.0, and the difference between y(2.0) and exp(2.0).

        dx            y(2.0)          y(2.0) - exp(2.0)

        0.2        __________          __________ 
        0.1        __________          __________ 
        0.05       __________          __________ 
        0.025      __________          __________ 
      

• the results of this computation (eight numbers);
• your guess for the rate of convergence of this method. (As dx descreases, does the error at x=2.0 behave like a constant, or dx, or dx² or dx³...?)