|Stiff Systems||Exercise 2|
|Direction Field Plots||Exercise 3|
|The Backward Euler Method||Exercise 4|
|Newton's method||Exercise 5|
|The Trapezoid Method||Exercise 6|
|Matlab ODE solvers||Exercise 7|
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.
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
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.
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.
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) );
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));
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
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).
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.
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
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:
function [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.
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.
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 ____________________
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.)
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
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
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.
function [f, dfdy]=stiff10000_ode(x,y) % [f, dfdy]=stiff10000_ode(x,y) % comments % your name and the dateand 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.
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
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
% 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
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.
The matrix generalization of the partial derivative in Newton's method is the Jacobian matrix:
function [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.)
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.
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:
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.
function [ x, y ] = trapezoid ( f, xRange, yInitial, numSteps ) % comments
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 lines
f_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);
stiff50 numSteps trapezoid error ratio 10 ____________________ __________ 20 ____________________ __________ 40 ____________________ __________ 80 ____________________ __________ 160 ____________________ __________ 320 ____________________
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.
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
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.
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|
|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.
ode45('vanderpol_ode',[0,10],[0;0])Please include this plot with your summary.
figure(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.
[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.)
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?
Back to MATH2071 page.