Math428: Problem Set #7 solution guide

1.

Use the trapezoid rule and BDF2 to solve

$$\begin{eqnarray*}{y^{(1)}}' & = & -\frac{30001}{40} y^{(1)} + \frac{9999}{40} \s... ...40} y^{(2)}, \\ y^{(1)}(0) & = & 1, \\ y^{(2)}(0) & = & 1, \\ \end{eqnarray*}$$

for $0 \leq t \leq 1$ with a step size of h=0.1.

If we let

$$A = \left[ \begin{array}{cc} -\frac{30001}{40} & \frac{9999}{... ...rac{9999}{40}\sqrt{3} & - \frac{10003}{40} \end{array}\right],$$

We can think of trapezoid rule as

$$\v u_{k+1} = \left(I-\frac{h}{2}A\right)^{-1}\left(I+\frac{h}{2}A\right) \v u_k.$$

Similarly, we can write BDF2 as

$$\v u_{k+1} = \left(I-\frac{h}{2}A\right)^{-1}\left(I+\frac{h}{2}A\right) \v u_k.$$

Solving the system, we see that trap rule is stable but oscillates like crazy while BDF2 does not. Both are second order.

As an exercize, it is good to check the eigenvalues of A, and if you do, you will find that the system is stiff.

2.

[Based on KC 8.12:4] Consider the multistep method

$$u_{k+1} = -\alpha u_k + (1+\alpha) u_{k-1} + \frac{h}{2} \left[-\alpha f_{k+1} + (4 + 3 \alpha) f_k \right].$$

For what value(s) of $\alpha$ is this method second order? Generate plots of the region of stability for several values of $\alpha$. Experimentally determine values of $\alpha$ so that the method is A-stable.

As usual, if you Taylor expand u_k+1, f_k+1, associate f_k's with u_k''s, and collect the remainder terms at each power of h, you will find that at the first three orders, the value of $\alpha$ does not matter. However, for h³, we find that the only way to make the coefficient zero is if $\alpha = 8/7$.

If you you fiddle around a bit with the region of absolute stability for this method, you find that reason to believe that the method is A-stable for all $\alpha>0$.

3.

[Double credit] Solve the heat equation

$$\begin{eqnarray*}v_t = v_{xx} + x(1-x), \ \ \ \ 0 \leq x \leq 1, \end{eqnarray*}$$

with initial conditions v(x,0) =x and boundary conditions v(0,t) = 0 and v(1,t) = 1 on the time interval $0 \leq t \leq 2$. First, use the explicit first order finite difference scheme with size $h=\frac{1}{10}$ and $k=\frac{1}{10}$, and then $h=\frac{1}{10}$and $k=\frac{1}{400}$. Next, use Crank-Nicholson with the same parameters, and compare the two methods.

For a modest problem like this, you do not need to worry much about the subtleties of the matrix operations. However, if the mesh were finer, many of the techniques you learned in Math426 would come into play for solving the linear system. Also, different implementations of the algorithm will run faster or slower depending upon how do it.

The first choice of parameters for the explicit method (left) will yield unstable oscillations like the following. Notice that the highest frequencies are amplified most, and be sure to look at the scale on the vertical axes! The higher temporal resolution calculation (right) satisfies the stability requirements and we obtain a more reasonable answer.

The implicit method (Crank-Nicholson) fares much better. While the smaller timestep is less accurate, it is nonetheless stable and reasonable for this PDE.