Math 566 Grady Wright Homework 5 Due April 16, 2015 1. Consider the following third-order accurate LMM yj+1 = yj−1 + h [7fj − 2fj−1 + fj−2 ] , 3 (a) Derive this method using any technique you desire. Hint: The most straightforward approach is to consider the appropriate approximation to the formal solution Z tj+1 y(tj+1 ) = y(tj−1 ) + f (t, y(t)dt tj−1 of the general IVP y 0 (t) = f (t, y). (b) Plot the stability domain for this LMM and include on your plot the stability domain for the third order Adams-Bashforth method (AB3). Compare and contrast the two stability domains in terms of the problems they might be good for solving. 2. Predictor-Corrector (a) Implement a Matlab function for numerically solving the general vector-valued IVP y0 = f (t, y), α ≤ t ≤ β, y(α) = y0 (y ∈ Rn ), using the fourth-order, predictor-corrector method of AB4 (4-step method) and AM4 (3-step method): ∗ yj+1 yj+1 h [55f (tj , yj ) − 59f (tj−1 , yj−1 ) + 37f (tj−2 , yj−2 ) − 9f (tj−3 , yj−3 )] (AB4 predictor) 24 h ∗ = yj + [9f (tj+1 , yj+1 ) + 19f (tj , yj ) − 5f (tj−1 , yj−1 ) + f (tj−2 , yj−2 )] (AM4 corrector) 24 j = 3, 4, . . . , N − 1 . = yj + Your function should be called abm4, and take as input a function that represents the vector-valued function f (t, y), a, b, the initial condition y0 , and the number of steps to take N . The output of your function should be a column vector t containing all of the time-steps, and a matrix y containing the numerical solution of all the components of the system y at each time-step. You will need to obtain the three values y1 , y2 , and y3 to get your abm4 function started. These values should be obtained with a fourth-order accurate method to keep the scheme consistent. The easiest approach is to use the standard fourth-order Runge-Kutta method (see Example 5.13 on p. 126 of the book). Turn in a printed listing and e-mail me your program(s). Note: Your function must be entirely general. This means that I should be able to create my own Matlab function for any f (t, y), the values a, b, and y0 , and use your program abm4 to try and solve the problem. Furthermore, your program should minimize the number of times the function f (t, y) needs to be called. Try to make your code use matrix-vector products to optimize it for speed. (b) Consider the initial value problem y0 = 1 + y y 2 , 1 ≤ t ≤ 3, + t t y(1) = 0 , where the exact solution is y(t) = t tan(log(t)). Use your abm4 method from part (a) to numerically solve this IVP for 1 ≤ t ≤ 3 for the different N values [50,100,200,400]. Compute the error in the numerical solution at t = 3 and make either a table or a graph clearly showing fourth-order convergence of the solution. On a single graph, plot the solution for N = 100 over the interval 1 ≤ t ≤ 3. 3. All 2-stage, second-order accurate, explicit RK methods have the following Butcher diagram 0 c2 a21 γ1 γ2 where b1 + b2 = 1 and c2 b2 = a21 b2 = 12 . (a) Show that all these methods have the same stability domain. Note: In general one can show that all p-stage explicit RK methods of order p have the same stability domain. (b) Plot the stability domain. (c) You will note from the plot in part (b) that it appears the stability domain does not cover any neighborhood of the imaginary axis in the complex ξ-plane. Show that this is indeed true. Hint: One way to proceed is use the result that you should have obtained from part (a), which is that the characteristic equation for the any RK2 scheme applied to y 0 = λy is r = 1+ξ + 21 ξ 2 , where ξ = λk with k the step-size. Letting r = exp(iθ) in this equation solve for the root with ξ(0) = 0 (this is the root closest to the origin in the complex ξ-plane). Show that this root satisfies ξ(θ) = iθ − i 3 1 2i 5 7 6 θ − θ4 + θ + θ + ... 6 8 15 48 Use this result to argue that no neighborhood of the imaginary axis can be included in the stability domain. 2