Alternating-direction implicit method

Remarks

The matrices $\delta_x^2$ and $\delta_y^2$ are symmetric, negative definite and commute, with common eigenvectors $\sin(p \pi j h) \sin(q \pi k h)$. On such an eigenvector, with eigenvalues $-\alpha$ and $-\beta$ ($\alpha, \beta > 0$), one step multiplies by $$g = \frac{1 - \frac{\Delta t}{2} \alpha}{1 + \frac{\Delta t}{2} \alpha} \cdot \frac{1 - \frac{\Delta t}{2} \beta}{1 + \frac{\Delta t}{2} \beta},$$ and each factor has absolute value less than $1$; this gives clause 2. Eliminating $u^*$ shows that the method equals the two-dimensional Crank–Nicolson method plus the term $\frac{\Delta t^2}{4} \delta_x^2 \delta_y^2 (u^{n+1} - u^n)$, which is $O(\Delta t^3)$ for smooth solutions; this gives clause 3. Each tridiagonal system is solved in linear time by the tridiagonal matrix algorithm, so a step costs $O(M^2)$ operations, against a sparse system with $M^2$ unknowns for Crank–Nicolson.

Going deeper

For a student

The method is an operator splitting of $\delta_x^2 + \delta_y^2$, and its analysis uses that the two pieces commute; for variable coefficients or non-rectangular domains this fails and the analysis is harder. The same iteration, used with varying parameters, is an iterative solver for the discrete Poisson equation and for Lyapunov and Sylvester matrix equations. In three dimensions one usually uses the Douglas form of the method instead.