Nine-point finite-difference stencil

A nine-point finite-difference stencil on a square grid couples a central value to its four axial neighbours and four diagonal neighbours. Ordering the grid by columns produces a block tridiagonal matrix whose diagonal and off-diagonal blocks are symmetric tridiagonal Toeplitz matrices.

For the Dirichlet discrete Laplacian on a square, let $T=h^{-2}\mathrm{tridiag}(1,-2,1)$. The standard axial five-point operator is $T\otimes I+I\otimes T$, and the axial-weight $2/3$, corner-weight $1/6$ nine-point operator additionally contains $h^{2}T\otimes T/6$. Its eigenvalues are therefore strictly larger than the five-point values for the same sine modes. Both lie above the corresponding negative continuum Laplacian eigenvalues, so the five-point operator has smaller spectral error. For fixed mode numbers the leading errors are respectively ${\pi }^4{h}^2(k^4+l^4)/12$ and ${\pi }^4{h}^2(k^2+l^2{)}^2/12$.

For the square-grid nine-point finite-difference stencil for the Laplacian,If $\Delta u=0$, both displayed correction terms vanish: the mixed sixth derivatives sum to $\Delta u_{xxyy}=0$, and the pure sixth derivatives sum to the negative of that sum. Expanding two orders further yields $D_{9}u=h^6{u}_{xxxxxxxx}/3024+O(h^8)$ for a smooth harmonic function. Thus the compact harmonic scheme has sixth-order normalized consistency of a numerical method and, with the discrete maximum principle inverse bound, sixth-order nodal convergence on a square with exact smooth Dirichlet boundary data. For example $\mathrm{Re}(x+iy{)}^8$ shows that the sixth-order term need not vanish.

Let ${\Gamma }_9$ have center weight $-10/3$, axial weights $2/3$ and diagonal weights $1/6$, and let $M_h$ have center weight $2/3$ and axial weights $1/12$. Taylor expansion gives $h^{-2}{\Gamma }_{9}u=\Delta u+h^2{\Delta }^{2}u/12+O(h^4)$ and $M_{h}u=u+h^2\Delta u/12+O(h^4)$. On a smooth solution of the constant-coefficient Helmholtz equation, the second-order term becomes $h^2\Delta (\Delta u+\lambda u)/12=0$. Thus the normalized PDE residual is fourth order, while the original unscaled stencil residual is sixth order. This local cancellation alone does not prove uniform global accuracy at or near a resonance.

The Jacobi method has axial neighbour weights $1/5$ and diagonal weights $1/20$. Interior row sums are one, boundary-adjacent sums are smaller, and the positive-weight graph is connected. If an eigenvalue had modulus one, a maximal-modulus eigenvector component would force all its neighbours to have the same modulus. Propagation reaches a deficient row, a contradiction. This proves strict spectral-radius control even though the infinity norm of the iteration matrix can equal one.

Extend the interior grid vector by zero. With axial weights $2/(3h^2)$ and diagonal weights $1/(6h^2)$, the negative central coefficient equals the total neighbour weight. Expanding the edge-square sum gives the displayed quadratic form. Equality forces all axial differences to vanish, and connection to the zero exterior forces the vector to be zero. Simultaneous row and column permutations preserve symmetry and negative definiteness.

For a smooth solution of $\Delta u=f$, the corrected formula ${\Gamma }_{9}u=h^{2}f+(h^2/12){\Gamma }_{5}f$ has unscaled residual $O(h^6)$ and normalized local truncation error $O(h^4)$. The inverse interior matrix has Euclidean operator norm $O(h^{-2})$; on $O(h^{-2})$ grid points this gives unweighted Euclidean nodal error $O(h^3)$. Uniform derivatives up to the boundary are required for a grid-independent remainder bound.

The nine-point finite-difference stencil for the Laplacian has expansion $D_{9}u=\Delta u+h^2{\Delta }^{2}u/12+h^4(u_{xxxxxx}+u_{yyyyyy}+5u_{xxxxyy}+5u_{xxyyyy})/360+O(h^6)$. If $\Delta u=f$, differentiating this equation gives $u_{xxxxyy}+u_{xxyyyy}=f_{xxyy}$ and $u_{xxxxxx}+u_{yyyyyy}={\Delta }^{2}f-3f_{xxyy}$. Hence the corrected compact equationhas sixth-order normalized consistency of a numerical method. Source derivatives can be supplied analytically; numerical approximation of $\Delta f$ must be fourth order and of the fourth derivatives second order. With exact square-grid Dirichlet boundary data and uniformly smooth solution, the discrete maximum principle inverse bound converts this residual into sixth-order nodal convergence. The unknown stencil stays compact even if more distant known source values are used.

On the unit square with exact Dirichlet boundary data, let $D_{9}U=B_{h}f$ and suppose the normalized residual on the exact smooth solution has supremum norm at most $\epsilon$. The error satisfies $D_{9}e=-\tau$ and vanishes on the boundary. The barrier $q=[x(1-x)+y(1-y)]/4$ has $D_{9}q=-1$, $q\ge 0$ on the boundary, and maximum $1/8$. Apply the discrete maximum principle to $e-\epsilon q$ and $-e-\epsilon q$, obtaining $|e|\le \epsilon q$. Thus $\Vert e{\Vert }_{\infty }\le \epsilon /8$, and fourth-order residual gives fourth-order supremum norm nodal convergence. The proof needs a uniform smoothness bound and a compatible boundary discretization.