Chapter 5 Initial Value Problems 5.1 Finite Difference Methods We don’t plan to study highly complicated nonlinear differential equations. Our first goal is to see why a difference method is successful (or not). The crucial questions of stability and accuracy can be clearly understood for linear equations. Then we can construct difference approximations to a great variety of practical problems. Another property we need is computational speed. That depends partly on the complexity of the equation u� = f (u, t). Often we measure speed by the number of times that f (u, t) has to be computed in each time step (that number could be one). When we turn to implicit methods and predictor-corrector methods, to improve stability, 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 then improves. Each time, we test stability on u� = a u. When a is negative, �t is often limited by −a �t ↑ C. This has an immediate effect: the equation with a = −99 requires a much smaller �t than a = −1. Let me organize the equations as scalar and vector, nonlinear in general or linear with constant coefficients a and A: N equations u � = f ( ) u � i = fi( ) u � = au u � = Au a � a � 0 Aij � i j Re �(A) � 0 1 equation u, tu, t�f/�u �f /�uFor an ordinary differential equation u� = f (u, t), good codes will increase the accuracy (and keep stability) far beyond the O(�t) error in Euler’s methods. You can 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 Strang� �  �  �  �  �  � c�2006 Gilbert Strang CHAPTER 5. INITIAL VALUE PROBLEMS Stiff Differential Equations First we call attention to one extra question: Is the equation stiff ? Let me begin with a made-up example to introduce stiffness and its effects: v(t) = e−t + e−99t controls decay controls �t The step �t is 99 times smaller because of e −99t that disappears anyway Those decay rates −1 and −99 could be the eigenvalues of A, as in Example 1. Example 1 dv −50 49 v v(0) 2 = with = . (1)dt w 49 −50 w w(0) 0 −99tThe solution has v(t) = e−t + e−99t and w(t) = e−t − e . The time scales are different by a factor of 99 (the condition number of A). The solution will decay at the slow time scale of e−t, but computing e−99t may 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 happens and 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: du −Ku dui ui+1 − 2ui + ui−1Method of Lines = has = (�x)2 . (2)dt (�x)2 dt This comes from the heat equation �u/�t = �2 u/�x2, by discretizing only the space derivative. 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 about a = −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 the space variables first produces a large system that can be given to an ODE solver. (We have ordinary differential equations in time.) If it wants to, that solver can vary the time step �t and even the discretization formula as u(t) speeds up or slows down.PSfragreplacementsc5.1. FINITE DIFFERENCE METHODS �2006 Gilbert Strang This method splits the approximation of a PDE into two parts. Finite differences/finite elements in earlier chapters produce the first part (discrete in space). The upcoming stability-accuracy analysis applies to the second part (discrete in time). This idea is very simple and useful, even if it misses the good methods later in this chapter that take full advantage of space-time. For the heat equation ut = uxx, the useful fact utt = uxxxx allows us to cancel space errors with time errors—which we won’t notice when they are separated in the semidiscrete method of lines. Forward Euler and Backward Euler The equation u� = f(u, t) starts from an initial value u(0). The key point is that the rate of change u� is determined by the current state u at any moment t. This model of reality, where all the history is contained in the current state u(t), is a tremendous success throughout science and engineering. (It makes calculus almost as important as linear algebra.) But for a computer, continuous time has to change to discrete time. One differential equation allows many difference equations! The simplest method (Euler is pronounced “Oiler”) uses a forward difference: Un+1 − Un Forward Euler = f(Un, tn) is Un+1 = Un + �t fn . (3)�t Over each �t interval, the slope of U doesn’t change. Figure 5.1 shows how the correct solution to u� = au follows a smooth curve, while U(t) is only piecewise linear. A better method (higher accuracy) will stay much closer to the curve by using more information than the one slope fn = f(Un, tn) at the start of the step. u�t u = e�t U t (0) = 1 = 1 + �(forward) t t �t u = e�t U = 1 1 − �t(backward) Figure 5.1: Forward Euler and Backward Euler for u� = u. One-step errors � 1 2 (�t)2 . Backward Euler comes from using fn+1 at the end of the step, when t = tn+1: Un+1 − Un Backward Euler = f(Un+1, tn+1) is Un+1 − �t fn+1 = Un . (4)�t This is an implicit method. To find Un+1, we have to solve equation (4). When f is