The 5 and 9-Point Finite Difference Scheme

We limit ourselves to the 2-D Poisson equation

$$\nabla^2 U=f\quad (x,y)\in \Omega$$

$$\nabla^2 = \frac{\partial ^2}{\partial x^2}+\frac{\partial ^2}{\partial y^2}$$

$$f=f(x,y) \mbox {is a known function}$$

and domain $\Omega \in {\Bbb R}^2$ is bounded, open, connected and has a piecewise smooth boundary $\partial \Omega$.

Take an arbitrary domain and grid using $\Delta x$ spacing in BOTH $x$ and $y$ direction. The domain grid is $\equiv \Omega_{\Delta x}$. Gridding will aligned parallel to the $x,y$ coordinate system. The grid is depicted in Figure 32

Figure 32: Picture of $\Omega$ with $\Omega _{\Delta x}$ indicated. Solid circles indicate ``interior'' points, hollow circle are ``near-boundary'' points, filled squares indicate ``boundary'' points.

We need to specify boundary conditions. Consider here Dirichlet boundary condition, for concreteness:

$$U(x,y)=\phi(x,y)$$   for $$\qquad (x,y)\in \partial \Omega$$

Take $U^{k,l}_{\Delta}\equiv U(x_0+k\Delta x, y_0+l\Delta x)$, the grid function. We approximate $\nabla^2$ operator acting on the field by center differences to ${\cal O}(\Delta x^2):$ the $5$-point scheme is thus defined as

$$\frac{1}{\Delta x^2}\left(\delta^2_x+\delta^2_y\right)u_{k,l}=f_{k,l}$$ (159)

where $u_{k,l}=U^{k,l}_{\Delta}+{\cal O}(\Delta x^2)$, an approximation to the grid function.

$$f_{k,l}\equiv f(x_0+k\Delta x, y_0+l\Delta x)$$

(160) is written as

$$u_{k-1,l}+u_{k+1, l}+u_{k,1-l}+u_{k, l+1} -4u_{k,l}=(\Delta x)^2 f_{k,l}$$ (160)

computational cell, molecule, or stencil

Of course $\Omega _{\Delta x}$ is comprised of interior, boundary, and near-boundary points. (161) approximates $U_{\Delta x}$ on the interior points. No finite difference approximation is needed for the boundary points. The remaining points, the near boundary points require a special approach since the computational stencil in (161) is not universally applicable. We'll defer discussion of the near-boundary point issue till later. Suppose all values of $u_{k,l}$ are either members of the set of interior or boundary points. Boundary values are known. Interior points are unknown and each is a linear combination defined by (161), i.e. by its nearest neighbors:

(161) can be written as the linear algebraic system of equations

$$A {\bf u}={\bf b} \qquad\qquad \mbox {(exercise, show this)}$$ (161)

${\bf b}$ includes $(\Delta x)^2 f_{k,l}$ and contributions from boundary values.

So we ask some basic questions about the resulting linear algebraic system:

(1)

Is the linear system nonsingular? $\Rightarrow$ finite difference solution ${\bf u}\equiv (u_{k,l})_{k,l, interior}$ integers exists and is unique, if so.

(2)

Suppose ${\bf u}$ exists and is unique. Does ${\bf u}\rightarrow {\bf U}_{\Delta}$ as $\nabla x\rightarrow 0$, and what is the magnitude of the error?

(3)

What efficient methods should be used to solve the linear system? This is crucial since likely to be a very large number of equations.

Take $u_{k,l}$ and arrange it into a 1-D vector of size $m^2$, say. Note that construction of (162) is not uniquely structured: there are $(m^2)!$ ways to arrange it.

Lemma

The matrix $A$ in (161) is symmetric and the set of its eigenvalues is

$$\sigma(A)=\left\{ \lambda_{\alpha, \beta} : \alpha, \beta = 1, 2, \ldots m \right\}$$

where

$$\lambda_{\alpha, \beta}=-4\left\{\sin^2\left[\frac{\alpha\pi}{2(m... ...1)}\right]\right\} \qquad \mbox{ where} \qquad \alpha, \beta = 1, 2, \cdots m.$$

Proof: One can prove symmetry by examinig the elements of the matrix in a general case. For us it would be more instructive to consider an example and from it see the symmetry: take $\nabla^2 U=0$ and use stencil with $m=3$. We'll ignore specific boundaries. Take a domain that's a rectangle and let $\Delta x = 1$. The 5-point approximation is then

$$u_{k-1, l}+u_{k+1, l}+u_{k,l-1}+u_{k, l+1}-4 u_{k, l}=0.$$

