Introduction

The explicit methods that we discussed last time are well suited to
handling a large class of ODE's. These methods perform poorly,
however, for a class of ``stiff'' problems, that occur all too
frequently in applications. We will examine *implicit methods*
that are suitable for such problems. We will find that the
implementation of an implicit method has a complication we didn't
see with the explicit method: solving the nonlinear equation that
generally arises.

The term ``stiff'' as applied to ODE's does not have a precise definition. Loosely, it means that there is a very wide range between the most rapid and least rapid (with ) changes in solution components. A reasonably good rule of thumb is that if Runge-Kutta or Adams-Bashforth or other similar methods require much smaller steps than you would expect, based on your expectations for the solution accuracy, then the system is probably stiff. The term itself arose in the context of solution for the motion of a stiff spring when it is being driven at a modest frequency of oscillation. The natural modes of the spring have relatively high frequencies, and the ODE solver must be accurate enough to resolve these high frequencies, despite the low driving frequency.

This lab will take three sessions. If you print this lab, you may prefer to use the pdf version.

Stiff Systems

I've warned you that there are problems that 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 system

whose solution is

for a constant. For the initial condition at , the constant can easily be seen to be

- For large and except very near , the solution behaves as if it were approximately , which has a derivative of modest size.
- Small deviations from the curve (because of initial conditions or numerical errors) cause the solution to have very large derivatives.

We will be using this stiff differential equation in the exercises below, and in the next section you will see a graphical display of stiffness for this equation.

Direction Field Plots

One way of looking at a differential equation is to plot its ``direction field.'' At any point , we can plot a little arrow equal to the slope of the solution at . This is effectively the direction that the differential equation is ``telling us to go,'' sort of the ``wind direction.'' How can you find ? If the differential equation is , and represents , then , or for a sensibly-chosen value . Matlab has some built-in functions to generate this kind of plot.

**Exercise 1**: In this exercise, you will see a graphical illustration of why a differential equation is ``stiff.''- Copy the following lines into a file called
`stiff2_ode.m`:function f = stiff2_ode ( x, y ) % f = stiff2_ode ( x, y ) % computes the right side of the ODE % dy/dx=f(x,y)=lambda*(-y+sin(x)) for lambda = 2 % x is independent variable % y is dependent variable % output, f is the value of f(x,y). LAMBDA=2; f = LAMBDA * ( -y + sin(x) );

- Also copy the following lines into a second file
called
`stiff2_solution.m`.function y = stiff2_solution ( x ) % y = stiff2_solution ( x ) % computes the solution of the ODE % dy/dx=f(x,y)=lambda*(-y+sin(x)) for lambda = 2 % and initial condition y=0 at x=0 % x is the independent variable % y is the solution value LAMBDA=2; y = (LAMBDA^2/(1+LAMBDA^2))*sin(x) + ... (LAMBDA /(1+LAMBDA^2))*(exp(-LAMBDA*x)-cos(x));

- Create new versions of these files, two using
`LAMBDA=50`and named`stiff50_ode.m`and`stiff50_solution.m`, and two using`LAMBDA=10000`and named`stiff10000_ode.m`and`stiff10000_solution.m`. - The solution to our stiff ODE is roughly , so we are
interested in values of between 0 and , and values of
between -1 and 1. The Matlab
`meshgrid`command is designed for that (it is kind of a two-dimensional`linspace`). To evaluate the direction vector , will be all 1's (use the Matlab`ones`function), and comes from our right hand side function. Finally, we use the special Matlab command`quiver`to display the vector plot. Use the following commands to plot the direction field for the range of`(x,y)`values:h = 0.1; % mesh size [x,y] = meshgrid ( 0:h:2*pi, -1:h:1 ); px = ones ( size ( x ) ); py = stiff2_ode ( x, y ); quiver ( x, y, px, py ) axis equal %this command makes equal x and y scaling

- Finally, to see how the direction
field relates to the approximate solution, plot the function
on the same frame.
hold on x1=(0:h:2*pi); y1=stiff2_solution(x1); plot(x1,y1,'r') % solution will come out red.

Send me a copy of your plot (`print -djpeg ex1.jpg`, or File Export and choose jpeg).

