E. Finite Difference Method

Crank-Nicolson Summary

To solve \(u_t = u_{xx}\), first discretize \(x\) to get the system of ODEs \begin{align*} \dot U = \delta_x^+\delta_x^- U, \end{align*} where the operators are defined as \begin{align*} \delta_x^+ &= \frac{S-I}{\Delta x},\\ \delta_x^- &= \frac{I-S^{-1}}{\Delta x},\\ \delta_x^+\delta_x^- &= \frac{S-2I+S^{-1}}{(\Delta x)^2}, \end{align*} and \(U\) is a vector function of \(t\):

\begin{align*} U(t) = \begin{pmatrix} U_1(t)\\ U_2(t)\\ U_3(t)\\ \vdots\\ U_{m-1}(t)\\ \end{pmatrix}. \end{align*} In matrix form, the operators are \begin{align*} SU &= \begin{pmatrix} U_2\\ U_3\\ U_4\\ \vdots\\ U_{m}\\ \end{pmatrix} = \begin{pmatrix} 0 & 1 &&&\\ & 0 & 1 &&\\ && 0 & 1 &\\ &&& \ddots \\ &&&& 0\\ \end{pmatrix}U + \begin{pmatrix} 0\\ 0\\ 0\\ \vdots\\ \beta_R\\ \end{pmatrix}, \\ S^{-1}U &= \begin{pmatrix} U_0\\ U_1\\ U_2\\ \vdots\\ U_{m-2}\\ \end{pmatrix} = \begin{pmatrix} 0 &&&&\\ 1 & 0 &&&\\ & 1 & 0 &&\\ &&& \ddots \\ &&& 1 & 0\\ \end{pmatrix}U + \begin{pmatrix} \beta_L\\ 0\\ 0\\ \vdots\\ 0\\ \end{pmatrix}, \end{align*} where \(\beta_L = U_0\) and \(\beta_R = U_m\) are Dirichlet boundary conditions. Now discretize \(t\) and apply the trapezoidal rule \begin{align*} U^{n+1} &= U^{n} + \int_{t_n}^{t_{n+1}} \dot U\,dt\\ &= U^{n} + \int_{t_n}^{t_{n+1}} \delta_x^+\delta_x^- U\,dt\\ &\approx U^{n} + \Delta t\delta_x^+\delta_x^- \left(\frac{U^{n}+U^{n+1}}{2}\right). \end{align*} So the Crank-Nicolson method is \begin{align*} \left(I - \frac{\Delta t}{2}\delta_x^+\delta_x^- \right)U^{n+1} = \left(I + \frac{\Delta t}{2}\delta_x^+\delta_x^- \right)U^{n}. \end{align*} To write in matrix form, expand the operator to get \begin{align*} \left(I - \alpha\left(S-2I+S^{-1}\right) \right)U^{n+1} = \left(I + \alpha\left(S-2I+S^{-1}\right) \right)U^{n}, \end{align*} where \begin{align*} \alpha = \frac{\Delta t}{2(\Delta x)^2}. \end{align*} Rearrange to get \begin{align*} \left(-\alpha S + (1+2\alpha)I -\alpha S^{-1} \right)U^{n+1} = \left(\alpha S + (1-2\alpha)I + \alpha S^{-1} \right)U^{n}. \end{align*} Thus the scheme in matrix form is \begin{align*} & \begin{pmatrix} 1+2\alpha & -\alpha&&&&\\ -\alpha & 1+2\alpha & -\alpha&&&\\ & -\alpha & 1+2\alpha & -\alpha&&\\ &&& \ddots &\\ &&& -\alpha & 1+2\alpha & -\alpha \\ &&&& -\alpha & 1+2\alpha \end{pmatrix}U^{n+1} + \begin{pmatrix} -\alpha\beta_L^{n+1}\\ 0\\ 0\\ \vdots\\ 0\\ -\alpha\beta_R^{n+1}\\ \end{pmatrix} \\ &= \begin{pmatrix} 1-2\alpha & \alpha&&&&\\ \alpha & 1-2\alpha & \alpha&&&\\ & \alpha & 1-2\alpha & \alpha&&\\ &&& \ddots &\\ &&& \alpha & 1-2\alpha & \alpha \\ &&&& \alpha & 1-2\alpha \end{pmatrix}U^{n} + \begin{pmatrix} \alpha\beta_L^{n}\\ 0\\ 0\\ \vdots\\ 0\\ \alpha\beta_R^{n}\\ \end{pmatrix}. \end{align*} Note that the solver can only go forward in time: One can only solve the unknown \(U^{n+1}\) given \(U^{n}\), not the other way around, or otherwise the scheme will become unconditionally unstable.

