MATH/APPM 4650
Class notes, April 16
Review of Chapter 5

Everything we have discussed in this chapter is about techniques for solving y'(t) = f(t, y(t)),     y(0) = α
For the most part y(t) is a number-valued function of t, but in Section 5.9 we allowed y(t) to be a vector-valued function of t.
This does not change anything at all about the techniques, since all we'd ever do is change number addition to vector addition.

• (5.1): As long as f is continuous (e.g., not f(t, y) = y/t) and differentiable (e.g., not f(t, y) = sqrt(y)), the differential equation will have a unique solution defined for some time interval (t₀ - ε, t₀ + ε).
  
  Thus when we're doing approximations, there is always a theoretical "exact" solution we can compare to (even if we don't know precisely what it is).
  
  
• (5.2): Euler's method. w_i+1 = w_i + h f(t_i, w_i), where t_i = h i and w₀ = α.
  • This is definitely the simplest method.
  • The error in the final approximation is O(h). (This is also the local truncation error.)
  • We proved explicitly that Euler's method converges as h approaches zero.
  • Euler's method is a special case of all other methods (it's the Taylor method of order one, the first-order Runge-Kutta method, and the one-step Adams-Bashforth method).
  
  
  
• (5.3): Higher-order Taylor methods.
  Basic idea of the method: y(t + h) ~ y(t) + h y'(t) + h²/2 y''(t) + ... + h^m/m! y^(m)(t) with local truncation error O(h^m).
  We already know y'(t) = f(t, y(t)), so we can compute y''(t) and all higher derivatives using the differential equation.
  
  For example, 2nd-order Taylor method is w_i+1 = w_i + h f(t_i, w_i) + h²/2 [ f_t(t_i, w_i) + f(t_i, w_i) f_y(t_i, w_i) ]
  
  
  • Error is always O(h^m).
  • Easy to remember and derive.
  • Fine if f is really simple, so computing derivatives is no big deal.
  • Really difficult for general functions, where derivatives are hard to compute.
  • Error coefficients usually not as good as Runge-Kutta. (Remember homework.)
  
  
  