*very rapidly*toward , and then follow that curve as it varies*slowly*. Some numerical methods overshoot the solution curve in their enthusiasm to reach it and avoiding this overshoot characterizes methods for stiff systems.- Copy the following lines into a file called

**Exercise 2**: In this exercise, you will be seeing a numerical illustration of why an ODE is ``stiff.''- Using
`LAMBDA=10000`to represent a stiff equation, how many points would be needed in the interval [0,2*pi] to make a pleasing plot of`stiff10000_solution(x)`?x=linspace(0,2*pi,10); % try 10, 20, 40, etc. plot(x,stiff10000_solution(x))

Try it, but I think you will agree with me that it takes about 40 evenly-spaced points to make a reasonable curve. Do not send me copies of your plots. - Now, choose either of
`euler.m`, or`rk3.m`from last lab and attempt to solve`stiff10000_ode.m`for 40 steps over the interval`[0,2*pi]`, starting from the initial condition`y=0`at`x=0`. Your solutions blow up (are unstable), don't they? (If you are plotting the solutions, look at the scale!) Multiply the number of steps by 10 repeatedly until you get something reasonable, i.e., try 40, 400, 4000, etc. steps. How many steps does it take to get a reasonable plot? It takes many more steps to get a reasonable solution than you would expect, based on solution accuracy, and this behavior is characteristic of stiff systems.

- Using

The Backward Euler Method

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

You might think there is no difference between this method and Euler's
method. But look carefully-this is *not* a ``recipe,'' the way
some formulas are. It is an equation that must be solved for
, *i.e.*, the equation defining is implicit.
It turns out that implicit methods are much better suited to stiff
ODE's than explicit methods.

If we plan to use Backward Euler to solve our stiff ode equation,
we need to address the method of solution of the implicit equation
that arises. Before addressing this issue in general, we can treat
the special case:

This equation can be solved for easily enough! The solution is

In the following exercise, we will write a version of the backward Euler method that implements this solution, and only this solution. Later, we will look at more general cases.

**Exercise 3**:- Write a function m-file called
`back_euler_lam.m`with signature linefunction [x,y]=back_euler_lam(xRange,yInitial,numSteps) % [x,y]=back_euler_lam(xRange,yInitial,numSteps) % comments % your name and the date.

that implements the above algorithm (with`10000`). You may find it helpful to start out from a copy of your`euler.m`file, but if you do, be sure to realize that this function does not have the input parameter`f`that`euler.m`does because the function is built-in. - Be sure to include comments following the signature line.
- Test
`back_euler_lam`by solving the system (1) for 40 steps over the interval , starting from the initial condition at`x=0`. Your solution should not blow up and its plot should look reasonable. Please include this plot with your summary. - For the purpose of checking your work, the first few values of y are
y=[ 0.0, 0.156334939127, 0.308919855800, 0.453898203657, ...].

- Write a function m-file called

In the next exercise, we will compare backward Euler with
(forward) Euler for accuracy. Of course, backward Euler gives results when the
stepsize is large and Euler does not, but we are curious about the
case that there are enough steps to get answers. In addition, we
would like to compare the computed and exact solutions over the whole
range, and we will use the Matlab `norm`
function to do it. Because it would require too many steps to run
this test with `1.e4`, you will be using
`50`.

**Exercise 4**:- Change the value to
`LAMBDA=50`in`back_euler_lam`. - Compute the error in the solution versus the exact solution as
abs( y(end) - stiff50_solution(x(end)) )

Compute the ratios of the errors for each value of`numSteps`divided by the error for the succeeding value of`numSteps`. Use`back_euler_lam`to fill in the following table, over the interval`[0,2*pi]`and starting from`yInit=0`.stiff50 numSteps back_euler_lam error ratio 40 ____________________ __________ 80 ____________________ __________ 160 ____________________ __________ 320 ____________________ __________ 640 ____________________ __________ 1280 ____________________ __________ 2560 ____________________

- Based on the ratios in the table, estimate the order of
accuracy of the method,
*i.e.,*estimate the exponent in the error estimate . is an integer in this case. - Repeat this experiment using
`euler`, and fill in the following table. Note that the errors using`euler`end up being about the same size as those using`back_euler_lam`, once`euler`starts to behave.stiff50 Error comparison numSteps back_euler_lam euler 40 ___________________ _________________ 80 ___________________ _________________ 160 ___________________ _________________ 320 ___________________ _________________ 640 ___________________ _________________ 1280 ___________________ _________________ 2560 ___________________ _________________

