Chapter 5Initial Value Problems5.1 Finite Difference MethodsWe don’t plan to study highly complicated nonlinear differential equations. Our firstgoal is to see why a difference method is successful (or not). The crucial questionsof stability and accuracy can be clearly understood for linear equations. Then wecan construct difference approximations to a great variety of practical problems.Another property we need is computational speed. That depends partly onthe complexity of the equation u0= f (u, t). Often we measure speed by the numberof times that f (u, t) has to be computed in each time step (that number could beone). When we turn to implicit methods and predictor-corrector methods, to improvestability, the cost per step goes up but we gain speed with a larger step ∆t.This chapter begins with basic methods (forward Euler, backward Euler) and thenimproves. Each time, we test stability on u0= a u. When a is negative, ∆t is oftenlimited by −a ∆t ≤ C. This has an immediate effect: the equation with a = −99requires a much smaller ∆t than a = −1. Let me organize the equations as scalarand vector, nonlinear in general or linear with constant coefficients a and A:1 equationN equationsu0= f(u, t)u0i= fi(u, t)u0= auu0= Aua ≈ ∂f /∂u a ≤ 0Aij≈ ∂fi/∂ujRe λ(A) ≤ 0For an ordinary differential equation u0= f(u, t), good codes will increase theaccuracy (and keep stability) far beyond the O(∆t) error in Euler’s methods. Youcan rely on freely available software like ode45 to make the two crucial decisions:1. to choose an accurate difference method (and change the formula adaptively)2. to choose a stable time step (and change ∆t adaptively).We will introduce accurate methods, and find the stability limits on ∆t.c2006 Gilbert Strangc2006 Gilbert StrangCHAPTER 5. INITIAL VALUE PROBLEMSStiff Differential EquationsFirst we call attention to one extra question: Is the equation stiff ? Let me beginwith a made-up example to introduce stiffness and its effects:v(t) = e−t+ e−99t↑ ↑controls decay controls ∆tThe step ∆t is 99 timessmaller because of e−99tthat disappears anywayThose decay rates −1 and −99 could be the eigenvalues of A, as in Example 1.Example 1ddtvw=−50 4949 −50vwwithv(0)w(0)=20. (1)The solution has v(t) = e−t+ e−99tand w(t) = e−t− e−99t. The time scales are differentby a factor of 99 (the condition number of A). The solution will decay at the slowtime scale of e−t, but computing e−99tmay require a very small ∆t for stability.It is frustrating to have ∆t controlled by the component that is decaying so fast.Any explicit method will have a requirement 99∆t ≤ C. We will see how this happensand how an implicit method (like ode15s and od23t in MATLAB) can avoid it.Trefethen [ ] points to these applications where stiffness comes with the problem:1. Chemical kinetics (reactions go at very different speeds)2. Control theory (probably the largest application area for MATLAB)3. Circuit simulation (components respond at widely different time scales)4. Method of Lines (large systems come from partial differential equations).Example 2 The N by N second difference matrix K produces a large stiff system:Method of Linesdudt=−Ku(∆x)2hasduidt=ui+1− 2ui+ ui−1(∆x)2. (2)This comes from the heat equation ∂u/∂t = ∂2u/∂x2, by discretizing only the spacederivative. Example 1 had eigenvalues −1 and −99, while Example 2 has N eigenvalues—but the difficulty is essentially the same ! The most negative eigenvalue here is abouta = −4/(∆x)2. So a small ∆x (for accuracy) will require a very small ∆t (for stability).This “semidiscrete” method of lines is an important idea. Discretizing thespace variables first produces a large system that can be given to an ODE solver. (Wehave ordinary differential equations in time.) If it wants to, that solver can vary the timestep ∆t and even the discretization formula as u(t) speeds up or slows down.5.1. FINITE DIFFERENCE METHODSc2006 Gilbert StrangThis method splits the approximation of a PDE into two parts. Finite differences/finiteelements in earlier chapters produce the first part (discrete in space). The upcomingstability-accuracy analysis applies to the second part (discrete in time). This idea is verysimple and useful, even if it misses the good methods later in this chapter that take fulladvantage of space-time. For the heat equation ut= uxx, the useful fact utt= uxxxxallows us to cancel space errors with time errors—which we won’t notice when they areseparated in the semidiscrete method of lines.Forward Euler and Backward EulerThe equation u0= f (u, t) starts from an initial value u(0). The key point is that therate of change u0is determined by the current state u at any moment t. This modelof reality, where all the history is contained in the current state u(t), is a tremendoussuccess throughout science and engineering. (It makes calculus almost as importantas linear algebra.) But for a computer, continuous time has to change to discretetime. One differential equation allows many difference equations!The simplest method (Euler is pronounced “Oiler”) uses a forward difference:Forward EulerUn+1− Un∆t= f(Un, tn) is Un+1= Un+ ∆t fn. (3)Over each ∆t interval, the slope of U doesn’t change. Figure 5.1 shows how the correctsolution to u0= au follows a smooth curve, while U(t) is only piecewise linear. Abetter method (higher accuracy) will stay much closer to the curve by using moreinformation than the one slope fn= f (Un, tn) at the start of the step.PSfrag replacementsttu(0) = 1∆t∆tu = e∆tU = 1 + ∆tu = e∆tU =11 − ∆t(forward)(backward)Figure 5.1: Forward Euler and Backward Euler for u0= u. One-step errors ≈12(∆t)2.Backward Euler comes from using fn+1at the end of the step, when t = tn+1:Backward EulerUn+1− Un∆t= f(Un+1, tn+1) is Un+1− ∆t fn+1= Un. (4)This is an implicit method. To find Un+1, we have to solve equation (4). Whenf is linear in u, we are solving a linear system at each step (and the cost is low ifc2006 Gilbert StrangCHAPTER 5. INITIAL VALUE PROBLEMSthe matrix is tridiagonal, like I − ∆t K in Example 2). We will comment later oniterations like Newton’s method or predictor-corrector in the nonlinear case.The first example to study is the linear scalar equation u0= au. Compare forwardand backward Euler, for one step and for n steps:Forward Un+1= (1 + a