Neumann Boundary Conditions

When Neumann boundary conditions \(\beta_L^\prime, \beta_R^\prime\) are given instead of Dirichlet, we set \begin{align*} &-\frac{3}{2\Delta x} U_0 + \frac{4}{2\Delta x} U_1 -\frac{1}{2\Delta x} U_2 = \beta_L^\prime, \\ &\frac{3}{2\Delta x} U_m - \frac{4}{2\Delta x} U_{m-1} +\frac{1}{2\Delta x} U_{m-2} = \beta_R^\prime, \end{align*} where the left hand side are 2nd order approximation of the \(u_x\) values at both boundaries, derived from the method of undetermined coefficients. This can be computed by, for example, the Finite Difference Coefficients Calculator. Equivalently, we can write \begin{align*} U_0 &= -\frac{2\Delta x}{3}\beta_L^\prime + \frac{4}{3} U_1 - \frac{1}{3} U_2, \\ U_m &= \frac{2\Delta x}{3}\beta_R^\prime + \frac{4}{3} U_{m-1} - \frac{1}{3} U_{m-2}. \end{align*} Replacing all \(\beta_L\) and \(\beta_R\) in the above linear system by these new expressions of \(U_0\) and \(U_m\), we obtain the following linear system with Neumann BCs: \begin{align*} & \begin{pmatrix} 1+2\alpha-\frac{4}{3}\alpha & -\alpha+\frac{1}{3}\alpha&&&&\\ -\alpha & 1+2\alpha & -\alpha&&&\\ & -\alpha & 1+2\alpha & -\alpha&&\\ &&& \ddots &\\ &&& -\alpha & 1+2\alpha & -\alpha \\ &&&& -\alpha+\frac{1}{3}\alpha & 1+2\alpha-\frac{4}{3}\alpha \end{pmatrix}U^{n+1} + \begin{pmatrix} \frac{2\alpha\Delta x}{3}(\beta_L^\prime)^{n+1}\\ 0\\ 0\\ \vdots\\ 0\\ -\frac{2\alpha\Delta x}{3}(\beta_R^\prime)^{n+1}\\ \end{pmatrix} \\ &= \begin{pmatrix} 1-2\alpha +\frac{4}{3}\alpha & \alpha-\frac{1}{3}\alpha&&&&\\ \alpha & 1-2\alpha & \alpha&&&\\ & \alpha & 1-2\alpha & \alpha&&\\ &&& \ddots &\\ &&& \alpha & 1-2\alpha & \alpha \\ &&&& \alpha-\frac{1}{3}\alpha & 1-2\alpha+\frac{4}{3}\alpha \end{pmatrix}U^{n} + \begin{pmatrix} -\frac{2\alpha\Delta x}{3}(\beta_L^\prime)^{n}\\ 0\\ 0\\ \vdots\\ 0\\ \frac{2\alpha\Delta x}{3}(\beta_R^\prime)^{n}\\ \end{pmatrix}. \end{align*}

Nonlinear Term