(In filling out this table, please include at least four significant figures in your numbers. If necessary, use`format long`to get extra precision printed.) - Compare the rate of convergence using
`euler`with that using`back_euler_lam`.

- Change the value to

Newton's method

You should be convinced that implicit methods are worth while.
How, then, can the resulting implicit equation (usually it is
nonlinear) be solved? The Newton (or Newton-Raphson) method is a good choice.
(We saw Newton's method last semester in Math 2070,
Lab 4,
you can also find it in Iserles, Section 7.2,
or in a Wikipedia article
`http://en.wikipedia.org/wiki/Newton's_method`.)
Briefly, Newton's method is a way to solve equations by successive iterations.
It requires a reasonably good starting point and it requires that the
equation to be solved be differentiable. Suppose you want to solve
the nonlinear (vector) equation

and you know a good starting place . The following iteration is called Newton iteration.

and it (usually) converges quite rapidly. (In the case that is a vector, the partial derivative must be interpeted as the Jacobian matrix, whose components are .) The superscript is used here to distinguish it from the step number,

There is an easy way to remember the formula for Newton's method. Write the finite difference formula for the derivative as

On each step, the backward Euler method requires a solution of the equation

When we have found a satisfactorily approximate solution to (5), then we take and proceed with the backward Euler method. We think of each of the values as successive corrections to .

I call your attention to a difficulty that has arisen. For the Euler, Adams-Bashforth and Runge-Kutta methods, we only needed a function that computed the right side of the differential equation. In order to carry out the Newton iteration, however, we will also a function that computes the partial derivative of the right side with respect to . Last semester we assumed that the derivative of the function was returned as a second variable, and we will use the same convention here.

To summarize, on each time step, start off the iteration by predicting the solution using an explict method. Forward Euler is appropriate in this case. Then use Newton's method to generate successive correction steps. In the following exercise, you will be implementing this method.

**Exercise 5**:- Change your
`stiff10000_ode.m`so that it has the signaturefunction [f, dfdy]=stiff10000_ode(x,y) % [f, dfdy]=stiff10000_ode(x,y) % comments % your name and the date

and computes the partial derivative of`stiff10000_ode`with respect to`y`. To do this, write the derivative of the formula for`stiff10000_ode`out by hand, then program it. - Similarly, change the signature for
`stiff50_ode`. There is no need to change`stiff2_ode`because you won't be using it again. - Copy the following function to a
file named
`back_euler.m`using cut-and-paste.function [x,y]=back_euler(f,xRange,yInitial,numSteps) % [x,y]=back_euler(f,xRange,yInitial,numSteps) computes % the solution to an ODE by the backward Euler method % % xRange is a two dimensional vector of beginning and % final values for x % yInitial is a column vector for the initial value of y % numSteps is the number of evenly-spaced steps to divide % up the interval xRange % x is a column vector of selected values for the % independent variable % y is a matrix. The k-th row of y is the transpose of % the approximate solution at x(k) % your name and the date % force x to be a column vector x=zeros(numSteps+1,1); x(1) = xRange(1); h = ( xRange(2) - xRange(1) ) / numSteps; y(1,:) = transpose(yInitial); for k = 1 : numSteps x(k+1) = x(k) + h; Y = transpose(y(k,:)) + ... h * feval( f, x(k), transpose(y(k,:))); [Y,isConverged]= ... newton4euler(f,x(k+1),transpose(y(k,:)),Y,h); if ~ isConverged error(['back_euler failed to converge at step ', ... num2str(k)]) end y(k+1,:) = transpose(Y); end

- Copy the following function to a file named
`newton4euler.m`(Read this name as ``Newton for Euler.'')function [Y,isConverged]=newton4euler(f,x,ytranspose,Y,h) % [Y,isConverged]=newton4euler(f,x,ytranspose,Y,h) % special function to evaluate Newton's method for back_euler % your name and the date TOL = 1.e-6; MAXITS = 10; isConverged= (1==0); % starts out FALSE for n=1:MAXITS [fValue fPartial] = feval( f, x, Y); F = ytranspose + h * fValue - Y; dFdY = h * fPartial - eye(length(Y)); increment=dFdY\F; Y = Y - increment; if norm(increment,inf) < TOL*norm(Y,inf) isConverged= (1==1); % turns TRUE here return end end

