1. Consider the following method for solving the heat equation ut = uxx:
.
(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(h2) accurate as h → 0
1
Math 232: Numerical Solution of Differential Equations II Homework 7
provided when k = 24h2 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
1
u(x,0) =
0 if x < 0, if x ≥ 0. With appropriate Dirichlet boundary conditions, the exact solution is
where erfc is the complementary error function
erfc(
(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 suitablysmall 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?