When there is a nonlinear term in the equation, say the PDE is \(u_t = Lu + N(u)\), then the corresponding Crank-Nicolson method is \begin{align*} U^{n+1} &= U^{n} + \int_{t_n}^{t_{n+1}} \dot U\,dt\\ &= U^{n} + \int_{t_n}^{t_{n+1}} L_h U + N(U)\,dt\\ &\approx U^{n} + \Delta t\left(\frac{L_h U^{n} + N(U^{n}) + L_h U^{n+1} + N(U^{n+1})}{2}\right). \end{align*} Thus \begin{align*} \left(I - \frac{\Delta t}{2}L_h \right)U^{n+1} - \frac{\Delta t}{2} N(U^{n+1}) = \left(I + \frac{\Delta t}{2}L_h \right)U^{n} + \frac{\Delta t}{2} N(U^{n}). \end{align*} Now on the right hand side apply the 2nd order approximation (in \(t\)) \begin{align*} N(U^{n+1}) = N(U^{n}) + (U^{n+1} - U^{n})\circ N^\prime(U^n), \end{align*} where \(\circ\) denotes the element-wise product. We obtain \begin{align*} &\left(I - \frac{\Delta t}{2}L_h \right)U^{n+1} - \frac{\Delta t}{2} \left(N(U^{n}) + (U^{n+1} - U^{n})\circ N^\prime(U^n)\right) = \left(I + \frac{\Delta t}{2}L_h \right)U^{n} + \frac{\Delta t}{2} N(U^{n})\\ \implies\quad&\left(I - \frac{\Delta t}{2}L_h \right)U^{n+1} - \frac{\Delta t}{2} N^\prime(U^n)\circ U^{n+1} = \left(I + \frac{\Delta t}{2}L_h\right)U^{n} - \frac{\Delta t}{2} N^\prime(U^n)\circ U^{n} + \Delta t N(U^n)\\ \implies\quad&\left(I - \frac{\Delta t}{2}L_h - \frac{\Delta t}{2} N^\prime(U^n)\right)U^{n+1} = \left(I + \frac{\Delta t}{2}L_h - \frac{\Delta t}{2} N^\prime(U^n)\right)U^{n} + \Delta t N(U^n). \end{align*} In most cases this requires reconstruction of the matrix at each time step. Note that the above argument is not specific to finite difference method. Same argument works for spectral methods too.

Alternating-Direction Implicit Method (No Cross Term)

To solve \(u_t = u_{xx} + u_{yy}\), first discretize \(x\) and \(y\) to get the system of ODEs \begin{align*} \dot U = (\delta_x^+\delta_x^- + \delta_y^+\delta_y^-) U, \end{align*} where \(U\) is a matrix function of \(t\):

\begin{align*} U(t) = \begin{pmatrix} U_{11}(t) & U_{12}(t) & \cdots & U_{1m_2}(t)\\ U_{21}(t) & U_{22}(t) & \cdots & U_{2m_2}(t)\\ \vdots & \vdots & \ddots & \vdots\\ U_{m_11}(t) & U_{m_22}(t) & \cdots & U_{m_1m_2}(t)\\ \end{pmatrix}. \end{align*} The Crank-Nicolson method is \begin{align*} \left(I - \frac{\Delta t}{2}(\delta_x^+\delta_x^- + \delta_y^+\delta_y^-) \right)U^{n+1} = \left(I + \frac{\Delta t}{2}(\delta_x^+\delta_x^- + \delta_y^+\delta_y^-)\right)U^{n}, \end{align*} where the operators can be approximated by \begin{align*} \left(I - \frac{\Delta t}{2}(\delta_x^+\delta_x^- + \delta_y^+\delta_y^-) \right) &\approx \left(I - \frac{\Delta t}{2}\delta_x^+\delta_x^- \right)\left(I - \frac{\Delta t}{2}\delta_y^+\delta_y^- \right), \\ \left(I + \frac{\Delta t}{2}(\delta_x^+\delta_x^- + \delta_y^+\delta_y^-)\right) &\approx \left(I + \frac{\Delta t}{2}\delta_x^+\delta_x^- \right)\left(I + \frac{\Delta t}{2}\delta_y^+\delta_y^- \right). \end{align*} These approximations have truncation error smaller than \(O(\Delta t^2 + \Delta x^2 + \Delta y^2)\). The system of equations we need to solve can be written as \begin{align*} \left(I - \frac{\Delta t}{2}\delta_x^+\delta_x^- \right)\left(I - \frac{\Delta t}{2}\delta_y^+\delta_y^- \right)U^{n+1} = \left(I + \frac{\Delta t}{2}\delta_x^+\delta_x^- \right)\left(I + \frac{\Delta t}{2}\delta_y^+\delta_y^- \right)U^{n}, \end{align*} which can be split into \begin{align*} &\left(I - \frac{\Delta t}{2}\delta_x^+\delta_x^- \right)V^{n} = \left(I + \frac{\Delta t}{2}\delta_y^+\delta_y^- \right)U^{n}, \\ &\left(I - \frac{\Delta t}{2}\delta_y^+\delta_y^- \right)U^{n+1} = \left(I + \frac{\Delta t}{2}\delta_x^+\delta_x^- \right)V^{n}, \end{align*} where \(V^n\) is an intermediate value. This is the Peaceman-Rachford scheme. Although not immediately obvious (as the operators are not commutative), it can be shown that thus obtained solution will be equivalent to the above. Now finding \(U^{n+1}\) is a two-step process. Since each step involves operator in only one direction, one can simply apply the Thomas algorithm \(m_j\) times, avoiding iterative solvers.

