Description
1. Consider the following method for solving the heat equation ut = uxx:
U
n+2
i = U
n
i +
2k
h
2
(U
n+1
i−1 − 2U
n+1
i + U
n+1
i+1 ).
(a) Determine the order of accuracy of this method (in both space and time).
(b) Suppose we take k = αh2
for some fixed α > 0 and refine the grid. For what values of
α (if any) will this method be Lax-Richtmyer stable and hence convergent?
Hint: Consider the MOL interpretation and the stability region of the time-discretization
being used.
(c) Is this a useful method?
2. The m-file heat_CN.m solves the heat equation ut = κuxx using the Crank-Nicolson method.
Run this code, and by changing the number of grid points, confirm that it is second-order
accurate. (Observe how the error at some fixed time such as T = 1 behaves as k and h go to
zero with a fixed relation between k and h, such as k = 4h.)
You might want to use the function error_table.m to print out this table and estimate the
order of accuracy, and error_loglog.m to produce a log-log plot of the error vs. h. See
bvp_2.m for an example of how these are used.
3. Modify heat_CN.m to produce a new m-file heat_trbdf2.m that implements the TR-BDF2
method on the same problem. Test it to confirm that it is also second order accurate. Explain how you determined the proper boundary conditions in each stage of this Runge-Kutta
method.
4. Modify heat_CN.m to produce a new m-file heat_FE.m that implements the forward Euler
explicit method on the same problem. Test it to confirm that it is O(h
2
) accurate as h → 0
1
Math 232: Numerical Solution of Differential Equations II Homework 7
provided when k = 24h
2
is used, which is within the stability limit for κ = 0.02. Note how
many more time steps are required than with Crank-Nicolson or TR-BDF2, especially on
finer grids.
5. Modify heat_CN.m to solve the heat equation for −1 < x < 1 with step function initial data
u(x, 0) =
1 if x < 0,
0 if x ≥ 0.
With appropriate Dirichlet boundary conditions, the exact solution is
u(x, t) = 1
2
erfc
x
√
4κt
where erfc is the complementary error function
erfc(x) = 2
√
π
Z ∞
x
e
−z
2
dz.
(a) Test this routine m = 39 and k = 4h. Note that there is an initial rapid transient decay
of the high wave numbers that is not captured well with this size time step.
(b) How small do you need to take the time step to get reasonable results? For a suitably
small time step, explain why you get much better results by using m = 38 than m = 39.
What is the observed order of accuracy as k → 0 when k = αh with α suitably small
and m even?
2
