Math 6610 • Analysis of Numerical Methods I

{Finite-Difference Methods}
Consider the initial boundary value problem
\[
  v_t=\beta v_{xx},\qquad 0<x<1,\quad t>0,\qquad v(x,0)=f(x)
\]
with Dirichlet boundary conditions
\[
  v(0,t)=0,\qquad v(1,t)=0.
\]
Suppose that the interval \([0,1]\) is partitioned by meshpoints \(x_j=jh\) for
\(j=0,1,\ldots,N\) and \(h=1/N\). A possible discretization of the above problem uses the
Crank-Nicolson scheme for the differential equation:
\[
  \left(I-\frac{\beta k}{2}D_+D_-\right)u_j^{n+1}
  =
  \left(I+\frac{\beta k}{2}D_+D_-\right)u_j^n
\]
for \(j=1,\ldots,N-1\), combined with the boundary conditions:
\[
  u_0^n=0\qquad\text{and}\qquad u_N^n=0
\]
for each time level \(n\). Initial data is obtained by setting \(u_j^0=f_j\) for
\(j=0,1,\ldots,N\). (Note that \(D_+\) and \(D_-\) are the usual forward and backward
difference operators.) Set up the matrix system that must be solved in order to
advance the numerical solution from time level \(n\) to time level \(n+1\). Show that the
suggested discretization is stable in the \(\ell_2\) norm
\[
  \lVert u\rVert_2=\left(\sum_{j=0}^N h\lvert u_j\rvert^2\right)^{1/2}
\]
for all choices of timestep \(k\) and spacestep \(h\). (Note that von Neumann analysis
is not enough to prove stability for the initial boundary value problem.) What other
issues are important to consider in deciding whether this scheme is a good choice
for the problem?

\textbf{Fact 1:} A real symmetric \(n\times n\) matrix \(A\) can be diagonalized by an
orthogonal similarity transformation, and \(A\)'s eigenvalues are real.

\textbf{Fact 2:} The \((N-1)\times(N-1)\) matrix \(M\) defined by
\[
  M_{ij}=
  \begin{cases}
    -2, & i=j,\\
    1, & \lvert i-j\rvert=1,\\
    0, & \text{otherwise}
  \end{cases}
  \qquad (1\le i,j\le N-1)
\]
has eigenvalues
\[
  \mu_l=-4\sin^2\left(\frac{\pi l}{2N}\right),\qquad l=1,2,\ldots,N-1.
\]

\textbf{Fact 3:} The \((N+1)\times(N+1)\) matrix \(M\) defined by
\[
  M_{ij}=
  \begin{cases}
    -1, & i=j\in\{0,N\},\\
    -2, & i=j,\ 1\le i\le N-1,\\
    1, & \lvert i-j\rvert=1,\\
    0, & \text{otherwise}
  \end{cases}
  \qquad (0\le i,j\le N)
\]
has eigenvalues
\[
  \mu_l=-4\sin^2\left(\frac{\pi l}{2(N+1)}\right),\qquad l=0,1,\ldots,N.
\]

\textbf{Fact 4:} For a real \(n\times n\) matrix \(A\), the Rayleigh quotient of a vector
\(x\in\mathbb{R}^n\) is the scalar
\[
  r(x)=\frac{x^T Ax}{x^T x}.
\]
The gradient of \(r(x)\) is
\[
  \nabla r(x)=\frac{2}{x^T x}(Ax-r(x)x).
\]
If \(x\) is an eigenvector of \(A\) then \(r(x)\) is the corresponding eigenvalue and
\(\nabla r(x)=0\).