Label elements ``lexicographically'', as in the Figure 33

Figure 33: Unit rectangle, gridded with $\Delta x = 1$ and $m=4$. Points labeled lexicographically.

Let $w_i\equiv u(P_i)$ so that $w_i=u_{k,l}$ and $P_i \equiv (x_k, y_l)$ where $k=1, 2, \ldots n$ and $l=1, 2\ldots m$. Finally, set $i=k+(m-l)n$ so that the lexicographic label $i$ is consistent with the position label $k,l$. Then, at each node the equations are:

$$P_1:\quad 4w_1-w_2-w_4 = u_{0,3}+u_{1,4}$$

$$P_2:\quad 4w_2-w_3-w_1-w_5 = u_{2,4}$$

$$P_3:\quad 4w_3-w_2-w_6- = u_{4,3}+u_{3,4}$$

$$P_4:\quad 4w_4-w_5-w_1-w_2 = u_{0, 2}$$

$$P_5:\quad 4w_5-w_6-w_4-w_2 -w_8= 0$$

$$P_6:\quad 4w_6-w_3-w_3-w_9 = u_{4,2}$$

$$P_7:\quad 4w_7-w_8-w_4= u_{0,1}+u_{1,0}$$

$$P_8:\quad 4w_8-w_9-w_7-w_5=u_{2,0}$$

$$P_9:\quad 4w_9-w_8-w_6=u_{3,0}+u_{4,1}$$

Why go through this trouble?? This generates a BANDED MATRIX Of bandwidth = $2n$ !! $\therefore$ MUST USE THIS ORDERING FOR LARGE SYSTEMS. As an exercise, construct the matrix problem using any other arrangement and compare the computational characteristics of the resulting matrix against the one we just worked out. So the matrix $A$ is

$$\left\{ \begin{array}{lllllllll} $4$ & $-1$ & $0$ &$-1$\... ...$0$ & $0$ &$-1$ & $0$ & $-1$ & $4$\\ \end{array}\right\}$$

and the resulting system to solve is

$$A{\bf w}={\bf b}$$

and ${\bf b}$ contains $f$ contributions as well as boundary contributions. So we see it's symmetric. One can show this is a general characteristic of the lexicographic arrangement for this PDE and this computational stencil.

$\Box$

Eigenvalues (general case): eigenvalues of $A$ are independent of how $A$ is formed $\rightarrow$ symmetric perturbations conserve eigenvalues.

The eigenvalues problem, in terms of the original values, is

$$u_{k-1, l}+u_{k+1, l}+u_{k,l-1}+u_{k,l+1}-4u_{k,l}=\lambda u_{k,l}\quad k,l=1,2 \ldots m$$ (162)

where $\lambda$ is an eigenvalue. This is a homogeneous equation, which is further constrained to have eigenfunctions such that $u_{k,0}=u_{k, m+1}=u_{0,1}=u_{m+1,l}=0$ by the Dirichlet boundary conditions.

Given $\alpha, \beta \in \left\{1, 2, \ldots m\right\}$ we have eigenfunctions

$$u_{k,l}=\sin\left(\frac{k\alpha \pi}{m+1}\right)\sin \left(\frac{l\beta\pi}{m+1}\right)\quad k,l=0,1\ldots m+1$$

which automatically satisfy boundary conditions.

Why this form? Check PDE references on harmonic functions and their connection to the equation $\nabla^2U=f$.

Substituting $u_{k,l}$ into (163) and exploiting identity $\sin(\phi-\psi)+\sin (\phi+\psi)=2\sin\phi$ we obtain $\lambda=\lambda_{\alpha, \beta}$

Corollary: The matrix $A$ is negative definite and, a fortiori, nonsingular.

Proof: Already established that $A$ is symmetric. Previous lemma showed that eigenvalues are negative and distinct $\Rightarrow$ nonsingular.

(Recall that all eigenvalues of a symmetric matrix are real, all eigenvalues of a skew-symmetric matrix are purely imaginary, all eigenvalues of a general real matrix are either real or form complex-conjugate pairs. Also, if all eigenvalues of symmetric matrix $> 0\Rightarrow$ matrix is positive definite. If all eigenvalues of symmetric matrix $< 0\Rightarrow$ matrix is negative definite.)

$\Box$

Remarks:

• $U=\sin(\alpha\pi x)\sin(\beta\pi y)$ is the general solution of $\nabla U=\lambda U$ on $\Omega =$ a unit square in $x$ and $y$, with $U=0$ on $\partial \Omega$. $\alpha, \beta$ are positive integers and $\lambda=-(\alpha^2+\beta^2)\pi$. If we sample $U$ on a grid of points $$\left\{\frac{k}{m+1}, \frac{l}{m+1}\right\}$$ for $\alpha, \beta=1, 2, \cdots m$, we obtain the five-point discretization formula for $\nabla^2$ on $\Omega$ = unit square with $(\cal O)$ boundary condition on $\partial \Omega$ and is finite dimensional, of course.