- Insert the following comments into either
`back_euler.m`or`newton4euler.m`where they belong. Some comments belong in both files. Include copies of both files when you send your summary to me.% f is a function whose signature is ??? % TOL = ??? (explain in words what TOL is used for) % MAXITS = ??? (explain in words what MAXITS is used for) % The Matlab "\" symbol means ??? % The Matlab function "eye" is used for ??? % The following for loop performs the spatial stepping (on x) % The following statement computes the initial guess for Newton

- The equation (4) includes the quantities , , and . What are the names of the variables in the code representing these mathematical quantities?
- Insert a comment in the code indicating which lines implement Equation (4).
- Insert a comment in the code indicating which lines implement
Equation (5). (
**Note:**This line is specific to the implicit Euler method, and will have to be changed when the method were changed.) - In the case that
`length(Y)>1`, is`Y`a row vector or a column vector? - If
`f='stiff10000_ode'`,`x=1.0`,`y=3.0`,`h=0.1`, and the initial guess for`Y=1`, write out by hand the (linear) equation that`newton4euler`solves. What is the value of the solution of this linear equation? - Verify that
`newton4euler`is correct by showing it yields the correct solution in the case that`f='stiff10000_ode'`,`x=1.0`,`y=3.0`,`h=0.1`, and the initial guess for`Y=1`. (If you discover that the Newton iteration fails to converge, you probably have the derivative in`stiff10000_ode`wrong.) - In order to check that everything is programmed correctly,
solve the same problem you solved for Exercise 3 using
`back_euler_lam`except use`back_euler`. Solve the ODE using`stiff10000_ode`, on the interval , with initial value 0, for 40 steps. Be sure that`LAMBDA=10000`in`back_euler_lam.m`You can compare all 40 values by comparing the two solutions generated by`back_euler_lam`and`back_euler`. The two solutions should agree to ten or more significant digits. (Testing code by comparison with previous code is called ``regression'' testing.)

- Change your

One simple nonlinear equation with quite complicated behavior is van der Pol's equation. This equation is written

As you know, van der Pol's equation can be written as a system of
first-order ODEs in the following way.

where and .

The matrix generalization of the partial derivative in Newton's method is the Jacobian matrix:

**Exercise 6**: In this exercise you will be solving van der Pol's equation with using`back_euler`.- Write an m-file to compute the right side of the van der Pol
system (6) and its Jacobian (7).
The m-file should be
named
`vanderpol_ode.m`and should have the signaturefunction [f, J]=vanderpol_ode(x,y) % [f, J]=vanderpol_ode(x,y) % more comments % your name and the date if length(y) ~=2 error('vanderpol_ode: y must be a vector of length 2!') end a=10; f = ??? df1dy1 = 0; df1dy2 = ??? df2dy1 = ??? df2dy2 = -a*(y(1)^2-1); J=[df1dy1 df1dy2 df2dy1 df2dy2];

be sure that`f`is a*column*vector and that the value of the parameter . (**Note:**If you get a message saying that convergence failed, you probably have an error in your calculation of the derivatives in`J`. Fix it.) - Solve the van der Pol system starting from
`y=[0;0]`over the interval from`x=0`to`x=2`using 40 intervals using`euler`and also using`back_euler`. You can see by plotting the solutions that`euler`does not give a good solution while`back_euler`does. Send me plots of both solutions with your summary. - Solve the system again, this time using 640 intervals. This
time both methods should remain stable and the answers should be
close. Please include plots of both solutions on a single frame with
your summary. (Use
`hold on`and`hold off`.)The reason that it takes so many intervals to achieve reasonable accuracy is because the the solution remains relatively flat for ``a long time'' before ``taking off.'' Small errors in the flat portion translate into fairly large errors in the location of the peak. The next peak (the solution is asymptotically periodic) occurs near

