3 Runge-Kutta Methods - IIT

1 3 Runge-Kutta MethodsIn contrast to the multistep Methods of the previous section, Runge-Kutta methodsare single-step Methods however, with multiplestagesper step. They are motivatedby the dependence of the Taylor Methods on the specific IVP. These new Methods donot require derivatives of the right-hand side functionfin the code, and are thereforegeneral-purpose initial value problem solvers. Runge-Kutta Methods are among themost popular ODE solvers. They were first studied by Carle Runge and Martin Kuttaaround 1900. Modern developments are mostly due to John Butcher in the Second-Order Runge-Kutta MethodsAs always we consider the general first-order ODE systemy (t) =f(t,y(t)).(42)Since we want to construct a second-order method, we start with the Taylor expansiony(t+h) =y(t) +hy (t) +h22y (t) +O(h3).

2 The first derivative can be replaced by the right-hand side of the differential equation(42), and the second derivative is obtained by differentiating (42), ,y (t) =ft(t,y) +fy(t,y)y (t)=ft(t,y) +fy(t,y)f(t,y),with Jacobianfy. We will from now on neglect the dependence ofyontwhen itappears as an argument tof. Therefore, the Taylor expansion becomesy(t+h) =y(t) +hf(t,y) +h22[ft(t,y) +fy(t,y)f(t,y)] +O(h3)=y(t) +h2f(t,y) +h2[f(t,y) +hft(t,y) +hfy(t,y)f(t,y)] +O(h3).(43)Recalling the multivariate Taylor expansionf(t+h,y+k) =f(t,y) +hft(t,y) +fy(t,y)k+..we see that the expression in brackets in (43) can be interpreted asf(t+h,y+hf(t,y)) =f(t,y) +hft(t,y) +hfy(t,y)f(t,y) +O(h2).Therefore, we gety(t+h) =y(t) +h2f(t,y) +h2f(t+h,y+hf(t,y)) +O(h3)or the numerical methodyn+1=yn+h(12k1+12k2),(44)56withk1= f(tn,yn),k2=f(tn+h,yn+hk1).

3 This is theclassical second-order Runge-Kutta method. It is also known asHeun smethodor theimproved Euler known asstagesof the Runge-Kutta method. Theycorrespond to different estimates for theslopeof the solution. Note thatyn+hk1corresponds to an Euler step with stepsizehstarting from (tn,yn). Therefore,k2corresponds to the slope of the solution one would get by taking one Euler stepwith stepsizehstarting from (tn,yn). The numerical method (44) now consistsof a single step with the average of the notation used here differs slightly from that used in the Iserles book. Therethe stages are defined differently. I find the interpretation in terms of slopes also saw earlier that the classical second-order Runge-Kutta method can beinterpreted as a predictor-corrector method where Euler s method is used as thepredictor for the (implicit) trapezoidal obtain general explicit second-order Runge-Kutta Methods by assumingy(t+h) =y(t) +h[b1 k1+b2 k2]+O(h3)(45)with k1=f(t,y) k2=f(t+c2h,y+ha21 k1).

4 Clearly, this is a generalization of the classical Runge-Kutta method since the choiceb1=b2=12andc2=a21= 1 yields that is customary to arrange the coefficientsaij,bi, andciin a so-calledRunge-KuttaorButcher tableauxas , the Butcher tableaux for the classical second-order Runge-Kutta methodis00 011 Runge-Kutta Methods are characterized by a strictly lower triangular ma-trixA, ,aij= 0 ifj i. Moreover, the coefficientsciandaijare connected by theconditionci= j=1aij, i= 1,2,.., .This says thatciis the row sum of thei-th row of the matrixA. This condition isrequired to have a method of order one, , for consistency. We limit our discussionto such Methods , for an explicit second-order method we necessarily havea11=a12=a22=c1= 0. We can now study what other combinations ofb1,b2,c2anda21in (45) giveus a second-order method.

5 The bivariate Taylor expansion yieldsf(t+c2h,y+ha21 k1) =f(t,y) +c2hft(t,y) +ha21fy(t,y) k1+O(h2)=f(t,y) +c2hft(t,y) +ha21fy(t,y)f(t,y) +O(h2).Therefore, the general second-order Runge-Kutta assumption (45) becomesy(t+h) =y(t) +h[b1f(t,y) +b2{f(t,y) +c2hft(t,y) +ha21fy(t,y)f(t,y)}] +O(h3)=y(t) + (b1+b2)hf(t,y) +b2h2[c2ft(t,y) +a21fy(t,y)f(t,y)] +O(h3).In order for this to match the general Taylor expansion (43) we wantb1+b2= 1c2b2=12a21b2= , we have a system of three nonlinear equations for our four unknowns. Onepopular solution is the choiceb1= 0,b2= 1, andc2=a21=12. This leads tothemodified Euler method(sometimes also referred to as themidpoint rule, see thediscussion in Section below)yn+1=yn+hk2withk1=f(tn,yn)k2=f(tn+ h2,yn+h2k1).Its Butcher tableaux is of the form00 0121200 choiceb1= 1,b2= 0 leads to Euler s method.