• $(\Delta x)^{-2}\lambda_{\alpha,\beta}$ is a good approximation of $-(\alpha^2+\beta^2)$ is provided $\alpha, \beta$ small compared $m$: expanding $\sin^2\theta$ in a series and using $(m+1)\Delta x=1$
  

  $$\frac{\lambda _{\alpha, \beta}}{(\Delta x)^2}=-4\left(\left\{\lef... ... \right\}+\left\{\left[\frac{\beta \pi}{2(m+1)}\right]^4+\cdots \right\}\right)$$

  $$=-(\alpha^2+\beta^2)\pi+\frac{1}{12}(\alpha^4+\beta^4)\pi^4(\Delta x)^2 + {\cal O}((\Delta x)^4)$$

  
  

Theorem (Approximation error):

$$\mbox {let } e_{k,l}=U_{\Delta x}-u_{k,l} \quad k,l=0,1,\ldots m+1$$

$$\mbox {let } {\bf e } \mbox { denote the vector after lexicographic rearrangement}$$

$$\mbox {in which approx is posed in terms of } {\bf u}.$$

Subject to $f$ being sufficiently smooth, $\exists c>0$ a number independent of $\Delta x$ such that

$$\vert\vert{\bf e}\vert\vert\le c(\Delta x)^2 \quad \Delta x\rightarrow 0$$

where $\vert\vert\cdot\vert\vert$ is the Euclidean norm $\left(i.e. \vert\vert x\vert\vert=[\langle {\bf x}, {\bf x}\rangle]^{\frac12}\right)$ or $l_2$ norm.

Proof: Homework exercise (hint: since $A$ is symmetric, the $l_2$ norm = spectral radius).

$\Box$

Near-Boundary Points:

Previous analysis works on rectangular domains, $L$- shaped domains, etc, provided ratios of all sides are rational numbers. In general, however, we expect near-boundary points in which the $5$-point formula cannot be implemented. Without loss of generality suppose we seek $\nabla^2$ approximation at point $P$ in Figure 34

Figure 34: Graphical construction of data at near-boundary point.

Let's ignore $y$-dependence for now. Given some $z(x)$, approximate $z''$ at $P\sim x_0$ as linear combination of the values of $z$ at $P, Q\sim x_0-\Delta_x, T\sim x_0+\tau\Delta x$. Expanding $z$ about $x_0$ in Taylor series,

$$\frac{1}{(\Delta x)^2}\left[\frac{2}{\tau+1}z(x_0-\Delta x)-\frac{2}{\tau}z(x_0)+\frac{2}{\tau(\tau+1)}x(x_0+\tau\Delta x)\right]$$

$$=z''(x_0)+\frac12(\tau-1)z'' (x_0)\Delta x+{\cal O}((\Delta x)^2)$$

unless $\tau=1$, where everything reduces to central differences, error is just ${\cal O}(\Delta x)$. To get ${\cal O}(\Delta x^2)$ error, we add the function value at $V\approx x_0-2\Delta x$ to linear combination, so that now

$$z''(x_0)=\frac{1}{(\Delta x)^2}\left[\frac{\tau-1}{\tau+2}z(x_0-2... ...a x)-\frac{2(2-\tau)}{\tau+1} z(x_0-\Delta x) -\frac{3-\tau}{\tau}z(x_0)\right.$$

$$+\left.\frac{6}{\tau(t+1)(\tau+2)}z(x_0+\Delta x)\right]+{\cal O}(\Delta x)^2).$$

A good approximation to $\Delta^2 U$ at $P\therefore$ involves 6 points: $P, Q, R, S, T, V$. Assuming $P=(k^0,l^0)$ in our coordinate system, we get the linear combination

$$\left[ \begin{array}{l} \frac{\tau-1}{\tau+2}u_{k^0-2,l^0} +... ...\circ}}=\Delta x^{2}f_{k^{\circ,l^{\circ}}} \end{array}\right]$$

where $u_{k^{\circ}+\tau, l^{\circ}}$ is the value of $U$ at $T$, given by boundary conditions.

Note: if $\tau=1$, we get 5-point formula and $P$ is an internal point, as it should be.

A similar treatment applies in the $y$-direction $\ldots$ note that $\Delta x$ should be small for ${\cal O}(\Delta x^2)$ to be small.