`x=17`and requires over 10,000 intervals to predict accurately. (At least twice the 320 intervals per unit of used here for the first.) More sophisticated algorithms use*variable*-sized steps and can solve such problems much more efficiently.

- Write an m-file to compute the right side of the van der Pol
system (6) and its Jacobian (7).
The m-file should be
named

The Trapezoid Method

You have looked at the (forward) Euler and backward Euler methods. These methods are simple and involve only function or Jacobian values at the beginning and end of a step. They are both first order, though. It turns out that the trapezoid method also involves only values at the beginning and end of a step and is second order accurate, a substantial improvement. This method is also called ``Crank-Nicolson,'' especially when it is used in the context of partial differential equations. As you will see, this method is appropriate only for mildly stiff systems.

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

from which you can see that this is also an ``implicit'' formula. The backward Euler and Trapezoid methods are the first two members of the ``Adams-Moulton'' family of ODE solvers.

In the exercise below, you will write a version of the trapezoid
method using Newton's method to solve the per-timestep equation, just
as with `back_euler`. As you will see in later exercises, the
trapezoid method is not so appropriate when the equation gets very
stiff, and Newton's method is overkill when the system is not stiff.
The method can be successfully implemented using an approximate
Jacobian or by computing the Jacobian only occasionally, resulting in
greater computational efficiency. We are not going to pursue this
alternative in this lab, however.

**Exercise 7**:- Make a copy of your
`back_euler.m`file, naming the copy`trapezoid.m`. Modify it to have the signature:function [ x, y ] = trapezoid ( f, xRange, yInitial, numSteps ) % comments

- Add appropriate comments below the signature line.
- Replace the lines
Y = transpose(y(k,:)) + ... h * feval( f, x(k), transpose(y(k,:))); [Y,isConverged]= ... newton4euler(f,x(k+1),transpose(y(k,:)),Y,h);

with the new linesf_xk_yk= feval( f, x(k), transpose(y(k,:))); Y = transpose(y(k,:)) + h * f_xk_yk; [Y,isConverged]=newton4trapezoid( ... f,x(k+1),transpose(y(k,:)),Y,h,f_xk_yk);

- Make a copy of your
`newton4euler.m`, naming the copy`newton4trapezoid.m`. - Change the signature of
`newton4trapezoid`tofunction [Y,isConverged]=newton4trapezoid(f,x,ytranspose,Y,h,f_xk_yk)