• (5.4): Runge-Kutta methods.
  Basic idea: If you expand y(t + h) in a Taylor series in h, and you also expand the following expression in a Taylor series in h, y(t) + h f[ t + h/2, y(t) + h/2 f(t, y(t)) ] they are the same up to the h² term.
  
  So we get a method that's at least as good as the Taylor method, but we only need to compute values of f (not its derivatives).
  
  So the algorithm is k₁ = f(t_i, w_i),     k₂ = f(t_i + h/2, w_i + h/2 k₁)
  w_i+1 = w_i + h k₂ This is one of the simplest second-order Runge-Kutta methods, but there are many other Runge-Kutta methods (even for second-order).
  
  • There are three popular second-order Runge-Kutta methods which use k₁ and some k₂: Midpoint, Modified Euler, and Heun's. All have error O(h²) but with different constants. (Heun's is best.)
  • The most common Runge-Kutta method is the fourth-order Runge-Kutta (FORK): it uses k₁, k₂, k₃, and k₄, and has error O(h⁴).
  • To get O(h⁵), you need six k's (not five!). In general, after fourth-order, Runge-Kutta methods start to get less efficient.
  • Runge-Kutta methods tend to be better than the Taylor methods they approximate, since you often part of the higher-order terms for free (but not all of them).
  • The main problem with Runge-Kutta is that you need to evaluate the function several times each step. If the function is complicated, this will be more time-consuming than an Adams method.
  • However Runge-Kutta tends to give better accuracy than Adams methods.
  
  Typical algorithm (Heun's method), given f, a, b, N, α
  h = (b-a)/N
  w = α
  t = a
  FOR i FROM 1 to N DO
     k1 = f(t, w)
     k2 = f(t+2*h/3, w+2*h/3*k1)
     w = w + h/4*(k1 + 3*k2)
     t = t + h
  END DO
  OUTPUT w
  
  
• (5.6): Adams-Bashforth methods.
  Basic idea: Instead of just using one previous step, use several previous steps.
  And instead of using Taylor series to derive the method, use interpolating polynomials to approximate an integral: y(t + h) = y(t) + ∫_t^t+h f(τ, y(τ)) dτ
  
  Typical Adams-Bashforth method (two-step): w_i+1 = w_i + h/2 [3f(t_i, w_i) - f(t_i, w_i)]
  
  • An m-step Adams-Bashforth method always has order O(h^m).
  • Since it's a multistep method, you need to get it started using some other method. Always use a method of at least the same order.
    For example, using a two-step Adams-Bashforth method, one needs to generate w₁. One could do this using Heun's method (order two) or standard Runge-Kutta (order four) but never Euler's method (order one).
  • The main benefit of Adams-Bashforth over Runge-Kutta is that you only need one new function evaluation per step to get the same order of error. If the function is very complicated, this is a big advantage.
  • However, a Runge-Kutta method will generally do better than an Adams-Bashforth method of the same order. (Smaller coefficient.) Recall your homework.
  • Unlike Runge-Kutta, there is only one Adams-Bashforth method for a given order.
  Typical algorithm (two-step Adams-Bashforth method, using Heun's method to start), given f, a, b, N, α
  h = (b-a)/N
  w = α
  t = a
  oldf = f(t, w)
  w = w + h/4 * (oldf + 3*f(t+2*h/3, w+2/3*h*oldf))
  t = t + h
  newf = f(t, w)
  FOR i FROM 2 to N DO
     w = w + h/2*(3*newf - oldf)
     t = t + h
     oldf = newf
     newf = f(t, w)
  END DO
  OUTPUT w
  
  
• (5.6): Adams-Moulton methods.
  Basic idea: Do the same thing as in Adams-Bashforth methods, but allow the right-hand side to depend on f(t_i+1, w_i+1). By doing this, we need one fewer step to get the same error.
  Typical Adams-Moulton method (one-step, trapezoid rule) w_i+1 = w_i + h/2 [f(t_i, w_i) + f(t_i+1, w_i+1)].
  
  • An m-step Adams-Moulton method always has error O(h^m+1) (one better than all other methods).
  • The drawback is having to solve for w_i+1. This can only be done in very special cases.
  • Typically, one would use an Adams-Moulton method as a "corrector" for an Adams-Bashforth method.
    For example, use the two-step Adams-Bashforth as a predictor, and a one-step Adams-Moulton as a corrector. (Both methods have error O(h²).)
  Typical algorithm (one-step Adams-Moulton predictor-corrector method, using Heun's method to start), given f, a, b, N, α
  h = (b-a)/N
  w = α
  t = a
  oldf = f(t, w)
  w = w + h/4 * (oldf + 3*f(t+2*h/3, w+2/3*h*oldf))
  t = t + h
  newf = f(t, w)
  FOR i FROM 2 to N DO
     wapprox = w + h/2 * (3*newf - oldf)
     t = t + h
     oldf = newf
     newf = f(t, wapprox)
     w = w + h/2 * (oldf + newf)
     newf = f(t, w)
  END DO
  OUTPUT w
  
  
• (5.5 and 5.7): Adaptive methods.
  Adaptive methods shrink the size of h when the solution is changing rapidly, and keep it the same when the solution is changing slowly. Thus they're always more efficient than non-adaptive methods (where the same h is used the whole time).
  
  There are two basic ingredients for an adaptive method:
  1. Using two different methods to generate the new value, so you can compare them and estimate the error.
  2. Checking whether the error is smaller than some given tolerance, and shrinking h and starting over if it is.
  
  The most common adaptive method is the Runge-Kutta-Fehlberg method, which uses a fourth-order Runge-Kutta for the actual prediction and a fifth-order Runge-Kutta to estimate the error. This is usually abbreviated "rkf45." It's discussed in Section 5.5.
  
  The most common Adams version of an adaptive method is the fourth-order Adams-Bashforth used as "predictor" and a third-order Adams-Moulton used as "corrector." Both methods are order four, but by analyzing it a little more carefully, you can also get an error estimate from this. It's discussed in Section 5.7. I don't think there's a standard abbreviation for it.
  
  Since these methods depend on making lots of heuristic choices, it's better to use prepackaged versions.
  
  Here is how you'd use a prepackaged method to solve y'(t) = y(t)² + t² with initial condition y(0) = 0, if what you want is y(0.5).
  
  Maple:
  soln := dsolve({diff(y(t),t) = y(t)^2 + t^2, y(0)=0}, numeric, method=rkf45);
  soln(0.5);
  Maple prefers not using Adams methods adaptively.
  
  MATLAB:
  Assume a file exists with the name F.m and containing the following:
  function dy=F(t,y)
  dy=y^2+t^2;
  
  Then at the command prompt, use (for RKF45):
  [times,soln]=ode45(@F, [0 0.5], [0]);
  [times(length(times)),soln(length(soln))]
  
  For Adams method, you'd use:
  [times,soln]=ode113(@F, [0 0.5], [0]);
  [times(length(times)),soln(length(soln))]
  
  Mathematica (Runge-Kutta):
  soln=NDSolve[{y'[t]==y[t]^2+t^2,y[0]==0},y,{t,0.5}, Method->ExplicitRungeKutta];
  y[0.5]/.soln
  
  Using an Adams method:
  soln=NDSolve[{y'[t]==y[t]^2+t^2,y[0]==0},y,{t,0.5}, Method->Adams];
  y[0.5]/.soln
  
  Mathematica has proprietary versions of these algorithms, so it doesn't use the standard fourth-order method, for example.