Math 128a - Spring 2002 - Notes on ODEs Solving Ordinary Differential Equations Read sections 8.0 to 8.6 in the textbook. Initial Value Problem (IVP): x'(t) = f(x(t),t) with initial condition (IC) x(t0) = x0 (note that x and f could be vectors of length n) In this case we want an approximate solution for t > t0, up to some t_final, if it exists. Ex: Most models of the physical world involve differential equations ordinary (ODE) or partial (PDE) Physics: Everyday mechanics: Newton's Laws (ODE) F = ma = m * x'' Quantum mechanics: Schroedinger's equation (PDE) Fluids: Navier-Stokes equation (PDE) Electricity and magnetism: Maxwell's equation (PDE) Electrical Circuits: Kirchoff's, Faraday's, Ohm's ... Laws (ODE) Chemistry: Chemical reactions (ODE) Biology: Population dynamics (ODE) etc etc etc Ex: x'(t) = c*(b(t) - x(t)), c>0, b(t) a given function = 0 if x(t) = b(t) > 0 if x(t) < b(t) < 0 if x(t) > b(t) i.e. x(t) wants to get closer to b(t) b(t) = sin(t), c = 1, t0=0, t_final = 20, x(0) = 2 hold off, clf, [t,y]=ode23(@chase,[0 20],2); plot(t,y), hold on, plot(t,b(t),'r') function f = chase(t,y) f = c*(b(t)-y); b(t) = sqrt(t)*sin(t), c = 5, t0=0, t_final = 20, x(0) = 2 b(t) = sqrt(t)*sin(5*t), c = 5, t0=0, t_final = 20, x(0) = 2 b(t) = sqrt(t)*sin(5*t), c = 20, t0=0, t_final = 20, x(0) = 2 Qualitatively: x(t) chases b(t), the faster the large c is Ex: Solve x'(t) = c*(b(t) - x(t)) = -c*x(t) + c*b(t) exactly: (review of calculus, works for c = c(t) too) rewrite as: x'(t) + c*x(t) = c*b(t) multiply by exp(c*t), simplify: exp(c*t)*x'(t) + c*exp(c*t)*x(t) = c*exp(c*t)*b(t) = d/dt ( exp(c*t)*x(t)) integrate from 0 to s: integral from 0 to s of d/dt (exp(c*t)*x(t) dt = exp(c*s)*x(s) - exp(0)*x(0) = integral from 0 to s of c*exp(c*t)*b(t) dt solve for x(s) x(s) = exp(-c*s)*x(0) + integral from 0 to s of c*exp(c*(t-s))*b(t) dt Pictures explain this behavior more simply ("x(t) chases b(t)") than the formula If, as above, the first derivative is the highest derivative of x(t) appearing in the ODE, then the ODE is first order. If the p-th derivative is the highest one appearing, it is p-th order. It is possible to convert any p-th order ODE into a first order ODE, with x of larger dimension. Ex: Consider Newton's Law F = ma, which governs mechanical systems, for the simple mass-spring system: ----- | \ / | m where m is a mass that can go up and down (without friction), x is the vertical position of m (x=0 is the rest position, where m hangs without moving, and x>0 means is higher, x<0 lower) and the spring has "spring constant k"; this means that moving m down a distance x>0 up from its resting state creates a force F=-k*x pushing m down, and moving m down to x<0 creates a force F=-k*x pullling m up. Thus F=ma becomes -k*x(t) = m*x''(t) To define the state of the system at t{0}=0, we need to specify the position x(0)=x0 and velocity x'(0)=v0 of m as ICs. We convert this into a first order system of 2 ODEs by creating a new vector function X(t) = [ x(t) ] [ x'(t) ] which satifies the first order ODE X'(t) = [ x'(t) ] = [ x'(t) ] = [ X_2(t) ] [ x''(t) ] [ -k/m*x(t) ] [ -k/m*X_1(t) ] = [ 0 1 ] * X(t) [ -k/m 0 ] with initial conditions X(t{0}) = [ x0 ] [ v0 ] Solve numerically: hold off, clf, [t,X]=ode23(@weight1,[0 30],[1;0]); plot(t,X(:,1),'b',t,X(:,2),'r'), grid function f = weight1(t,X) k=1; m=2; f = [ X(2) ; -k/m * X(1) ]; What does it look like? What is relationship of blue (position) and red (velocity)? What happens if you double k? We can also write down the solution of this ODE explicitly as follows: x''(t) = -k/m*x(t) has two general solutions, x(t) = sin(f*t) and cos(f*t), where f = sqrt(k/m). Thus x(t) = a*sin(f*t) + b*cos(f*t) is a solution for any constants a and b. We will choose a and b to satisfy the initical conditions: a*sin(f*0) + b*cos(f*0) = b = x(0) = x0 f*a*cos(f*0) - f*b*sin(f*0) = f*a = x'(0) = v0 so x(t) = v0/f*sin(f*t) + x0*cos(f*t) Ex: Newton's Laws with 2 masses: ----- | \ spring constant k1 / force on m1 = -k1*x1 from top spring | +k2*(x2-x1) from bottom spring m1 rest position x1 so m1*x1'' = -k1*x1 + k2*(x2-x1) | \ spring constant k2 / | force on m2 = -k2*x2 from bottom spring m2 rest position x2 so m2*x2'' = -k2*x2 Get 2 2nd order equations => 4 first order equations solve numerically on homework using matlab's built in functions Most cars are designed by doing such a simulation of what happens when a car crashes, and what happens to the occupants, in addition to actual crash tests with dummies (the numerical simulations are much faster and cheaper than actual car crashes!) For design of turbine blades in jet engines, simulating the effect of birds hitting them is used by the FAA to certify the safety of jet engine design for commercial airliners. The ship certification industry uses simulation to decide on whether to insure ocean-going ships. The FDA requires simulations of any body implants (eg heart valves). Boundary Value Problems (BVP) (we won't have time for methods but you should know about them) x''(t) = f(x'(t),x(t),t) with boundary conditions (BCs) x(t{0}) = x{0} and x(t{1}) = x{1} In this case we want x(t) for t{0} < t < t{1} This describes heat distribution in a bar, the shape of a vibrating string, and many other things: Ex: a thin bar stretches of length 1 has conductivity k(x) (0 <= x <= 1), and the ends of the bar are held at temperatures T(0) = T0 and T(1) = T1, then the temperature T(x) along the bar satisfies k(x)*T''(x) + k'(x)*T'(x) = 0, T(0) = T0 and T(1) = T1. If k(x) = k is constant, this simplifies to k*T''(x) = 0, or T(x) = a*x+b for some a and b. plugging in BCs yields T0 = a*0 + b = b, and T1 = a*1 + b => b = T1-T0, so the temperature changes linearly from one end to the other. Delay differential equation x'(t) = f(x(t),x(t-a),t) with x(t{0}) = x{0} or x'(t) = f(x(t),x(t-a),x(t-b),...,t) with x(t{0}) = x{0} This describes any system where the reaction of some component is delayed (eg any mechanical system controlled by sensors and a computer of finite speed) Differential-Algebra equations (DAE) Suppose x has n components 0 = f(x'(t),x(t),y(t),t) (m equations) 0 = g(x(t),y(t),t) (n-m equations) This is a natural way to describe a "constrained" system, e.g. an object sliding downhill on a surface (then x(t) is the position and velocity of the object, the first equation is "F=ma" for the unconstrained components of x(t), and the second equation says that x(t) is constrained to lie on the surface. The standard equations describing an electric circuit are DAEs. The entire circuit industry in Silicon Valley and elsewhere runs simulations of circuits on their chips before they actually build them to make sure they work first. The software they use for this is called SPICE, which was originally developed at Berkeley. Back to ODES: A large body of excellent software has been developed for all of the above categories of ODEs. You will almost certainly just want to use some of this existing software to solve an ODE that you get in practice, especially if you only need to solve it a few times. Some of this software is in Matlab, but a much more extensive set is in Netlib (www.netlib.org, in C or Fortran). The only exceptions would be 1) if your ODE is especially structured so that you can come up with a much faster method specialized to your problem, and you need to solve your problem many many times, or 2) your problem is so difficult that existing software can't solve it at all, and you need to develop your own method. Overview of Theory of Existence and Uniqueness of ODEs We can't expect to be able to show that a numerical method works well unless we that the true solution exists and is unique Ex 1: Consider x'(t) = [ 0 if t <= 1 , x(0) = 0 [ 1 if t > 1 Can this have a solution? Clearly, x(t) must be constant for t<=1, and this constant must then be x(0) = 0. Next x(t) must be a straight line with slope 1 for t>1. If x'(1) is to exist, so that the ODE is really solvable everywhere, then x(t) must certainly be continuous at 1, so x(t) would have to be x(t) = [ 0 if t <= 1 [ t-1 if t > 1 This function, sometimes called a "ramp" is clearly not differentiable at t=1, so the ODE does not have a solution everywhere. (Sometimes having a solution everywhere except a few points is ok.) This illustrates that we want f(x,t) in x'(t) = f(x(t),t) to at least be a continuous function. Ex 2: Consider x'(t) = f(x) = x^2, x(0) = 1. f(x) is a smooth function. But it is easy to see that the solution x(t) is an increasing function, since its derivative is always positive. We can solve by rewriting dx/dt = x^2 as dx/x^2 = dt or -1/x = t + C or x(t) = 1/(-C-t). Plugging in x(0) = 1 yields C=-1 or x(t) = 1/(1-t) which blows up at t=1. Thus we need f(x,t) to be bounded somehow to expect a solution for large t. Ex 3: Consider x'(t) = x^(1/2), x(0) = 0. It is easy to see that x(t) = t^2/4 is a solution. But for any t0>0, another solution is given by x(t) = [0 if t < t0 [ (t-t0)^2/4 if t>= t0 so the solution is not unique. Note that the slope of f(x) = sqrt(x) is infinite at x=0; this is the problem, and so to guarantee a unique solution (for otherwise proving that a numerical method converges will be hard!) we will assume that f(x) cannot have infinite slope. All these examples motivate the following theorems: Theorem 1 (Existence of solution to x'=f(x,t)). Suppose f(x,t) is continuous on the rectangle R = {(x,t): x0-b <= x <= x0+b and t0-a <= t <= t0+a} with an upper bound of |f(x,t)| <= M on R. Then the ODE has a solution x(t) for t0 - c <= t <= t0 + c where c = min(a,b/M) Ex 4: Consider Ex 2 again. Apply the theorem with f(x) = x^2, x0 = 1 and t0 = 0. f(x) is continuous on R for any a and b, with bound M = (1+b)^2. This means that a solution exists for |t| <= c = min(infinity, b/(1+b)^2) = b/(1+b)^2. It is easy to see that b/(1+b)^2 is maximized for b=1, yielding that a solution exists for |t|<.25. In fact we know that a solution exists for t<1, so this is a little pessimistic; Theorem 1 is too weak to say exactly when and if the solution blows up! Ex 3 shows that we need more than f being continuous to guarantee a unique solution: Def: If there is a constants L >= 0 such that g(x) satisfies |g(x1)-g(x2)| <= L*|x1-x2| for all x1 and x2, then we say that g(x) is Lipschitz continuous, with Lipschitz constant L. This says that if g(x) has a derivative, that derivative is bounded by L. Conversely, if g(x) has a derivative which is bounded in absolute value by L, then g(x) is Lipschitz continuous. But g(x) need not be differentiable (it can have corners, for example). Theorem 2 (Uniqueness of solution to x'=f(x,t)). Suppose f(x,t) is Lipschitz continuous in x for t0-a <= t <= t0+a. Then x'(t)=f(x,t), x(t0) = x0 has a unique solution for t0-a <= t <= t0+a and for any x0. Algorithms for ODEs: Idea for solving x'(t) = f(x(t),t), x(t0) = x0, suppose we have approximated x up to some t, so we know x(t): how to get next value?: approximate x() near t by a polynomial p(), evaluate p(t+h) for a small step h>0 to get x(t+h) Simplest possible polynomial: a straight line through x(t) with slope x'(t) x(t+h) = x(t) + h*x'(t) = x(t) + h*f(x(t),t) This is called Euler's method. Here is the simplest possible Matlab code, which computes an array of values of t and of x, with x(i) being the computed solution at t(i) = t0 + (i-1)*h: t(1) = t0 x(1) = x0 for i=2:n t(i) = t(i-1) + h x(i) = x(i-1) + h*f(x(i-1),t(i-1)) end Here is another way to look at it: Write f(x(t),h) = x'(t) ~ ( x(t+h)-x(t) )/h, and solve for x(t+h). Let's try it, and see how well it works. Ex: First an example where we can evaluate everything by hand: x'(t) = x(t), x(0) = 1, The exact solution is x(t) = exp(t) Euler's method yields x(t+h) = x(t) + h*x'(t) = x(t) + h*x(t) = (1+h)*x(t) or x(0+n*h) = (1+h)^n*x(0) = (1+h)^n Write t=n*h, or x(t) = (1+h)^(t/h) = ((1+h)^(1/h))^t As h -> 0 this approaches e^t, so it converges Ex: x'(t) = sin(t) - x(t) (earlier, we showed that the solution x(t) "chased" sin(t) options = odeset('AbsTol',1e-12,'RelTol',1e-12); % compute answer to within 1e-12 hold off, clf, [t,y]=ode45(@chase,[0 10],4,options); plot(t,y), hold on Compute answer with step h=.1, .01, .001: [T1,Y1]=euler(@chase,[0,10],4,.1); [T01,Y01]=euler(@chase,[0,10],4,.01); [T001,Y001]=euler(@chase,[0,10],4,.001); plot(T1,Y1,'c'), plot(T01,Y01,'m'), plot(T001,Y001,'g'), Look at final error: yfinal = y(length(y)); Y1final = Y1(length(Y1)); Y01final = Y01(length(Y01)); Y001final = Y001(length(Y001)); [yfinal-Y1final; yfinal-Y01final; yfinal-Y001final] Result (may depend on platform): -2.1e-2, -2.1e-3 -2.1e-4 In other words, decreasing the step size h by a factor of 10 decreases the final error by a factor of 10, or that the error is proportional to h. If we wanted error 1e-8, how small would h have to be? How many steps would it take? Two things to do next: 1) Find "higher order methods", where error proportional to O(h^2), O(h^3) or a higher power, so we can get a smaller error for the same h, or use a larger h, take fewer steps, and run faster, or both 2) How to choose h to make the final error as small as desired, but not too small, to save time Approach: 1) figure out error in taking one step of method, 2) see how errors "accumulate" over many steps Let's apply this approach to Euler's method: True x(t+h) = x(t) + h*x'(t) + h^2/2*x''(z) where t <= z <= t+h Euler's x(t+h) = x(t) + h*x'(t) Subtract to get the error from one step: True - Euler = h^2/2*x''(z) = h^2/2*(df/dx (x(z),z) * f(x(z),z) + df/dt (x(z),z)) = O(h^2) We call this one-step error the Local Truncation Error (LTE), because it comes from truncation the Taylor expansion. In the case of Euler's method, we say that the LTE is O(h^2) or second order. Here is one way to motivate (not prove!) why the total error might be O(h): To get to time t_final, we need to take n = t_final/h steps, each of which has LTE O(h^2), so we might expect these errors to add up, for a total error of t_final/h * O(h^2) = O(h), i.e. proportional to h. Later we will formalize this idea. So we could expect that if another method had LTE proportional to h^3, (or more generally h^k), we would expect that the total error would be (1/h) times larger, i.e. O(h^2) (or more generally O(h^{k-1})). This is indeed the case, provided that the polynomial (or other) formula we use to approximate x(t+h) has a reasonable property called "stability" (the most overused term in numerical analysis). We next present a sequence of methods with LTE O(h^k), all of which are stable. Later, we will investigate stability, and when these methods actually converge as h -> 0. Taylor Series Method: Euler's method x(t+h) = x(t) + h*f(x(t),t) = x(t) + h*x'(t) is just the first two terms of the Taylor expansion of x around t, so the LTE is just the error term in the Taylor series, which is proportional to h^2. So we could make the LTE smaller by taking more terms in the Taylor expansion. For example, x(t+h) = x(t) + h*x'(t) + h^2/2*x''(t) + ... + h^k/k!*x^(k)(t) would have LTE O(h^(k+1)). But how would we get the derivatives x''(t) etc to implement this? If f(x(t),t) is given to us by a formula, we can differentiate it using the chain rule: x' = f(x,t) x'' = f_x*x' + f_t where f_x means partial f / partial x (x(t),t) and f_t means partial f / partial t (x(t),t) = f_x*f + f_t x''' = (f_xx*f + f_xt)*f + f_x*(f_x*f + f_t) + f_tx*f + f_tt where f_xx means partial^2 f / partial x^2 (x(t),t) and f_tx means partial^2 f / partial t partial x (x(t),t) etc So a method with LTE = O(h^3) would be x(t+h) = x(t) + h*x'(t) + h^2/2*x''(t) (*) = x + h*f + h^2/2*(f_x*f + f_t) Ex: x'(t) = sin(t) - x(t) (earlier, we showed that the solution x(t) "chased" sin(t) options = odeset('AbsTol',1e-12,'RelTol',1e-12); % make answer accurate hold off, clf, [t,y]=ode45(@chase,[0 10],4,options); plot(t,y), hold on Compute answer with step h=.1, .01, .001 [T1,Y1]=taylor2(@chase2,[0,10],4,.1); [T01,Y01]=taylor2(@chase2,[0,10],4,.01); [T001,Y001]=taylor2(@chase2,[0,10],4,.001); plot(T1,Y1,'c'), plot(T01,Y01,'m'), plot(T001,Y001,'g'), Look at final error: yfinal = y(length(y)); Y1final = Y1(length(Y1)); Y01final = Y01(length(Y01)); Y001final = Y001(length(Y001)); [yfinal-Y1final; yfinal-Y01final; yfinal-Y001final] Result: 4.5e-04 4.5e-06 4.5e-08 Advantages: High order and so high accuracy possible Disadvantages: May be hard to get all derivatives Derivatives may not exist (f need not be differentiable to have a unique solution) Runge-Kutta Methods Runge-Kutta methods use a clever trick to get more terms in Taylor expansion, requiring only that we evaluate f(x,t) at some intermediate points between t and t+h. No formulas for derivative of f(x,t) are needed; all we need is a subroutine to evaluate f(x,t). To get a second order method (i.e. one with LTE = O(h^3), consider a method of the form x(t+h) = x(t) + a*h*f(x(t),t) + b*h*f( x(t)+c*f(x(t),t)*h, t+d*h ) or x(t+h) = x + a*h*f + b*h*f( x +c*f*h , t+d*h ) where a,b,c,d are constants to be determined to make the LTE = O(h^3) Do a Taylor expansion to get x(t+h) = x + a*h*f + b*h*(f + f_x*(c*f*h) + f_t*(d*h)) + O(h^3) = x + h*(a+b)*f + h^2*(f_x*f*(b*c) + f_t*(b*d)) + O(h^3) Comparing to the above expressions (*) if we pick a,b,c,d so that a+b = 1 (**) b*c = 1/2 b*d = 1/2 Then we will have x(t+h) = x + h*x' + h^2/2*x'' + O(h^3) = true value + O(h^3) and the LTE = O(h^3) It remains to solve (**). Two solutions are (1) a=0, b=1, c=d=1/2, also called the "modified Euler Method": F1 = h*f( x(t), t ) F2 = h*f( x(t) + F1/2, t + h/2 ) x(t+h) = x(t) + F2 (2) a=b=1/2, c=d=1, also called "Heun's Method" G1 = h*f( x(t), t ) G2 = h*f( x+G1, t+h ) x(t+h) = x(t) + (G1+G2)/2 Ex: Same ODE as above with Modified Euler: Compute answer with step h=.1, .01, .001 [T1,Y1]=RK2_ModEuler(@chase,[0,10],4,.1); [T01,Y01]=RK2_ModEuler(@chase,[0,10],4,.01); [T001,Y001]=RK2_ModEuler(@chase,[0,10],4,.001); plot(T1,Y1,'c'), plot(T01,Y01,'m'), plot(T001,Y001,'g'), Look at final error: yfinal = y(length(y)); Y1final = Y1(length(Y1)); Y01final = Y01(length(Y01)); Y001final = Y001(length(Y001)); [yfinal-Y1final; yfinal-Y01final; yfinal-Y001final] Result: 6.8e-04 6.4e-06 6.3e-08 Ex: Same ODE as above with Heun: Compute answer with step h=.1, .01, .001 [T1,Y1]=RK2_Heun(@chase,[0,10],4,.1); [T01,Y01]=RK2_Heun(@chase,[0,10],4,.01); [T001,Y001]=RK2_Heun(@chase,[0,10],4,.001); plot(T1,Y1,'c'), plot(T01,Y01,'m'), plot(T001,Y001,'g'), Look at final error: yfinal = y(length(y)); Y1final = Y1(length(Y1)); Y01final = Y01(length(Y01)); Y001final = Y001(length(Y001)); [yfinal-Y1final; yfinal-Y01final; yfinal-Y001final] Result: 8.7e-04 8.2e-06 8.2e-08 4th order Runge-Kutta Methods (LTE = O(h^5)) are among the most popular. We omit the detail derivation, but just state the method: H1 = h*f( x(t), t ) H2 = h*f( x(t)+H1/2, t+h/2 ) H3 = h*f( x(t)+H2/2, t+h/2 ) H4 = h*f( x(t)+H3, t+h ) x(t+h) = x(t) + (H1 + 2*H2 + 2*H3 + H4)/6 Ex: Same ODE as above with 4th order RK Compute answer with step h=.1, .01, .001 [T1,Y1]=RK4(@chase,[0,10],4,.1); [T01,Y01]=RK4(@chase,[0,10],4,.01); [T001,Y001]=RK4(@chase,[0,10],4,.001); plot(T1,Y1,'c'), plot(T01,Y01,'m'), plot(T001,Y001,'g'), Look at final error: yfinal = y(length(y)); Y1final = Y1(length(Y1)); Y01final = Y01(length(Y01)); Y001final = Y001(length(Y001)); [yfinal-Y1final; yfinal-Y01final; yfinal-Y001final] Result: 3.0e-7 2.8e-11 -1.3e-13 Why is this not 1e4 times smaller? Resolve more accurately: options = odeset('RelTol',2.22e-14,'AbsTol',2.22e-14); hold off, clf, [t,y]=ode45(@chase,[0 10],4,options); plot(t,y), hold on New results: 3.0e-7 2.8e-11 -4.7e-16 as expected Choosing h adaptively to control error and cost ode45 does not use the same h at every step; sometimes it is larger and sometimes smaller: hold off, plot(t,y), hold on, n = length(t); hh = t(2:n) - t(1:n-1); plot(t(1:n-1),1000*hh,'r'), grid ode45 does this in order to control the error as specified by the user (in options). Making h smaller makes the method more expensive, because more steps are needed to get to the same t_final. Making h smaller also makes the method more accurate. So how do we choose h? So ideally, given a user-specified error bound, we would choose h just small enough so that the error is less than the user-specified error bound, but no smaller. In other words, do the least possible work needed to achieve the error bound. We need to talk about 1) How to estimate the LTE at each step 2) How to choose h to keep LTE at a desired level 3) What the global error bound is (next section) We will mimic the algorithm we used for "adaptive quadrature", where the algorithm estimates its own error, and chooses h just small enough to keep the LTE less than the desired level. We need some notation. We will let X(t,t0,x0) denote the true solution of X'(t) = f(X(t),t) starting at X(t0)=x0. 1) There are two methods to estimate the LTE. Here is the first method to compute the LTE. Suppose we have two methods to compute x(t+h), one with LTE C*h^k, and another more accurate method with LTE C1*h^(k+1). Write x(t+h) = X(t+h,t,x(t)) + C*h^k + O(h^(k+1)) x1(t+h) = X(t+h,t,x(t)) + C1*h^(k+1) + O(h^(k+2)) Subtracting yields x(t+h) - x1(t+h) = C*h^k + O(h^(k+1)) = LTE When x(t+h) is given by the 4th order Runge-Kutta method, and x1(t+h) is given by a 5th order Runge-Kutta method, the method is called the Runge-Kutta-Fehlberg method, after its inventors. Here is the second method to compute the LTE. Suppose we are given a method where the LTE is O(h^k). In other words x(t+h) = true value of ODE solution at t+h starting at x(t) + C*h^k + O(h^(k+1)) = X(t+h,t,x(t)) + C*h^k + O(h^(k+1)) Now suppose we compute the solution at t+h another way, by taking two steps of the same method, each with step h/2, yielding x1(t+h). We get two local truncation errors, each of size C*(h/2)^k + O(h^(k+1)) so that x1(t+h) = X(t+h,t,x(t)) + C*h^k*2^(1-k) + O(h^(k+1)) Subtracting yields x(t+h)-x1(t+h) = C*h^k*(1-2^(1-k)) + O(h^(k+1)) or (x(t+h)-x1(t+h)) / (1-2^(1-k)) = C*h^k + O(h^(k+1)) = LTE 2) Now that we can estimate the LTE, what do we do with it? How do we choose h to make it just as small as desired? Here is the simplest way to do it: t = t0 h_min = smallest value of h to use, to keep from taking too tiny steps h_max = largest value of h to use, to get enough solution points. LTE_bnd = upper bound on LTE desired by user k = order of LTE (i.e. expect LTE = O(h^k)) repeat tt = t + h compute xx = value of solution at t+h using step h compute LTE_estimate using formulas above if LTE_estimate > LTE_bnd ... LTE too big, so make h smaller h = h/2 else if LTE_estimate < LTE_bnd/2^k ... LTE too small, so make h bigger, but not too big h = min( 2*h, h_max ) else ... h ok as is, take a step t = t + h x(t+h) = xx endif until t >= t_final or h < h_min ... quit if can't make h small enough to make LTE_estimate <= LTE_bnd Alternatively, we might want to keep the sum of the LTE less than a bound LTE_sum. By making sure that LTE_estimate <= LTE_sum * h / (t_final - t0) instead of LTE_estimate <= LTE_bnd we can be sure that the sum of all the LTE_estimates is at most LTE_sum / (t_final - t0) * sum of all the h's = LTE_sum / (t_final - t0) * (t_final - t0) = LTE_sum as desired. We will use this below. Proof of convergence Now we prove that if we pick h (or each h(i)) small enough, then the final error will be as small as desired. All of the methods we have seen so far are of the form (*) x(t+h) = x(t) + h*something and have a LTE proportional to h^k for k=2,3,4,... Methods like (*) are called "one step methods" because x(t+h) only depends on x(t), not on earlier values like x(t-h); we will study such methods later. Now we will prove that, under some reasonable assumptions, that these one-step methods actually converge, in other words that the difference between the computed x(t_final) and true x(t_final) goes to 0 as h goes to 0. We will prove more: that if the LTE is O(h^k), then the error at t_final is O(h^(k-1)). We need some notation. We will let X(t,t0,x0) denote the true solution of X'(t) = f(X(t),t) starting at X(t0)=x0. We will let x(t) denote the computed solution at t, as given by (*). We want to show that the final error | X(t_final,t0,x0) - x(t_final) | = O(h^(k-1)) To see how the error made at one step (the LTE) can contribute to the final error, let's imagine that we start with x(t0) = x0 and compute x(t0+h) from (*), and then make no more errors. In other words, suppose that x(t_final) = X(t_final,t0+h,x(t0+h)), the exact solution at t_final but starting at X(t0+h) = x(t0+h), i.e. a slightly wrong value at t0+h. Thus we want to bound (**) | X(t_final,t0,x0) - X(t_final,t0+h,x(t0+h)) | = | X(t_final,t0+h,X(t0+h,t0,x0)) - X(t_final,t0+h,x(t0+h)) | or the difference between the true solutions of the ODE at t_final, but starting at slightly different initial conditions at t0+h: the true solution X(t0+h,t0,x0) and the computed solution x(t0+h). So what we need to do is bound the difference between two true solutions of the ODE starting at different initial conditions. We need to make some assumption about how "nice" f(x,t) is, and will make the simplest one possible: Lemma: Let x'(t) = f(x(t),t) be an ODE, and denote the solution that starts at x(t0)=x0 by X(t,t0,x0). Suppose that |f_x(x,t)| <= M for some constant M (recall f_x = partial f / partial x). (This guarantees that f is Lipschitz continuous in x, so that a unique solution exists by Theorem 2 above.) Then | X(t,t0,x0) - X(t,t0,x1) | <= exp(M*(t-t0)) * |x0-x1| In other words, the difference between the two solutions starting at x0 and x1 is proportional to |x0-x1| and grows at most exponentially in t. Proof: | X'(t,t0,x0) - X'(t,t0,x1) | = | f(X(t,t0,x0),t) - f(X(t,t0,x1),t) | by the definition of the ODE = | f_x(z(t),t) * (X(t,t0,x0) - X(t,t0,x1)) | by Rolle's Theorem <= M * | X(t,t0,x0) - X(t,t0,x1) | since |f_x| <= M Now let D(t) = log | X(t,t0,x0) - X(t,t0,x1) | so D'(t) = +- |X'(t,t0,x0)-X'(t,t0,x1)| / |X(t,t0,x0)-X(t,t0,x1)| <= M and integrating we get D(t) <= M*t + constant. To figure out the constant, we know we start at D(t0) so and so D(t) <= M*t + D(t0) - M*t0 = M(t-t0) + D(t0) and exponentiating we get exp(D(t)) = | X(t,t0,x0) - X(t,t0,x1) | <= exp(M*(t-t0)) * exp(D(t0)) = exp(M*(t-t0)) * | X(t0,t0,x0) - X(t0,t0,x1) | = exp(M*(t-t0)) * | x0 - x1 | as desired. Going back to (**) we get | X(t_final,t0+h,X(t0+h,t0,x0)) - X(t_final,t0+h,x(t0+h)) | <= exp(M*(t_final-t0-h)) * |X(t0+h,t0,x0) - x(t0+h) | <= exp(M*(t_final-t0) * |LTE| So if we only made the first LTE and no subsequent errors, the total error would be exp(M*(t_final-t0) * |LTE|. But of course we make a LTE at every step, bounded the same way, so we need to add these errors up. (draw picture) To do this, we need some more notation. Let h(1), h(2),... be the steps, and t(i) = t(0) + sum_{k=1}^i h(k), with t(N) = t_final. Then x(t(i)) is the i-th computed solution, and x(t(N)) is the final computed solution. We want to bound x(t(N)) - X(t(N),t(0),x(0)), where t(0) = t0 and x(0) = x0. We use the following telescoping series: x(t(N)) - X(t(N),t(0),x(0)) = x(t(N)) - X(t(N),t(N-1),x(t(N-1)) +X(t(N),t(N-1),x(t(N-1)) - X(t(N),t(N-2),x(t(N-2))) +X(t(N),t(N-2),x(t(N-2)) - X(t(N),t(N-3),x(t(N-3))) + ... +X(t(N),t(1),x(t(1)) - X(t(N),t(0),x(t(0))) A term in this series looks like X(t(N),t(i),x(t(i)) - X(t(N),t(i-1),x(t(i-1))) = X(t(N),t(i),x(t(i)) - X(t(N),t(i),X(t(i),t(i-1),x(t(i-1)))) (picture of what this means) Our Lemma bounds the absolute value of this by exp(M*(t(N)-t(i)))*|x(t(i)) - X(t(i),t(i-1),x(t(i-1)))| <= exp(M*(t(N)-t(0)))*|x(t(i)) - X(t(i),t(i-1),x(t(i-1)))| = exp(M*(t(N)-t(0)))*LTE(i) where LTE(i) is the local truncation error at the i-th step. Applying the triangle inequality to our telescoping series and using the above bound for each term yields: | x(t(N)) - X(t(N),t(0),x(0)) | <= exp(M*(t_final-t0)) * sum_{i=1}^N LTE(i) In the simplest case, with all steps h(i)=h, LTE(i) = O(h^k), | x(t(N)) - X(t(N),t(0),x(0)) | <= exp(M*(t_final-t0)) * (t(N)-t(0))/h * O(h^k) = exp(M*(t_final-t0)) * O(h^(k-1)) as desired. Alternatively, if the user asked the algorithm to choose the step size h(i) to keep sum_{i=1}^N LTE(i) less than a bound LTE_sum as discussed above, then the error satisfies | x(t(N)) - X(t(N),t(0),x(0)) | <= exp(M*(t_final-t0)) * LTE_sum In practice one probably doesn't know M, and if one did, this bound may still be too large. But it is enough to show that we converge, i.e. that the final computed answer x(t_final) approaches the true solution X(t_final,t0,x0) as the sum of all the LTEs goes to zero. Multistep Methods All the methods we have discussed so far are called "one-step methods", because they all have update formulas that look like x(t+h) = x(t) + h * something that depends only on t and x(t) In particular they do use any values of x(s) for s<t, say x(t-h), x(t-2*h), ... It turns out that if we use this available information about x(s), then we can be much faster than Runge-Kutta methods, at least when f(x,t) is smooth enough. Recall that to make the LTE smaller in a Runge-Kutta method, we needed to evaluate f(x,t) at more and more values of (x,t). Here is table of the number of function evaluations needed to make the LTE = O(h^(k+1)) # function evaluations: 1 2 3 4 5 6 7 8 ... order k in LTE = O(h^(k+1)) 1 2 3 4 4 5 6 6 ... (Note that if the LTE = O(h^(k+1)), then the global error is O(h^k)). In other words, smaller LTE => higher k => more function evaluations => more work at each step. If evaluating f(x,t) is the most expensive part of the algorithm (and it usually is), then we see that to get LTE = O(h^(k+1)) we need at least k function evaluations, and more when k is 5 or more. Multistep methods will let us make k larger with just 1 function evalution per step. So they can be much faster than Runge-Kutta methods. However, they only converge when f(x,t) is smooth enough everywhere (at least k derivatives). In contrast, Runge-Kutta methods, properly designed, can handle corners and other non-smooth f(x,t) with less difficulty. Also, their convergence proofs are different than Runge-Kutta methods, and it is easy to write down reasonable looking multistep methods that actually do not work, even though their LTEs are very small. Adams-Bashforth Methods To derive these, we write x(t+h) = x(t) + (x(t+h)-x(t)) = x(t) + integral from t to t+h x'(s) ds = x(t) + integral from t to t+h f(x(s),s) ds so the problem is to approximate the above integral from t to t+h. The information we have available is x(t), x(t-h), x(t-2*h),...,x(t-m*h) along with f0 = f(x(t),t), f1 = f(x(t-h),t-h), f2 = f(x(t-2*h),t-2*h), ... , fm = f(x(t-m*h),t-m*h) if we save m old function values. As we have done before, the idea is to 1) Approximate the function that passes through (t,f0), (t-h,f1), (t-2*h,f2), ... , (t-m*h,fm) by a polynomial p(s) 2) Compute the integral from t to t+h p(s) ds = sum_{i=0}^m C_i*fi This is very similar to the way we approximated integrals before (integrate a polynomial approximation p(s)) but differs because the the interval [t,t+h] of integration does not contain the points of interpolation t, t-h, t-2*h, ..., t-m*h. We can still use the method of undetermined coefficients to figure out the coefficients C_i, by asking that the value of the integral is exact for polynomials p(s) of degree at most m. For example, when m=1 we choose C_0 and C_1 so that the formula is exact when p(s) = 1: integral from t to t+h 1*ds = h = C_0 + C_1 p(s) = s: integral from t to t+h s*ds = h^2/2 + t*h = C_0*t + C_1*(t-h) Solving these two linear equations in two unknowns yields C_0 = h*3/2, C_1 = -h/2 Continuing this way, the first few methods are m=0 => x(t+h) = x(t) + h*f0 (Euler's Method) m=1 => x(t+h) = x(t) + h*[1.5*f0 - .5*f1] m=2 => (homework) m=3 => x(t+h) = x(t) + (h/24)*[55*f0 - 59*f1 + 37*f2 - 9*f3] m=4 => x(t+h) = x(t) + (h/720)*[1901*f0 - 2774*f1 + 2616*f2 - 1274*f3 + 251*f4] The way we compute the LTE and then the global error is very similar to before, and we just outline it here. 1) From our theory about the error in polynomial interpolation, we get that x'(s) = f(x(s),s) = p(s) + Error = p(s) + x^(m+2)(z)/(m+2)! (s-t)(s-t+h)(s-t+2h)...(s-t+mh) = p(s) + O(h^(m+1)) where we assume that t <= s <= t+h 2) Integrating, we get that x(t+h)-x(t) = integral from t to t+h x'(s) ds = integral from t to t+h p(s) ds + O(h^(m+2)) or x(t+h) = x(t) + integral from t to t_h p(s) ds + LTE where LTE = O(h^(m+2)) 3) From our proof of convergence, we conclude that the global error is O(h^(m+1)) Adams-Moulton Methods This is very similar to Adams-Bashforth, except that the data we fit with a polynomial for integration includes fm1 in the list: fm1 = f(x(t+h),t+h) f0 = f(x(t),t), f1 = f(x(t-h),t-h), f2 = f(x(t-2*h),t-2*h), ... , fm = f(x(t-m*h),t-m*h) Now the procedure looks similar as before: (1) find a polynomial p(s) that interpolates (t+h,fm1), (t,f0), (t-h,f1), (t-2*h,f2), ... , (t-m*h,fm) (2) Compute the integral from t to t+h p(s) ds = sum_{i=-1}^m C_i*fi (3) Update x(t+h) = x(t) + sum_{i=-1}^m C_i*fi Since we are using one more interpolation point than before (m+2, not m+1) the interpolation error is smaller (O(h^(m+2)), not O(h^(m+1)), the LTE is smaller (O(h^(m+3)), not O(h^(m+2)) and the global error is smaller (O(h^(m+2)), not O(h^(m+1)) For example, when m=0, we get the following method with O(h^2) global error: x(t+h) = x(t) + h/2*(f0 + fm1) = x(t) + h/2*(f(x(t),t) + f(x(t+h),t+h)) But how would be actually implement this, since both sides depend on x(t+h)? It look like we need to solve a (nonlinear) equation for x(t+h): x(t+h) - h/2*f(x(t+h),t+h) = x(t) + h/2*f(x(t),t) Since Adams-Moulton requires us to solve a (nonlinear) equation for x(t+h), it is called an implicit method, whereas Adams-Bashforth is called an explicit method, since it does not. Using an implicit method seems so much harder than an explicit method that it seems like a really bad idea. But in fact it is quite useful, because we get a really good approximation for x(t+h) from an Adams-Bashforth method: Predictor-Corrector Methods Here is how we combine Adams-Bashforth and Adams-Moulton to get something really useful. We let Adams-Bashforth(m) and Adams-Moulton(m) denote the methods using values back to t-m*h, as above. at each time step Compute x(t+h) using Adams-Bashforth(m), so its LTE is O(h^(m+2)) Compute x1(t+h) using Adams-Moulton(m), using the x(t+h) just computed by Adams-Bashforth(m) to evaluate fm1 = f(x(t+h),t+h). Thus there is no nonlinear equation to solve, and the LTE is still O(h^(m+3)) Use x(t+h) - x1(t+h) as an approximate LTE for controlling h This method is called "predictor-corrector" because Adams-Bashforth is used to predict x(t+h) and Adams-Moulton is used to correct it. This has all the ingredients we want: (1) explicit method (no nonlinear equation solving) (2) one function evaluation f(x(t),t) per time step (3) simple formula to estimate LTE in order to control step size h One last question: how do we get started from x(t0)=x0, if the formulas need m old values x(t), x(t-h), ..., x(t-m*h)? Answer: we use a Runge-Kutta method with the same LTE to generate x(t0+h), x(t0+2*h),..., x(t0+m*h) and then switch to Predictor-Corrector. Convergence of Multistep Methods When the method is of the form of Adams-Bashforth or Adams-Moulton: x(t+h) = x(t) + h*something then our earlier convergence proof does not change. But one can imagine other methods, all of which have the following form: We use the notation t_i = t-i*h, x_i = x(t_i) and f_i = f(x_i,t_i), so that x_(-1) = x(t+h) is the next value to compute: a_(-1)*x_(-1) + a_0*x_0 + a_1*x_1 + ... + a_m*x_m = h*[ b_(-1)*f_(-1) + b_0*f_0 + b_1*f_1 + ... + b_m*f_m ] Here the a_i and b_i are constants that determine the method. We can assume a_(-1) = 1 since we can divide the formula through by a_(-1) If b_(-1) = 0 then the method is explicit, else it is implicit. The general question is what values of a_i and b_i yield methods that converge as h -> 0? It is intuitively clear that the LTE must go to zero as h -> 0 (all our analysis has assumed that taking a smaller h makes the error at each step smaller), but more is required. To see why, consider the method x_(-1) + 4*x_0 - 5*x_1 = h*(4*f_0 + 2*f_1) or x(t+h) = -4*x(t) + 5*x(t-h) + h*(4*f(x(t),t) + 2*f(x(t-h),t-h)) This seems like a reasonable method, because LTE = -h^4/8*x''''(z) = O(h^4). But let's try running it on x' = 1, x(0)=0, whose solution should just be x(t)=t. To get started we need two values, so we use the exact solution values x(0)=0, x(-h) = -h: hold off, [t,x]=ode_boom(@unit,[0 1],[0 -.01],.01); plot(t,x),semilogy(t,x), plot(x(2:end)./x(1:end-1)) So the ratio between consecutive computed solution x(t+h)/x(t) approaches the constant -5. Here is a small table of results (the exponential growth is apparent, a far cry from the true solution x(t)=t: t x(t) 5.0000e-02 -1.0370e+01 6.0000e-02 5.2140e+01 1.0000e-01 3.2552e+04 1.1000e-01 -1.6276e+05 1.5000e-01 -1.0173e+08 1.6000e-01 5.0863e+08 Let's try a smaller h (.005 instead of .01), maybe that will help: hold off, [t,x]=ode_boom(@unit,[0 1],[0 -.005],.005); plot(t,x),semilogy(t,x), plot(x(2:end)./x(1:end-1)) Same result (x(t+h)/x(t) -> -5). Let's try a different ode: hold off, [t,x]=ode_boom(@chase,[0 1],[0 -.01],.01); plot(t,x),semilogy(t,x), plot(x(2:end)./x(1:end-1)) Same solution (x(t+h)/x(t) -> -5) People do use these multistep methods, but obviously a small LTE is not enough to guarantee convergence, as was the case with one-step methods. Let's analyze the above method applied to x'(t) = 0, x(0) = 1, so the true solution is x(t)=1. Thus h*(4*f_0 + 2*f_1) = 0 exactly. This means we want to solve the "difference equation" x(t+h) + 4*x(t) - 5*x(t-h) = 0 yielding the sequence x(t+2*h), x(2+3*h),... and so on. Let's get the most general solution of this equation. To simplify notation, let start at t0=0. To get started, we need two initial values, say x(0) and x(h), from which we compute x(2*h), x(3*h) and so on. Let's "guess" that the solution x(k*h) = r^k for some value of r that we need to solve for. Plug into (*) x((k+1)*h) + 4*x(k*h) - 5*x((k-1)*h) = 0 to get r^(k+1) + 4*r^k - 5*r^(k-1) = 0 or, dividing through by r^(k-1) r^2 + 4*r - 5 = 0 We can solve this, yielding r = 1 or r = -5. Thus x(k*h) = (-5)^k and x(k*h) = (1)^k = 1 both solve (*), and since (*) is linear, so does x(k*h) = a*(-5)^k + b for any constants a and b. But we have not used our initial conditions x(0) and x(h) yet: we will use a and b to make sure our solution matches them: x(0) = a*(-5)^0 + b = a+b x(h) = a*(-5)^1 + b = -5*a +b Solving for a and b yields a = (x(0) - x(h))/6 b = 5*x(0)/6 + x(h)/6 or x(k*h) = (x(h)-x(0))/6*(-5)^k + (5*x(0)/6+x(h)/6) We did all this to see if the method converges as h -> 0. What this means more precisely is that as long as the LTE is small enough, then the error goes to zero for any initial conditions. We are given x(0)=1, so the question is whether the error -> 0 as h->0 for any small enough LTE x(h)-1 = x(h)-x(0). Let's evaluate the error at t=1, that is k = 1/h, yielding x(1) = (x(h)-1)/6*(-5)^(1/h) + (5/6 + x(h)/6) No matter how small the LTE is, O(h^m) for any m, we get that x(1) = O(h^m*(-5)^(1/h)) -> infinity as h -> 0, which does not converge! Furthermore, it is easy to see that x(t+h)/x(t) -> -5 as h -> 0, as we saw in our numerical examples. So we see that some possible multistep methods cannot possibly converge. To see what the general case looks like, we need to talk about solving difference equations like (*) above. Solving Difference Equations Read section 1.3 of the text. A general difference equation is of the form (*) a_k*y_{k+n} + a_{k-1}*y_{k-1+n} + a_{k-2}*y_{k-2+n} + ... + a_0*y_n = 0 where a_k is nonzero. We are given values of y_0, y_1, ... , y_{k-1} and use (*) to compute y_k, y_{k+1},... We want to compute an explicit solution of (*). Here is the general approach: We "guess" that the solution y_n = r^n for some constant r that we don't know. Plug in to get a_k*r^(k+n) + a_{k-1}*r^(k-1+n) + ... + a_0*r^n = 0 divide through by r^n to get (**) a_k*r^k + a_{k-1}*r^(k-1) + ... + a_0 = 0 a polynomial of degree k. Let's consider the simplest case, when all k solutions r_1,...,r_k of the polynomial (**) are distinct. This means that y_n = r_i^n is a solution of (*) for any i, as is any linear combination (***) y_n = a_1*r_1^n + a_1*r_2^n + ... + a*k*r_k^n for any constants a_1,...,a_k. We use these k free paramters to match the initial condition y_0,..., y_{k-1}: (****) y_0 = a_1*r_1^0 + a_1*r_2^0 + ... + a*k*r_k^0 y_1 = a_1*r_1^1 + a_1*r_2^1 + ... + a*k*r_k^1 ... y_{k-1} = a_1*r_1^(k-1) + a_1*r_2^(k-1) + ... + a*k*r_k^(k-1) This is a system of k linear equations for the a_1,...,a_k. The coefficient matrix is the Vandermonde that we have seen before, which is nonsingular. We see that if any |r_i|>1, then the solution of (*) will go to infinity like r_i^n unless the initial conditions are very special so that a_i = 0. This motivates us to make the following definition: Def: When all its roots are distinct, the difference equation is Stable if all |r_i| <= 1. Otherwise, if at least one |r_i|>1, we call it unstable (yet another use of the word stable...). When the polynomial (**) has multiple roots, the idea is similar, but a little more complicated. If r_i is a multiple root of (**), occuring t times (i.e. (r-r_i)^t divides (**) evenly), then we get the following different solutions of (*): y_n = r_i^n as before y_n = n*r_i^n y_n = n^2*r_i^n ... y_n = n^(t-1)*r_i^n We can still write down the general solution of (*) like (***), and solve for the unknown coefficients a_i as in (****). (We get a "confluent Vandermonde" matrix to solve that is still nonsingular.) Stability still depends on the |r_i| being at most 1 (and if |r_i|=1, then it must not be multiple). Convergence of Multistep Methods, again Now we can return to the convergence of the multistep method (*) a_(-1)*x_(-1) + a_0*x_0 + a_1*x_1 + ... + a_m*x_m = h*[ b_(-1)*f_(-1) + b_0*f_0 + b_1*f_1 + ... + b_m*f_m ] We know enough to state and understand (but not prove) the general Theorem about when such a method converges as h -> 0. Def: (*) is "consistent" if the LTE is O(h^k) for k >= 2. Def: (*) is "stable" if the difference equation a_(-1)*x_(-1) + a_0*x_0 + a_1*x_1 + ... + a_m*x_m = 0 is stable. Theorem: (*) converges to the answer of the ODE if and only if it is consistent and stable. We may apply this to Adams-Bashforth and Adams-Moulton methods, where 1*x_{-1} - 1*x_0 = h*something. In this case the polynomial is r-1, which has a single root at 1, and so is stable. This theorem is important because people regularly use multistep methods besides Adams-Bashforth and Adams-Moulton, which have the simple polynomials r-1. We'll see some in the next section. Stiff ODE Solvers Read section 8.12 of the text. To motivate this material, let's look at a sequence of numerical examples. We will return to our old example x'(t) = c*(sin(t) - x(t)), c > 0, x(0) = 4, The qualitative behaviour of this is for x(t) to "chase" sin(t), For when c is positive, the ODE tells us that x(t) decreases if x(t) > sin(t) and x(t) increases if x(t) < sin(t). x(t) "chases" sin(t) more closely the larger c is. So for large c, x(t) will quickly drop from its initial value x(0) = 4, and then very closely follow the smooth curve sin(t). We want to figure out how our numerical methods that we have studied so far behave as c increases, and see if we need some new ones (we will!). % set options to control accuracy of ode45 to get "right answer" opts = odeset('AbsTol',1e-12,'RelTol',1e-12); % last argument to ode45 is value of c=10 in c*(sin(t)-x(t)) c = 10; [t,y]=ode45(@chasec,[0 10],4,opts,c); % plot "right answer", along with sin(t) curve being "chased" plot(t,y,'b',t,sin(t),'r') Now let's try solving using our old favorite Euler's method x(t+h) = x(t) + h*f(x(t),t), for h=.1 and c=10, and plotting with the true answer h = .1; [T,Y]=eulerc(@chasec,[0 10],4,h,c); plot(t,y,'b',t,sin(t),'r',T,Y,'k') The blue (ode45) and black (euler) curves are very close past the initial steep region, and indeed the euler error at t=10 is 2e-3. Now let's start raising c to see what happens: c = 15 [t,y]=ode45(@chasec,[0 10],4,opts,c); [T,Y]=eulerc(@chasec,[0 10],4,h,c); plot(t,y,'b',t,sin(t),'r',T,Y,'k') Euler looks "unstable" for small t, but the error at t=10 is still 1.5e-3 c = 18 [t,y]=ode45(@chasec,[0 10],4,opts,c); [T,Y]=eulerc(@chasec,[0 10],4,h,c); plot(t,y,'b',t,sin(t),'r',T,Y,'k') The "ringing" looks worse, but the error at t=10 is still 1.3e-3 c = 19 [t,y]=ode45(@chasec,[0 10],4,opts,c); [T,Y]=eulerc(@chasec,[0 10],4,h,c); plot(t,y,'b',t,sin(t),'r',T,Y,'k') Now the error is still just 1.2e-2 at t=10, even though it is much larger for smaller t. If we keep doing this, as long as c < 20, the solution will eventually converge to the true solution, but take longer and longer: c = 19.7 [t,y]=ode45(@chasec,[0 10],4,opts,c); [T,Y]=eulerc(@chasec,[0 10],4,h,c); plot(t,y,'b',t,sin(t),'r',T,Y,'k') If we take c > 20, the solution blows up entirely, with growing oscillations: c = 20.3 [t,y]=ode45(@chasec,[0 10],4,opts,c); [T,Y]=eulerc(@chasec,[0 10],4,h,c); plot(t,y,'b',t,sin(t),'r',T,Y,'k') And finally, at c=20, the solution just rings constantly, without the amplitude shrinking or growing: c = 20 [t,y]=ode45(@chasec,[0 10],4,opts,c); [T,Y]=eulerc(@chasec,[0 10],4,h,c); plot(t,y,'b',t,sin(t),'r',T,Y,'k') The solution of the ODE seems quite nice, close to sin(t), but in fact it has an interesting property called "stiffness" which makes it hard to solve. Roughly speaking, the ODE is "stiff" if the solution looks "smooth," so that not many points (not a tiny h) are needed to represent it accurately, but where Euler requires many more points (tiny h) to converge. Now we will do a simpler example that behaves exactly the same way, but where we can do all the analysis by hand to see what is going on: x'(t) = -c*x(t), x(0)=1. The solution is x(t) = exp(-c*t), which "chases the curve x=0": c = 10; [t,y]=ode45(@const,[0 10],1,opts,c); [T,Y]=eulerc(@const,[0 10],1,h,c); plot(t,y,'b',T,Y,'k') c = 19; [t,y]=ode45(@const,[0 10],1,opts,c); [T,Y]=eulerc(@const,[0 10],1,h,c); plot(t,y,'b',T,Y,'k') c = 20; [t,y]=ode45(@const,[0 10],1,opts,c); [T,Y]=eulerc(@const,[0 10],1,h,c); plot(t,y,'b',T,Y,'k') c = 21; [t,y]=ode45(@const,[0 10],1,opts,c); [T,Y]=eulerc(@const,[0 10],1,h,c); plot(t,y,'b',T,Y,'k') Now let's do this one by hand, to see exactly what is happening: Euler's method says: x(t+h) = x(t) + h*f(x(t),t) = x(t) + h*(-c*x(t)) = (1-h*c)*x(t) = (1-h*c)^2 * x(t-h) = ... Letting t_i = i*h yields x(t_i) = x(i*h) = (1-h*c)^i*x(0) Under what conditions will this decay to 0 as desired? Clear if and only if | 1-h*c | < 1, or, when c is a positive real number, 0 < h < 2/c For our example, with h=.1, this means we only expect Euler to converge to zero when .1 < 2/c or c < 20, which is exactly what we saw, in both examples. When c=20, we get x(t+h) = -x(t), which is precisely the "ringing" we saw: the oscillations will never blow up or decay. Conversely, given c, we have to pick h < 2/c just to get an answer that doesn't blow up, let alone is accurate. As c gets larger (which only makes the answer "flatter", we keep having to take tinier and tinier h, and do more and more work. For example, here is a table of how many steps ode45 had to take to get the answer for x'(t) = c*(sin(t)-x(t)) to within 1e-12, as a function of c c # steps ---- ----- 1 2389 10 11053 15 14037 20 16621 40 24917 80 37385 100 42617 200 64161 The # steps is growing with c (roungly like c^.6), so if we use our current algorithms, we will have to do more work with larger c. We want methods that will let us take a much larger step h without going unstable in the ways we have seen. In other words efficiency becomes the most important criterion. We will be able to do this, but the price will be using "implicit methods" where we have to solve a nonlinear equation (or system of equations) at every time step. These methods will be called "stiff solvers" because they are designed for such ODEs. We recall how we derived Euler, and modify it to get the simplest stiff solver, which we will call Backward Euler: exact x(t+h) = x(t) + h*x'(t) + h^2/2*x''(t) + ... Forward Euler x(t+h) = x(t) + h*x'(t) = x(t) + h*f(x(t),t) with LTE = O(h^2) from truncating the Taylor expansion exact x(t) = x(t+h) + (-h)*x'(t+h) + (-h)^2/2*x''(t+h) + ... Backward Euler x(t) = x(t+h) - h*x'(t+h) = x(t+h) - h*f(x(t+h),t+h) with LTE = O(h^2) from truncating the Taylor expansion or x(t+h) = x(t) + h*f(x(t+h),t+h) Taking one step requires solving this (possibly nonlinear) equation for x(t+h), so it is "implicit". Let's apply this to x'(t) = -c*x(t), x(0)=1, since we can solve for x(t+h) by hand:: x(t+h) = x(t) + h*f(x(t+h),t+h) = x(t) + h*(-c*x(t+h)) or x(t+h) = x(t)/(1+c*h) or x(t_i) = x(i*h) = x(0)/(1+c*h)^i Now x(t_i) -> 0 for any positive c and h, just like the true solution x(t_i) = exp(-c*t_i) = exp(-c*h*i). This attractive property is named in the following definition: Def: An algorithm for solving x'(t) = -c*x(t) is A-stable (A stands for "Absolutely") if the computed solution converges to 0 when the true solution converges to 0. And when does the true solution converge to 0? When c is complex too? Trying this on x'(t) = 30*(sin(t)-x(t)), x(0)=4, h=.1, yields a stable solution with error -7.6e-4 at t=10: c=30;,h=.1; [t,y]=ode45(@chasec,[0 10],4,opts,c); [T,Y]=back_eulerc([0 10],4,.1,c); plot(t,y,'b',T,Y,'k',t,sin(t),'r') y(end)-Y(end) % = -7.6e-4 Just as with prior methods (Runge-Kutta, Adams-Backforth-Moulton) we want a smaller truncation error than O(h^2), preferably O(h^m) for larger m. Here is how we derive one. We use ideas from an earlier part of the course on numerical differentiation. These methods are called BDF (Backward Difference Formula) Methods. The idea is to 1) approximate x'(t+h) by a linear combination of x(t+h), x(t), x(t-h),..., x(t-m*h) using numerical differentiation 2) Set this approximation equal to the true x'(t+h) = f(x(t+h),t+h) Let us recall how to do step 1, approximating x'(t+h): 1.1) Find a polynomial p(s) interpolating x(s) at (t+h,x(t+h)), (t,x(t)), (t-h,x(t-h)), ... , (t-m*h, x(t-m*h)) The error, as many times before, is O(h^#points) = O(h^(m+2)) 1.2) Use the approximation x'(t+h) ~ p'(t+h), with error O(h^(m+1)) When we finally write down the formula, we multiply through by h again, so the LTE becomes O(h^(m+2)) (see below). The simplest example is m=0, or 2 points, We get p(s) a straight line with slope (x(t+h)-x(t))/h. Proceeding with our outline we get the method (x(t+h)-x(t))/h = f(x(t+h),t+h) where Error = O(h) or, solving for x(t+h), x(t+h) = x(t) + h*f(x(t+h),t+h)) which is just Backward Euler, with LTE = O(h^(m+2)) = O(h^2) as before. The next example is m=1, or 3 points. We get p(s) a parabola whose slope at t+h is p'(t+h) = (3/(2*h))* x(t+h) - (2/h)*x(t) + (1/(2*h))*x(t-h) = x'(t+h) + O(h^2) rewrite (3/(2*h))* x(t+h) - (2/h)*x(t) + (1/(2*h))*x(t-h) = f(x(t+h),t+h) as x(t+h) - (4/3)*x(t) + (1/3)*x(t-h) = (2*h/3)*f(x(t+h),t+h) which has error O(h^3) (since we multiplied through by h) Do these methods converge as h -> 0? Backward Euler does, because it is in the form x(t+h) = x(t) + h*something for which our original convergence proof works. But the m=1 method x(t+h) - (4/3)*x(t) + (1/3)*x(t-h) = (2*h/3)*f(x(t+h),t+h) is the kind of multstep algorithm involving a linear combination of several x(t+i*h) that we analyzed in the last section. Our general theorem tell us that we need to look at the roots of the polynomial r^2 - (4/3)*r + 1/3 = 0 and ask whether they all satisfy |root| < = 1. The roots can be seen to be 1 and 1/3, both of which satisfy the condition. So this method converges. We can keep doing this up to m=5, and we keep getting stable methods that converge. For m >= 5, we get instability and no convergence. These methods are built into the Matlab function ode15s (s for "stiff"), which we try below: the ode15s is almost 15 times faster than ode45 (1429 times steps for ode15s versus 21061 for ode45): c = 30; [t,y]=ode45(@chasec,[0 10],4,opts,c); [ts,ys]=ode15s(@chasec,[0 10],4,opts,c); plot(t,y,'b',ts,ys,'k',t,sin(t),'r') [length(y);length(ys);y(end)-ys(end)] % = [21061; 1429; 1e-12] Stiff equations arise regularly in elctronic circuit simulation, chemical kinetics, molecular dynamics, and any other application where some components of the system change much more rapidly than other parts, and it would require too many time steps to use a step size h small enough to accurately compute the "fast" components. Many open problems remain in this area. ODE Software At the beginning of our section on ODEs, we listed several standard kinds of ODEs for which I said good software existed, and I strongly recommended that you use this software when you encountered an ODE that you had to solve, because the software was very carefully written, with many years of development put into it to make it fast, accurate and reliable. Just as you would not write your own C compiler or spreadsheet or internet browser (except as a classroom exercise, or to do research), you would generally not write your own ODE solver, unless your problem were very difficult and special. So here is the same list of ODEs we gave at the beginning of the lectures on ODEs, along with pointers to software that we recommend. Roughly speaking, the software comes in 2 categories: Matlab (and systems like it) which have nice user interfaces for defining inputs and outputs, especially for certain application areas using ODEs, like controls. You have already used the Matlab software (ode23, ode45, ode15s) and I will not say more about it. Libraries written in C or Fortran, which you have to call from your own C or Fortran program. These generally have more complicated user interfaces than Matlab, but solve more general and more difficult problems, and can be much faster. As before a single web site that guides you to this (and other) software available on-line is GAMS = Guide to Available Mathematical Software, or gams.nist.gov. I summarize the software it points to below. There are a lot of packages with similar functionality available, mostly on Netlib (www.netlib.org) in the subdirectories ode (www.netlib.org/ode) odepack (www.netlib.org/odepack) and sodepack (www.netlib.org/sodepack) For simplicity I will suggest one code, but there are variations on these codes in these directories that you might find useful. Initial Value Problems (IVP): x'(t) = f(x(t),t) with initial condition (IC) x(t0) = x0 f may be a vector Runge-Kutta methods are available in www.netlib.org/ode/rksuite Predictor-corrector methods are available at www.netlib.org/odepack These include algorithms for dealing with stiff problems. Boundary Value Problems (BVP) (if we have time) x''(t) = f(x'(t),x(t),t) with boundary conditions (BCs) x(t{0}) = x{0} and x(t{1}) = x{1} In this case we want x(t) for t{0} < t < t{1} www.netlib.org/ode/colnew.f Delay differential equation x'(t) = f(x(t),x(t-a),t) with x(t{0}) = x{0} or x'(t) = f(x(t),x(t-a),x(t-b),...,t) with x(t{0}) = x{0} Wille and Baker, "DELSOL - a numerical code for delay differential equations", Applied Numerical Mathematics, vol 9, 1992 Differential-Algebra equations (DAE) Suppose x has n components 0 = f(x'(t),x(t),y(t),t) (m equations) 0 = g(x(t),y(t),t) (n-m equations) www.netlib.org/ode/ddassl.f