MA332 lecture notes, Fall 1998 Week 4, Tuesday Reading: F&B Chapter 5 Initial-value problems 1) Taylor methods -- we always start with Taylor a) Euler is a special case of Taylor 2) Runge-Kutta 3) Predictor-Corrector -- Adams Bashforth 4) Adaptive Runge-Kutta-Fehlberg 5) Systems of equations (maybe) IVPs in general First order IVP dy/dt = f(t, y), with initial value y(a) = y0 Generalize to a system of equations: dy1/dt = f1(t, y1, y2); dy2/dt = f2(t, y1, y2) Or to a higher order IVP: y^(n) = f(t, y, y', y'', ... , y^(n-1)); But for now let's stick with one first-order IVP: dy/dt = 1 + t * sin (ty) (this is not the example from the book) Check out the handout. 1) at each point in the plane (t, y) we can find a slope, which we can indicate with an arrow 2) the initial condition tells us where to start 3) the arrows tell us where to go Immediately, without really trying, we have our first numerical method, Euler's method. 1) start with t0 = a, w0 = f(a), where f(a) is the initial condition 2) choose a step size h 3) w1 = w0 + h * f(t0, w0) 4) repeat until you have answered the question, whatever it was Is there any reason to think that's a good idea? Well, remember Taylor's theorem, which says y(t+h) = y(t) + h y'(t) + h^2/2 y''(t) + ... = y(t) + h y'(t) + h^2/2 y''(xi) for some xi in (x, x+h) Well, that looks real familiar, since we know that y'(t) = f(t, y). Furthermore, that indicates that the error term is h^2/2 y''(xi). Local and global error --------------------- Of course, that's a local error (the error associated with a single step). The global error is the accumulated error after many steps. If the functions are well-behaved (meaning that the direction of the arrows changes slowly as we move around the plane), then the global error can be bounded. 1) good case: global error is one order worse than local error (which makes sense, since we are taking 1/h steps with h^2 error each 2) bad case: global error is unbounded Handout shows an example of the bad case. On page 182, the two conditions are: df/dt (t, y) < L and y''(t) < M If both are true then the gloval error is O(h), which is only one order worse than the local error, which is really the best we can hope for considering that the number of steps increases linearly with h. Generalized Taylor method ------------------------ Given y' we can find y'' and evaluate: y(t+h) = y(t) + h y'(t) + h^2/2 y''(t) + K h^3 y'''(xi) Local error is h^3, global is h^2, given same assumptions, plus bounded y'''. Increasing the order calls for more analytic differentiation, and more continuity requirements. Runge-Kutta methods ------------------- Similar to Gaussian integration in the sense that they look for magic places to evaluate the function, in order to make errors cancel. Various derivations of this kind lead to... Midpoint method: 1) eval f(ti, wi) and go h/2 in that direction to get t^ and w^ 2) eval f(t^, w^) and go h in that direction Modifier Euler method: 1) eval f(ti, wi) and go h in that direction to get t^ and w^ 2) take the average of f(ti, wi) and f(t^, w^) and go h in that direction or Heuh's method, which is getting silly now because there is a free parameter here and we can all have our own method named after us. All of these are O(h^3) local error, yielding O(h^2) globally, maybe. 4th Order Runge Kutta is a little more work to derive, but it works pretty well in practice: 1) k1 = h * the slope at ti, wi 2) follow k1 out to the halfway mark and k2 = h * slope there 3) follow k2 out to the halfway mark and k3 = h * slope there 4) follow k3 all the way across and k4 = h * slope there 5) the next w = w + (k1 + 2k2 + 2k3 + k4) / 6 Local error O(h^5), global 4th order. Interesting that the middle two are most reliable. Multistep methods ----------------- So far, estimates of the next w depend only on the current value of w and function evaluations at or to the right of t. Multistep method use previous values of w as well. Adams-Bashforth explicit methods frac i-4 i-3 i-2 i-1 i i+1 local error h/2 -1 3 h^3 h/12 5 -16 23 h^4 h/24 -9 37 -59 55 h^5 h/720 251 -1274 2616 -2774 1901 h^6 How do we get started? Usually use R-K method with the same error order to get the first n values. Adams-Moulton implicit methods frac i-4 i-3 i-2 i-1 i i+1 local error h/12 -1 8 5 h^4 h/24 1 -5 19 9 h^5 h/720 -19 106 -246 646 251 h^6 In order to use the implicit methods, we have to solve algebraically for the next w in terms of its predecessors. When this can be done analytically, the improvement in performance is often worthwhile. When it cannot, we might be stuck solving the equation numerically, which is likely to be silly. On the other hand, we could use an explicit method to approximate the next w, and then use the implicit method to improve the estimate. This approach is called predictor-corrector.