6 However, since nowc2b26=12anda21b26=12this method does not have second-order explicit Runge-Kutta Methods are of the formyn+1=yn+h j=1bjkjwithk1=f(tn,yn)k2=f(tn+c2h,yn+a21 hk1)..k =f(tn+c h,yn+h 1 j=1a ,jkj).Determination of the coefficients is rather complicated. We now describe (withoutderivation) the most famous Runge-Kutta Fourth-Order Runge-Kutta MethodsThe classical method is given byyn+1=yn+h[k16+k23+k33+k46](46)withk1=f (tn,yn)k2=f(tn+h2,yn+h2k1)k3=f(tn+h2,yn+ h2k2)k4=f(tn+h,yn+hk3).Its Butcher tableaux is of the form00 0 0 012120 0 0120120 010 0 1 local truncation error for this method isO(h5). It is also important to notethat the classical fourth-order Runge-Kutta method requires four evaluations of thefunctionfper time saw earlier that in each time step of the second-order Runge-Kuttamethod we need to perform two evaluations off, and for a fourth-order method thereare four evaluations.

7 More generally, one can observe the situation described in Table offper time step234567891011maximum order achievable2344566778 Table 4: Efficiency of Runge-Kutta data imply that higher-order (>4) Runge-Kutta Methods are relativelyinefficient. Precise data for higher-order Methods does not seem to be known. However,certain higher-order Methods may still be appropriate if we want to construct a Runge-Kutta method which adaptively chooses the step size for the time step in order to keepthe local truncation error small (see Section 5). Connection to Numerical Integration RulesWe now illustrate the connection of Runge-Kutta Methods to numerical before, we consider the IVPy (t) =f(t,y(t))y(t0) =y0and integrate both sides of the differential equation fromttot+hto obtainy(t+h) y(t) = t+htf( ,y( ))d.

8 (47)Therefore, the solution to our IVP can be obtained by solving the integral equation(47). Of course, we can use numerical integration to do the left endpoint method baf(x)dx b an =hn 1 i=0f(xi)on a single interval, , withn= 1, anda=t,b=t+hwe get t+htf( ,y( ))d t+h t1f( 0,y( 0))=hf(t,y(t))since 0=t, the left endpoint of the interval. Thus, as we saw earlier, (47) isequivalent to Euler s the trapezoidal rule baf(x)dx b a2[f(a) +f(b)]witha=tandb=t+hgives us t+htf( ,y( ))d h2[f(t,y(t)) +f(t+h,y(t+h))].60 The corresponding IVP solver is thereforeyn+1=yn+h2f(tn,yn) +h2f(tn+1,yn+1).Note that this isnotequal to the classical second-order Runge-Kutta methodsince we have ayn+1term on the right-hand side. This means that we have animplicitmethod. In order to make the method explicit we can use Euler s methodto replaceyn+1on the right-hand side byyn+1=yn+hf(tn,yn).

9 Then we end up with the methodyn+1=yn+h2f(tn,yn) +h2f(tn+1,yn+hf(tn,yn))oryn+1=yn+h[12k1+ 12k2]withk1=f(tn,yn)k2=f(tn+h,yn+hk1), , the classical second-order Runge-Kutta midpoint integration rule baf(x)dx (b a)f((a+b)/2)witha=t,b=t+ 2hgives us t+2htf( ,y( ))d 2hf(t+h,y(t+h)).Thus, we have theexplicit midpoint ruleyn+2=yn+ 2hf(tn+1,yn+1).This is not a Runge-Kutta method. It is an explicit 2-step method. In thecontext of PDEs this method reappears as theleapfrog method. As mentionedabove, sometimes the modified Euler methodyn+1=yn+hf(tn+h2,yn+h2f(tn,yn)).is called the midpoint rule. This can be explained by applying the midpointintegration rule witha=tandb=t+h. Then we have t+htf( ,y( ))d hf(t+h2,y(t+h2)).If we representy(t+h2) by its Euler approximationy(t)+h2f(t,y), then we arriveat the modified Euler method stated s rule yields the fourth-order Runge-Kutta method in case there is nodependence quadrature leads to so-calledGauss-Runge-KuttaorGauss-Legend remeth-ods.

10 One such method is the implicit midpoint ruleyn+1=yn+hf(tn+h2,12(yn+yn+1))encount ered earlier. The Butcher tableaux for this one-stage order two methodis given that the general implicit Runge-Kutta method is of the formyn+1=yn+h j=1bjkjwithkj=f(tn+cjh,yn+hj i=1aj,iki)for all values ofj= 1,.., . Thus, the implicit midpoint rule corresponds toyn+1=yn+hk1withk1=f(tn+h2,yn+h2k1) obviously an implicit generalimplicit Runge-Kutta methodsexist. However, their construction ismore difficult, and can sometimes be linked tocollocation Methods . Some detailsare given at the end of Chapter 3 in the Iserles