- Increase the iteration limit
`MAXITS=30;` - In order to implement the trapezoid method, you need to
write the function that appears in (4)
and (5) so it is valid for the
trapezoid rule. It is currently written for the backward
Euler method (5). To do this, consider
(8), repeated here
right sideModify
`newton4trapezoid.m`to solve this function. Do not forget to modify the Jacobian`dFdY`to reflect the new function`F`. Remember, if you have the Jacobian wrong, the Newton iteration will fail. - Test your
`newton4trapezoid`code for the case`f='stiff50_ode'`,`x=1.0`,`y=3.0`,`h=0.1`, =`f_xk_yk=1`, and the initial guess for`Y=1`. - Write out by hand the (linear) equation that
`newton4trapezoid`solves in the case`f='stiff50_ode'`,`x=1.0`,`y=3.0`,`h=0.1`,`f_xk_yk=1`. Does the solution to this equation agree with the one from`newton4trapezoid`? - Test the accuracy of the trapezoid method by computing the numerical
solution of
`stiff50_ode.m`(starting from over the interval , and fill in the following table. Compute ``error'' as the difference between the computed and exact values at .stiff50 numSteps trapezoid error ratio 10 ____________________ __________ 20 ____________________ __________ 40 ____________________ __________ 80 ____________________ __________ 160 ____________________ __________ 320 ____________________

- Are your results consistent with the theoretical convergence rate of the trapezoid rule? If not, you have a bug you need to fix.

- Make a copy of your

In the following exercise, you are going to see that the difference in accuracy between the trapezoid method and the backward Euler method can improve the solution of the vander Pol equation.

**Exercise 8**:- Use
`trapezoid.m`to solve the van der Pol equation starting from the value`[0;0]`on the interval`[0,10]`using 100, 200, 400, and 800 steps, and plot all four solutions on one plot. You should see that the solution for 400 steps is pretty much the same as the one for 800 steps, and the others are different. We will assume that the final two approximate solutions represent the correct solution. Please send me this plot with your summary. Hint: you can generate this plot easily with the following sequence of commands.hold off [x,y]=trapezoid('vanderpol_ode',[0,10],[0;0],100);plot(x,y) hold on [x,y]=trapezoid('vanderpol_ode',[0,10],[0;0],200);plot(x,y) [x,y]=trapezoid('vanderpol_ode',[0,10],[0;0],400);plot(x,y) [x,y]=trapezoid('vanderpol_ode',[0,10],[0;0],800);plot(x,y) hold off

- Use
`back_euler.m`to solve the van der Pol equation starting from the value`[0;0]`on the interval`[0,10]`using 400, 800, 3200, and 12800 steps, and plot all four solutions on one plot. The larger number of points may take a minute. You should see a progression of increasing accuracy in the second ``pulse.'' Please send me this plot with your summary.

- Use

The trapezoid method is unconditionally stable, but this fact does not mean that it is good for very stiff systems, however. In the following exercise, you will apply the trapezoid method to a very stiff system so you will see that numerical errors arising from the initial rapid transient persist when using the trapezoid rule but not for backwards Euler.

**Exercise 9**:- Solve the
`stiff10000_ode`system starting from`yInitial=0.1`on the interval`[0,10]`using 100 steps. Note that the initial condition is not zero. Use both the trapezoid method and backwards Euler, and plot both solutions on the same plot. Send me this plot with your summary file. - Use the trapezoid method to solve the same case using 200, 400, 800, and 1600 steps. Plot each solution on its own plot. You should see that the effect of the initial condition is not easily eliminated. Please include the 1600 step case when you send your summary file.

- Solve the

Matlab ODE solvers

Matlab has a number of built-in ODE solvers. These include:

Matlab ODE solvers | |

ode23 | non-stiff, low order |

ode113 | non-stiff, variable order |

ode15s | stiff, variable order, includes DAE |

ode23s | stiff, low order |

ode23t | trapezoid rule |

ode23tb | stiff, low order |

ode45 | non-stiff, medium order (Runge-Kutta) |

The ODE solvers we have written have four parameters, the function name, the solution interval, the initial value, and the number of steps. The Matlab solvers use a good adaptive stepping algorithm, so there is no need for the fourth parameter.

**Exercise 10**:- You can watch a solution evolve if you call the solver
without any output variables. Use
`ode45`to solve the same van der Pol problem we have been solving on the interval`[0,10]`, starting from`[0;0]`. Use the following Matlab command:ode45('vanderpol_ode',[0,10],[0;0])

Please include this plot with your summary. - You can see the difference between a stiff solver and a non-stiff
solver on this moderately stiff problem by comparing
`ode45`with`ode15s`. (These are my personal favorites.) You can have both plots on the screen at once using the following commandsfigure(1) ode45('vanderpol_ode',[0,10],[0;0]) figure(2) ode15s('vanderpol_ode',[0,10],[0;0])

You should be able to see that the step density (represented by the little circles on the curves) is somewhat less for`ode15s`. The difference would be much larger if you were to try a stiffer problem, such as`a=1000`in`vanderpol_ode.m`. You do not need to send me these plots. - If you wish to examine the solution in more detail, or
manipulate it, you need the solution values, not a picture. You
can get the solution values with commands similar to the following:
[x,y]=ode45('vanderpol_ode',[0,10],[0;0]);

If you wish, you can plot the solution, compute its error, etc. For this solution, what is the value of at`x=10`? How many steps did it take? (The length of`x`is one more than the number of steps.) - Suppose you want a very accurate solution. The default
tolerance is .001 relative accuracy in Matlab, but suppose you
want a relative accuracy of 1.e-8? There is an extra variable
to provide options to the solver. It works in the following manner:
myoptions=odeset('RelTol',1.e-8) [x,y]=ode45('vanderpol_ode',[0,10],[0;0],myoptions);

(The semicolon is left off the odeset call so you can see some of the other options and their defaults. Use`help odeset`for more detail.) How many steps did it take this time? What is the value of at`x=10`?

- You can watch a solution evolve if you call the solver
without any output variables. Use

Back to MATH2071 page.

Mike Sussman 2009-02-04