1. Derivation of linear multistep (LMS) methods Consider the following third-order accurate LMS method

.

• Derive this method using any technique you desire

Hint: The most straightforward approach is to consider the appropriate approximation to the formal solution Z ^tn+3

u(t_n+3) = u(t_n+1) + f(t,u(t))dt

tn+1

of the general IVP u⁰(t) = f(t,u).

• Draw the stencil for this method
• Determine if this method is zero-stable
• Determine if this method is consistent
• Determine if this method converges
• Determine the order of accuracy of this method
• Plot the stability domain for this method 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. Would you ever want to use this method?

2. Heat equation: Crank-Nicolson

(a) Implement a function for numerically solving the solving the 1-D heat equation

.

using the Crank-Nicolson scheme (i.e. second order finite difference for the spatial derivative and trapezoidal rule for the time derivative). Your function should take as input: functions representing the initial condition and boundary conditions, the number of points for the uniform spatial discretization (m + 1), the time span to solve the problem over, and the number of time steps (N). It should return the approximate solution at each time step (a matrix), a vector containing all the time steps, and a vector containing the spatial discretization. If using Matlab, then a possible function declaration is

[u,t,x] = cnhteq(f,g0,g1,tspan,alp,N,m)

The function should use a sparse matrix library or the fast discrete sine transform (DST) to solve the linear system that results from the implicit discretization. Hand-in a copy of your function and e-mail it to me.

If possible, reuse your code from homework 3 for the differentiation matrix of u_xx.

(b) Test your function on 1-D heat equation

,

which has the exact solution

Use a waterfall plot to plot the Crank-Nicolson solution of the problem over the interval 0 t 1 for k = 1/16 and h = 1/16. Illustrate the second order accuracy of the method by computing the relative errors in the solution at t = 1 for h = k = 2ⁿ, n = 4,5,6,7,8. Produce a table or plot clearly showing the second order accuracy.

3. Heat equation: BDF2 Repeat problem 3, but use BDF2 for the time-integrator instead of trapezoidal rule. To bootstrap this method use one initial step with Crank-Nicolson. You do not need to produce a waterfall plot, but you do need to illustrate the second order accuracy. Which method is more accurate?
4. Linear stability analysis To numerically solve the equation

u_t + u_xxx = 0, 0 x 1, t 0,

IC: u(x,0) = g(x), BC: periodic

with periodic boundary conditions, a naive numerical analyst proposes to use the simple scheme of forward Euler in time and centered, second-order differences in space:

.

• Using some general principles about numerical stability of IVP solvers, explain to the numerical analyst that this scheme is unconditionally unstable, and thus should never be used.
• Propose an alternative scheme that can be used to successfully approximate this PDE (i.e. give a scheme that is conditionally or unconditionally stable).

5. Wave equation

• Show that the two-way wave equation

u_tt = c²u_xx, 0 x < 2, t > 0, u(x,0) = f(x), u_t(x,0) = g(x), (1)

can be transformed into the equivalent first order system of two equations, given by

u_t + v_x = 0

(2)

v_t + c²u_x = 0

subject to u(x,0) = f(x), v(x,0) = G(x) = ^R g(x)dx for functions u(x,t) and v(x,t). This is a linear, constant coefficient hyperbolic system of the form

q_t + Aq_x = 0 (3)

where A is a 2 2 matrix and q(x,t) = (u(x,t),v(x,t)).

• Write a function for numerically solving (2) using fourth-order centered finite differences in space and the standard fourth order Runge-Kutta (RK4) method in time. For reference, the fourth-order accurate centered finite difference formula is

.

As in the previous problem, compute the solution on the equispaced grid x_j = jh, j = 0,,m, where h = 2/(m + 1).

2

• Use your function to solve the wave equation with m = 200, wave speed c = 1, and the initial conditions^[1]

.

Solve the problem for 0 t 2 using a time-step of k = 2/400. Produce a waterfall plot of the solution. Also plot the error in the solution at t = 2. You can do this by simply computing the difference between the numerical solution at t = 2 and the initial condition.

• Extra credit: Since the domain is periodic, we can greatly increase the accuracy of the spatial derivative approximation without much worry about stability. For this extra credit problem, replace the fourth-order finite difference approximation with a Fourier spectral method approximation of u⁰(x) and repeat part (b).

[1] The initial condition f(x) is not technically periodic, but it looks periodic to machine precision.