Chapter 9

1. Chapter 9 Ordinary Differential Equations: Initial-Value Problems Lecture (II)1 1 Besides the main textbook, also see Ref.: “Applied Numerical Methods with MATLAB for Engineers and Scientists”, Steven Chapra, 2nd ed., Ch. 20, McGraw Hill, 2008. Dr. Jie Zou PHY 3320

2. Outline • Error Analysis for and Stability of Euler’s Method • Improvements of Euler’s method • (I) Heun’s Method • (II) The Midpoint Method (next lecture) Dr. Jie Zou PHY 3320

3. Error analysis for Euler’s method • Two types of errors involved in the numerical solution of an ODE: • (1) Truncation errors: Errors due to the nature of the numerical techniques used to approximate values of y. Two parts: • (I) Local truncation error: Error due to an application of the numerical method over a single step. • (II) Propagated truncation error: Error due to the approximations produced during the previous steps. • Global truncation error (total error): The sum of (I) and (II). • (2) Roundoff errors: Errors due to the limitation of the computer to keep an infinite number of significant digits. Dr. Jie Zou PHY 3320

4. Analysis of the magnitude and properties of the truncation error • Derive Euler’s method from the Taylor series expansion: • Taylor series expansion of y(t) about a point ti: • h = ti+1 – ti: step size; Rn = y(n+1)()hn+1/(n+1)!: the reminder term;  lies somewhere between ti and ti+1. • Evaluating y(t) at point ti+1, we have: • Using the general form of an ODE, dy/dt = f(t, y), we have: • O(hn+1) : For the nth-order approximation, the local truncation error  hn+1.

5. True and approximate local truncation error • Comparing the following two equations: (1) Taylor series expansion (2) Euler’s method True local truncation error When h is sufficiently small, the higher-order terms are usually negligible. Approximate local truncation error Dr. Jie Zou PHY 3320

6. Global truncation error • It can be shown that the global truncation error is O(h), i.e.,  step size h. • Note: • The global error can be reduced by decreasing the step size h. • If the solution of the ODE is linear, there is no error involved in Euler’s method – Euler’s method is a first-order method. Dr. Jie Zou PHY 3320

7. Stability of Euler’s method • “Unstable”: A numerical solution is said to be unstable if errors grow exponentially for a problem for which there is a bounded solution. • The stability of a particular numerical method depends on three factors: • (1) Differential equation • (2) Numerical method • (3) Step size • An example: Let’s look at a simple ODE dy/dt = -ay, y(0) = y0 • Calculus: y = y0e-at • Euler’s method: yi+1 = yi + (dyi/dt)h = yi(1 – ah) • If h > 2/a, then |1 – ah| > 1 and |yi|   as i   - Euler’s method is conditionally stable. Dr. Jie Zou PHY 3320

8. Improvement of Euler’s method: (I) Heun’s method • Heun’s method: a predictor-corrector approach • (1) Predicator: Use the slope at the beginning of an interval to predict that at the end of the interval. • (2) Corrector: Use the averaged slope to obtain an improved estimate of y at the end of the interval. • Predictor: yi+10 = yim + f(ti,yi)h • Corrector: yi+1j= yim + {[(f(ti,yim) + f(ti+1, yi+1j-1)]/2}h, for j = 1, 2, … • The corrector equation is iterative; j indicates the jth iteration. • Stopping criterion: |a| = |(yi+1j – yi+1j-1)/yi+1j| x 100% (a) Predictor Dr. Jie Zou PHY 3320 (b) Corrector

9. Example: Heun’s method • Example 20.2 (Ref.): Use Heun’s method with iteration to integrate y’ = 4e0.8t – 0.5y from t = 0 to 4 with a step size of 1. The initial condition at t = 0 is y = 2. Employ a stopping criterion of 0.00001% to terminate the corrector iterations. • (1) By hand. • (2) Write an M-file. A copy of the code will be handed out later. Dr. Jie Zou PHY 3320

10. Results With Iteration Without Iteration Dr. Jie Zou PHY 3320 Table 20.2 (Ref. by Chapra)