HW05 - Boise State University

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