Example 3.36. Consider the IVP (3.88) y' =y y(0) =1 We know that the exact solution is y(x) = e". Let us run Euler's method, say on the interval [0, 1]: (3.89) Un+1 = Un + hF(In, y(x.)) (3.90) Un+1 = yn + hyn (3.91) = (1 + h)yn. Iterating this relation we obtain Un = (1 + h)yn-1 = (1+ h) yn-2= ... = (1+ h)"yo- Since yo = 1, we conclude (3.92) Un = (1 + h)" (n =1, 2, ...) This simple example allows us to express the approximation error in the nth step explicitly: (3.93) En = lyn - y(In)| = (1+h)" _ ehnj (Recall we want that the global approximation error max,=0....,N |Un - y(T,)| converges to zero.) In particular, if we recall that h = N, (3.94) EN = 1(1 + N-1 )N - el The statement of Theorem 3.34 thus in particular shows that the global approximation error does converge to zero, in particular ex - 0, so we have recovered the familiar formula (3.95) lim (1 + ) =e N-+00 as a special case of the convergence of Euler's method. Remark 3.37 (*). If you haven't done it before, it is a good exercise to show by explicit calculation using the binomial theorem that (1 + )~converges to e = = Exercise 3.38. Consider the IVP (3.96) y's 2y + 1 y(0) = 1 The exact solution is y(@) = ;ex - (i) Implement Euler's method for this IVP on the interval [0, 1] and plot the approximate solutions for h = 1/10, h = 1/20, h = 1/50, h = 1/100; also plot the exact solution for comparison. (ii) Plot the approximation error EN = lyn - y(IN)| against h as h tends to 0. (iii*) Calculate yN explicitly and show that EN - 0 as N + co(without using a generalLet us introduce some abbreviations for the rest of this computation: (3.139) T= In, y = y(on), y' = y'(a.), y" = y"(x.), F = F(In, Un) (whenever we use arguments evaluated at ~ = 2,, then arguments are omitted from nota- tion) Then substituting the definition of d, (3.140) Th = y'+ thy"-bF - byF(x + ch,y + ahF) + 0(h?) To analyze this further we notice that by the ODE, (3.141) y =F and by the chain rule, (3.142) y" = OF + y'd,F = OF + Fa, F. For a smooth function F of two variables we have by the two-dimensional Taylor formula: (3.143) F(x + Ax, y + Ay) = F(x, y) + AxOF(x, y) + Ayay F(x, y) + O(h2), whenever [Ax| + |Ay) 0. (iv*) Make your implementation from (i) more general so that solves a general IVP y' = F(x, y), y(ro) = yo on an interval [To, To + a] (i.e. you should write a function/method whose inputs are the function F, and values co, yo, a, as well as the parameter N for the number of steps). The output should be the list of computed approximations yn. You can also write routines for plotting and testing against a given exact solution. Test it on various examples of IVPs seen in class. Hint 1: In implementations, it is more practical to make N the independent parameter rather than h. One can then use h = Hint 2: You can use Python 3 with numpy for vectorization and matplotlib for plotting. (Alternatively, Matlab or Octave.) Solution to (iii*). Euler's method reads (3.97) Un+1 = Un + h(2yn + 1) = (1 +2h)yn + h. This is a linear recurrence relation and we want to solve for y,. There is a systematic way of doing this. Note that if the term +h would not be there, then y. would be (1 + 24)". This motivates the substitution (3.98) an = Un(1+ 2h)-". The recurrence becomes (3.99) Entl = anth(1+2h)-(+1) Subtract 2, and sum over n = 0, . . ., N - 1: N-1 N-1 (3.100) [( anti - an) = [h(1 + 2h) -(n+1) 71=0 n=0 On the left we have a telescoping sum, where only the last and the first terms survive: N-1 (3.101) 2N - 20 = h(1+ 2h )-(r+1) n=0 Substitute back and simplify (recall 20 = yo = 1): N-1 (3.102) UN = (1+2h) +h > (1+ 2h ) N-n-1 71=0 Reindexing by replacing n with (N - 1) - n we get N-1 (3.103) yN = (1+24) ~ +h > (1+2h)" 1=0Summing the geometric sum, (3.104} \"in +2h)\" = W. Plugging this back in, we get chFinal answer (3.105} ya; = {1 + 211)\" + till + 2h)\" _ 1) = go + H" . Recalling that lim,-.'r_,m[1 + %}'V = ei, we have (3.106} Ems. = as = you. III Euler's method is rarely used in practice, because much more accurate methods are available. We will look in more detail at two classes of such methods: onestep methods and linear multistep methods. Euler's method is a special case of both. 3.4. One-step methods. Recall that we are trying to solve IVPs of the form if = Fix, 1;), 3.10? l l { Bil-T0] = so- with F smooth. We assume that a unique solution 3; = y{:r]l exists on an interval [3.1], xnIe] .6 Denition 3.39. A general onestep method takes the form {3-108} yn+l : yrs + h(D{-Inayn1 h]! where n = U,...,N, at,, = rn+nh and h 2: U, N = (ll. The map (D denes what happens from one step to the next. In principle, (D could be any function, but the method is only useful if it actually gives an approximate solution to the IVP [3.1'07}. Thus, 'I! should involve the function F. For simplicity, we'll always assume that ii is smooth in all variables. Intuitiv'ely, if y = y[::r:} is assumed to be the exact solution, then @[13 y, h) represents an estimator for the value of the derivative of the solution y across the interval [3,513 + h]. Example 3.40 (Euler's method]. If @[I,y,h) = F{;r,y], then we recover Euler's method. FI[:.r._y]I is the value of the derivative of y at the left endpoint of the interval [33,1 + h] this may be an underestimate, or overestimate, depending on whether the derivative will increase or decrease (Le. depending on the second derivative of y]. Example 3.41 {Heun's method}. The one-step method given by (3-109) so. y. m = am, y) + Fir + h. y + we. ym is called Heart's method. Here the idea is to take an average of derivatives at the left- and right endpoints of the interval [33,1' + h]. At the left endpoint, we use F[3,y} (which (*] \"'0 also make a silent technical assumption: we x a compact set K containing the solution curve (2', M23\". Constants throughout may depend on K. equals y'{s'j if y is exact} and at the right endpoint we use F[:t: + lay + hF{.c,y)) as an approximation of 3411? + h] {using one step of Euler's method to estimate the value of y(:e + h) ). Emmpie 3.42 {Midpoint method). The onestep method given by {3.110} on, y, 3:} = Fe. + is. y + %hF(s, y}} is called midpoint method. 'We can interpret this as an estimate of the value of (3.111} y'[:c+h,t2}=F{I+h,'2.y{.r+h,f2}}, which is the derivative of y at the midpoint of the interval [1,1 4 11]. One step of Euler's method is used to estimate the value of y(;:: | iii/2). As we will see, Heun's method and the midpoint method are both more accurate than Euler's method. In order to analyze one-step methods, it is again helpful to consider the error made by the method in onestep assuming the solution was exact up to that point. Like in our analysis of Euler's method (see Theorem 3.34 above), this leads to the denition of the truncation error 313 1 _ In. (3112} Tn = W _ @(xn,y[xn},h}, which again can be viewed as i "v N times the error made in the {n + ljth step, assuming the solution was exact up to the nth step. Remark 3.43. Some authors prefer to dene the (local) truncation error Without the fac L tor h' Denition 3.44. we say that a one-step method given by (D is consistent if the maximal truncation error max,l |Tn| converges to zero as h } D. A bit of calculus shows that it is not diH'icult to determine when a onestep method is consistent. Fact 3.45 (Consistencyr criterion}. A onestep method given by i is consistent if and onty tf@[3,y.j = F[:t:._ y} on sotntion curves {any} = [I,y[1'}}. Why is this true? This is because as h. } I]. (3.113} w @(I,y(1'),hj converges to (3.114} in) animal (since 0 so that (3.118) max | En) 0. Example 3.50. Euler's method has order 1. Every consistent one-step method has at least order 1, so in some sense Euler's method is worst-possible. Exercise 3.51 (*). Use Taylor expansion to show that Heun's method and the midpoint method both have order 2. Hint: It is enough to show that the truncation error satisfies |T,| bjkj, j= 1 where K1 = F(c, y) k2 = F(x + cgh, y + azhki) ka = F(x + cah, y + azhki + aazhk2) ka = F(I + cgh, y + ashkit .. . + as,a_ 1hks-1) and where the bi, ci, a are real coefficients. We refer to s as the number of stages of the method. A Runge-Kutta method is completely determined by the coefficients bi, Ci, dij.Runge-Kutta methods are often visualized using a Butcher tableau: C2 021 ca 031 032 (3.129) . . Ca . . . 08,8-1 b1 b2 Example 3.55. Euler's method (s = 1), the midpoint method (s = 2) and Heun's method (s = 2) are explicit Runge-Kutta methods that can be represented by Butcher tableaus: (3.130) 0 0 Example 3.56. Method (C) from Exercise 3.52 is a 3-stage Runge-Kutta method represented by the Butcher tableau GOINTOIN O (3.131) 0 As we have seen in Exercise 3.52, this method has order 3. (Whereas the method (D), which is also a consistent 3-stage Runge-Kutta method only has order 1.) Example 3.57. The classical Runge-Kutta method of order 4 (also known as RK4) is the 4-stage Runge-Kutta method given by (3.132) 0 1 The ease of implementation and high accuracy of this method have made it very popular in practical applications. Exercise 3.58. Implement this 4-stage Runge-Kutta method and confirm empirically using the initial value problems from Exercises 3.38 and 3.52 that is has order 4. From Fact 3.45 we see that a Runge-Kutta method is consistent if and only if (3.133) [bj = 1. j=1 Consistent Runge-Kutta methods have at least order 1. Note carefully that the order of an s-stage Runge-Kutta is not necessarily given by s. It is not obvious at all how to chooseK $ 4 S = K 8 11 ? TABLE 1. Minimum number of stages s for an explicit Runge-Kutta method of order k [B]. For large & the answer is not known. the coefficients bi, ci, a; so that the resulting method has a large order. In fact, it is an unsolved problem to determine for a given & what the minimum value of s is so that there exists an s-stage Runge-Kutta method of order k. For k 2 5 it is well-known that s > k is necessary (see Butcher [B]). Determining the order of a given Runge-Kutta method is typically achieved by analyzing the truncation error (3.112) using appropriate Taylor expansion. Exercise 3.59. (i) Write down the Butcher tableau and the map o for a general 2-stage Runge-Kutta method. (ii) Determine all 2-stage Runge-Kutta methods of order at least 1. (iii*) Determine all 2-stage Runge-Kutta methods of order 2. Solution (*). (i) (3.134) a (3.135) D(x, y, h) = bF(x, y) + b2F(x + ch, y + ahF(x,y)). (ii) A method has order at least 1 if and only if it is consistent. Thus we need by + by = 1. The values of c, a are arbitrary. (iii*) We need to study the truncation error (3.136) In = =(y(In th) - y(In)) - D(In, y(In ), h). By Taylor's theorem we have (3.137) y(In th) = y(In) +hy'(In) + jhly"(In) + O(h3). Here O(h ) is the usual "Big-O" notation: it is a placeholder for an "error term" that is bounded by a constant times h3 in absolute value. With this we have (3.138) In = y'(In) + shy"(In) - D(In, y(In), h) + O(h?)1 (Graded). (i) Determine the exact solution to the IVP y' = siWy 11(1) = 2 (ii) Implement the 4stage RungeKutta method given by the Butcher tableau Test it on the IVP from (i J on the interval I [1 6] {so {11 5): Generate a scatter plot of the approximate solution for step size 31 51(plot each point of the approximation, do not connect the points) and also plot the curve of the exact solution (in the same plot). Output the decimal value of the global approximation error IEN = M5) 1m where y 2 at) is the exact solution from (i) and gun; is the approximation that your program produced (note N = afh. = 25). (iii) Conrm that the order of the method is 4 by computing the ap proximation error for various small step sizes say 171 22 for 1' _2. .,Q. then plotting log |E1v against logUi} and estimating the slope of the line (as shown 1n class or differently}. Note: To receive full credit for (ii), (iii) it sufces to include the source code and its output including plots. Include brief comments in your source code to explain what you are doing