ADI With Cross Term

To solve \(u_t = u_{xx} + u_{yy} + u_{xy}\), first discretize in space to get the system of ODEs

\[\dot U = (\delta_x^+\delta_x^- + \delta_y^+\delta_y^- + \delta_x^0\delta_y^0) U,\]

where \begin{align*} \delta_x^0 &= \frac{S_x-S_x^{-1}}{2\Delta x}, \\ \delta_y^0 &= \frac{S_y-S_y^{-1}}{2\Delta y}. \end{align*} No way to split the operator into separate ones in different directions, we treat the mixed derivative fully explicitly and write \begin{align*} \left(I - \frac{\Delta t}{2}(\delta_x^+\delta_x^- + \delta_y^+\delta_y^-) \right)U^{n+1} = \left(I + \frac{\Delta t}{2}(\delta_x^+\delta_x^- + \delta_y^+\delta_y^-) + \Delta t\delta_x^0\delta_y^0\right)U^{n}. \end{align*} By splitting the opertor only on the left, we obtain \begin{align*} \left(I - \frac{\Delta t}{2}\delta_x^+\delta_x^- \right)\left(I - \frac{\Delta t}{2}\delta_y^+\delta_y^- \right)U^{n+1} = \left(I + \frac{\Delta t}{2}(\delta_x^+\delta_x^- + \delta_y^+\delta_y^-) + \Delta t\delta_x^0\delta_y^0\right)U^{n}. \end{align*} A variant of this scheme is the Douglas-Rachford scheme \begin{align*} &\left(I - \frac{\Delta t}{2}\delta_x^+\delta_x^- \right)V^n = \left(I + \frac{\Delta t}{2}\delta_x^+\delta_x^- + \Delta t\delta_y^+\delta_y^- + \Delta t\delta_x^0\delta_y^0\right)U^{n}, \\ &\left(I - \frac{\Delta t}{2}\delta_y^+\delta_y^- \right)U^{n+1} = V^n - \frac{\Delta t}{2}\delta_y^+\delta_y^- U^n. \end{align*} Since the cross term is not integrated with the trapezoidal rule, the convergence order in time of the Douglas-Rachford reduces to 1. To recover 2nd order convergence in time, apply a predictor-corrector scheme:

• Predictor: \begin{align*} \left(I - \frac{\Delta t}{2}\delta_x^+\delta_x^- \right)V_1^{n} &= \left(I + \frac{\Delta t}{2}\delta_x^+\delta_x^- + \Delta t\delta_y^+\delta_y^- + \Delta t\delta_x^0\delta_y^0\right)U^{n}, \\ \left(I - \frac{\Delta t}{2}\delta_y^+\delta_y^- \right)V_2^{n} &= V_1^n - \frac{\Delta t}{2}\delta_y^+\delta_y^- U^n. \end{align*}

• Corrector: \begin{align*} &\left(I - \frac{\Delta t}{2}\delta_x^+\delta_x^- \right)V_3^{n} = \left(I + \frac{\Delta t}{2}\delta_x^+\delta_x^- + \Delta t \delta_y^+\delta_y^- + \frac{\Delta t}{2}\delta_x^0\delta_y^0\right)U^{n} + \frac{\Delta t}{2} V_2^{n}, \\ &\left(I - \frac{\Delta t}{2}\delta_y^+\delta_y^- \right)U^{n+1} = V_3^n - \frac{\Delta t}{2}\delta_y^+\delta_y^- U^n. \end{align*}

Here the Douglas-Rachford scheme is first applied and the result \(V_2^{n}\) is a predictor of the solution at time \(t_{n+1}\), to be corrected. This scheme has 2nd order convergence both in time and space.