MATH/APPM 4650
Class notes, April 10
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 usually 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. There are two possibilities:
     1. The methods have different orders of error. Then you use the higher-order method as the approximation of the exact solution to estimate the error of the lower-order method. (This is what's usually done for Runge-Kutta methods.)
     2. The methods have the same order of error with different coefficients. (You need to know the coefficients for this to work.) Then knowing (w₂-y)/(w₁-y) = C, you can solve for y given w₁ and w₂ and C. Thus you can estimate the errors in both w₁ and w₂. (This is what's usually done for Adams methods.)
  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 computes six values of the function, uses a fourth-order Runge-Kutta for the actual prediction, and uses a fifth-order Runge-Kutta to estimate the error. 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 the coefficients are relatively easy to compute. It's discussed in Section 5.7.
  
  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).
  
  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.
• (5.9): Systems
  Every method we have discussed works exactly the same way for a differential equation x'(t) = f(t, x(t) where x = (x₁, . . ., x_d) is a d-dimensional vector and f(t, x) = (f₁(t, x₁, . . ., x_d), . . ., f_d(t, x₁, . . ., x_d).
  If you have higher-order differential equations, you can always add new variables to get a first-order system. For example,
  x''(t) = -f(t, x(t)) would become the system x'(t) = v(t),     v'(t) = -f(t, x(t)).
• (5.10) Stability and convergence
  Definition: A difference scheme for a differential equation is convergent if the maximum difference between the scheme-produced function and the actual solution goes to zero as the number n of steps goes to infinity.
  The Lax equivalence theorem says that a difference scheme is convergent if and only if it is both consistent and zero-stable.
  Check consistency by computing the local truncation error: assuming y(t_k) = w_k, . . ., y(t_k-m+1) = w_k-m+1, find a Taylor series for
  1/h [ y(t_k+1) - w_k+1 ].
  You want this to go to zero as h goes to zero.
  Check zero-stability by finding the characteristic polynomial of the method and making sure it satisfies the root condition: all roots less than or equal to one in absolute value, and no root with absolute value one is repeated.
  To say more, you can distinguish strong stability (the only root of size one is λ=1, which has to be a root) vs. weak stability (there is at least one root of size one besides λ=1).
• (5.11) Stiff equations
  A rapidly decaying solution is hard to approximate numerically, unless you take a much smaller step size than you'd like.
  To understand whether it's happening in your method, you try a test equation y' = λy for a complex number λ.
  You get a stability function w_k+1 = Q(hλ)w_k. The method is A-stable if |Q(z)| < 1 whenever Re(z) < 0.
  Among methods we've studied, only the implicit Euler and the implicit trapezoid method are A-stable.