Derivation of an exact difference formula

Consider the linear parabolic PDE

$$\frac{\partial u}{\partial t}=L(t,x,D,D^{2})u$$ (1)

where the operator L is linear and $$D=\frac{\partial}{\partial x}$$
The difference formulae involving two adjacent time levels are obtained from the Taylor expansion

$$u(x,t+k) = \left(1+k\frac{\partial}{\partial t} + \frac{1}{2}K^{2}\frac{\partial^{2}}{\partial t^{2}}+ .......\right) u(x,t)$$

$$= exp\left(K\frac{\partial}{\partial t}\right)u(x,t)$$

If we now put $x=mh,$    $t=nk$ and     $u(mh,nk)=u^{n}_{m},$ then

$$u^{n+1}_{m}=exp(KL)u^{n}_{m}$$ (2)

Now an exact formula connecting D and $$\delta_{x},$$ the central difference operator in the x-direction, is

$$D=\frac{2}{h}\sinh^{-1}\frac{\delta x}{2}=\frac{1}{h}\left(\delta... ...}.3!}\delta^{3}_{x}+ \frac{1^{2}.3^{2}}{2^{4}. 5!}\delta^{5}_{x}+ ..... \right)$$ (3)

If equation (3) to eliminate D in terms of $\delta x$ in equation (2), the exact difference replacement is given by

$$u^{n+1}_{m}=exp KL\left(mh,nk,\frac{2}{h}sinh^{-1}\frac{\delta x}{2},\left(\frac{2}{h}sinh^{-1} \frac{\delta x}{2}\right)\right)u^{n}_{m}$$ (4)

All difference formulae in common use for solving equation (1) are approximations of equation (4).

Explicit Formulae
An explicit formula involves only one grid point at the advanced time level

$$t=(n+1)k$$

Consider the heat equation given by

$$\frac{\partial u}{\partial t}=\frac{\partial^{2}u}{\partial x^{2}}$$ (5)

Here $L=D^{2}$ and equation(2) becomes

$$u^{n+1}_{m}=exp (KD^{2})u^{n}_{m}$$ (6)

and from equation (3), we have

$$D^{2}=\frac{1}{h^{2}}\left(\delta^{2}_{x}-\frac{1}{12}\delta^{4}_{ x}+\frac{1}{90}\delta^{6}_{x}......\right)$$ (7)

Substituting this value of $D^{2}$ in equation (6) following by expansion leads to

$$u^{n+1}_{m}=\left[1+r\delta^{2}_{x} + \frac{1}{2}r\left(r-\frac{1}{6}\right)\delta^{4}_{x} + ...... \right]u^{n}_{m }$$ (8)

Where $$r=\frac{k}{h^{2}}$$ is the mesh ratio. From equation (8), if we retain only second order central differences, the forward difference formula

$$U^{n+1}_{m}=(1+r\delta^{2}_{x})U^{n}_{m}$$ (9)

is obtained which, on substitution for $\delta^{2}_{x}$ leads to

$$U^{n+1}_{m}=(1-2r)U^{n}_{m}+r(U^{n}_{m+1}+U^{n}_{m-1})$$ (10)

where $U^{n}_{m}$ is an approximation to $u^{n}_{m}.$
Truncation Error:
Let us investigate the local accuracy of the finite difference formula(10). Introduce the difference between the exact solution of the differential and difference equations at the grid point $(mh, nk)$ as

$$Z^{n}_{m}=u^{n}_{m}-U^{n}_{m}$$ (11)

Using Taylor's theorem

$$u^{n+1}_{m}=u^{n}_{m}+k\left(\frac{\partial u}{\partial t }\right... ...{1}{2}k^{2}\left(\frac{\partial^{2}u}{\partial t^{2}}\right)^{n}_{m}+ .........$$

$$u^{n}_{m\pm1}=u^{n}_{m}\pm h\left(\frac{\partial u}{\partial x}\r... ...}{24}h^{4}\left(\frac{\partial^{4}u}{\partial x^{4}}\right)^{n}_{m}\pm.........$$

and so

$$u^{n+1}_{m}-(1-2r)u^{n}_{m}-r\left(u^{n}_{m+1}+u^{n}_{m-1}\right)... ... ^{2}}-\frac{1}{6r}\frac{\partial^{4}u}{\partial x^{4}}\right)^{n}_{m}+........$$ (12)

>From equation (5), (10), (11) and (12), the result

$$Z^{n+1}_{m}=(1-2r)Z^{n}_{m}+r(Z^{n}_{m+1}+Z^{n}_{m-1})+\frac{1}{2... ...}}-\frac{1}{6 r}\frac{\partial^{4}u}{\partial x^{4}}\right)^{n}_{m}+...........$$

is obtained. The quantity

$$\frac{1}{2}K^{2}\left(\frac{\partial ^{2}u}{\partial t^{2}}-\frac{1}{6r}\frac{\partial^{4}u}{\partial x^{4}}\right)^{n}_{m}+..........$$ (13)

is defined as the local truncation error of formula (10) and the principal part of the truncation error is

$$\frac{1}{2}K^{2}\left(\frac{\partial^{2}u}{\partial t^{2}}-\frac{1}{6r}\frac{\partial^{4}u}{\partial x^{4}}\right)^{n}_{m}$$ (14)

Implicit Formulae:
An implicit formula involves more than one grid point at the advanced time level $t =(n+1)k$. These formulae can often be obtained from equation (2) written in the central form

$$exp\left(-\frac{1}{2}KL\right)u^{n+1}_{m}=exp\left(\frac{1}{2}KL\right)u^{n}_{m}$$ (15)

For the heat equation

$$\frac{\partial u}{\partial t}=\frac{\partial^{2}u}{\partial x^{2}}$$

the equation (15) becomes

$$exp\left(-\frac{1}{2}KD^{2}\right)u^{n+1}_{m}=exp\left(\frac{1}{2}KD^{2}\right)u^{n}_{m}$$ (16)

correct to second differences

$$D^{2}=\frac{1}{h^{2}}\delta^{2}_{x}$$

and substitution in equation (16) followed by expansion leads to the central difference formula

$$\left(1-\frac{1}{2}r\delta^{2}_{x} \right)U^{n+1}_{m}=\left(1+\frac{1}{2}r\delta^{2} x\right)U^{n}_{m}$$ (17)

with a principal truncation error of $O(k^{3}+kh^{2})$. This is the crank-Nicolson formula and may be written in the form

$$(1+r)U^{n+1}_{m}-\frac{1}{2}r\left(U^{n+1}_{m+1}+U^{n+1}_{m-1}\right)=(1-r)U^{n}_{m}+\frac{1}{2}r\left(U^{n}_{m+1}+u^{n}_{m-1}\right)$$ (18)

Solution of Tridiagonal systems:
The implicit difference formula given above has involved three unknown values of U at the advanced time level $t =(n+1)k$. The system of linear algebraic equations arising from the implicit difference formulae which must be solved at each time step is a special case of the tridiagonal system

$$-\alpha_{m}U_{m-1}+\beta_{m}U_{m}-r_{m}U_{m+1}=b_{m}$$

for $1\leq m \leq M-1,$ where $U_{0}$ and $U_{M.}$ are known from the boundary condition. If

$$\alpha_{m}>0, \quad \beta_{m}>0, \quad r_{m}>0$$

and $\beta_{m}\geq (\alpha_{m}+r_{m})\quad 1\leq m\leq M-1$, a highly efficient method is known for solving the tridiagonal system. The method is given as follows:
consider the difference relation

$$U_{m}=W_{m}U_{m+1}+ g_{m}$$

for $0\leq m \leq M-1$, from which it follows that

$$U_{m-1}=W_{m-1}U_{m}+ g_{m-1}$$

If this is used to eliminate $U_{m-1}$ from the original difference formula defining the tridiagonal system, the result

$$U_{m}=\frac{r_{m}}{\beta_{m}-\alpha_{m}W_{m-1}}U_{m+1}+\frac{b_{m}+\alpha_{m} g_{m-1}}{\beta_{m}-\alpha_{m} W_{m-1}}$$

is obtained, and so

$$W_{m}=\frac{r_{m}}{\beta_{m}-\alpha_{m}W_{m-1}},g_{m}=\frac{b_{m}+\alpha_{m}g_{m-1}}{\beta_{m}-\alpha_{m}W_{m-1}}$$

If $U_{0}=0$, then $W_{0}=g_{0}=0$, in order that the difference relation

$$U_{0}=W_{0}U_{1}+g_{0}$$

holds for any $U_{1}$. The remaining $W_{m}, g_{m}$ $(m=1,2,---,M-1)$ can now be computed as

$$W_{1}=\frac{r_{1}}{\beta_{1}}\quad \qquad g_{1}=\frac{b_{1}}{\beta_{1}}$$

$$W_{2}=\frac{r_{2}}{\beta_{2}-\alpha_{2}W_{1}}\quad \qquad g_{2}=\frac{b_{2}+\alpha_{2}g_{1}}{\beta_{2}-\alpha_{2}W_{1}}$$

.
.
.
.
.
.
.

$$W_{M-1}=\frac{r_{M-1}}{\beta_{M-1}-\alpha_{M-1}W_{M-2}}\quad \qquad g_{M-1}=\frac{b_{M-1}+\alpha_{M-1} g_{M-1}}{\beta_{M-1}-\alpha_{M-1}W_{M-2}}$$

If $U_{M}=Y,$ then $U_{1},U_{2},......., U_{M-1}$ are computed as

$$U_{M-1}=W_{M-1} Y + g_{M-1}$$

$$U_{M-2}=W_{M-2}U_{M-1}+g_{M-2}$$

.
.
.
.
.
.
.

$$U_{1}=W_{1}U_{2} + g_{1}$$

In using this method, substantial errors will appear in the computed values of $U_{1},U_{2},........,U_{M-1}$ unless

$$\vert W_{m}\vert\leq 1\qquad(m=1,2---,M-1)$$

Now

$$W_{1}=\frac{r_{1}}{\beta_{1}}\leq 1$$

$$W_{2}=\frac{r_{2}}{\beta_{2}-\alpha_{2}W_{1}}\leq \frac{r_{2}}{\beta_{2}-\alpha_{2}}\leq 1$$

and so on, since $\alpha_{ m}>0,\quad \beta_{m}>0,\quad r_{m}>0,$ $\beta_{m}\geq(\alpha_{m}+r_{m}) \quad for \quad 1\leq m\leq M-1$. This leads to $0<W_{m}\leq 1 \quad (1\leq m\leq M-1)$.
Convergence:
The problem of convergence of a finite difference method for solving equation (1) consists of finding the condition under which

$$u(X,T)-U(X,T),$$

The difference between the exact solutions of the differential and difference equations at a fixed point $(X,T)$, tends to zero uniformly, as the net is refined in such a uniformly, as the net is refined in such a way that $h, k\longrightarrow 0$ and $m,n \longrightarrow\infty$, with $mh(=X)$ and $nk(=T)$ remaining fixed. The fixed point $(X,T)$ is anywhere within the region $R$ under consideration, and it is sometimes convenient in the convergence analysis to assume that $h,k$ do not tend to zero independently but according to zone relationship like

$$k=rh^{2}$$ (19)

where $r$ is a constant.
As an example of a convergence analysis for difference formula (10), we introduce

$$Z^{n}_{m}=u^{n}_{m}-U^{n}_{m}$$

the difference between the theoretical (exact) solutions of the differential and difference equations at the grid point X=mh, T=nk. From equation (12), this satisfies the equation

$$Z^{n+1}_{m}=(1-2r)Z^{n}_{m}+r\left(Z^{n}_{m+1}+Z^{n}_{m-1}\right)+O\left(k^{2}+kh^{2}\right)$$ (20)

$$\quad 0<r\leq \frac{1}{2}$$ (21)

the coefficients on the right hand side of equation (20) are all non-negative and so

$$\vert Z^{n}_{m+1}\vert\leq(1-2r)\vert Z^{n}_{m}\vert+r\vert Z^{n}_{m+1}\vert+r\vert Z^{n}_{m-1}\vert+A(K^{2}+Kh^{2})$$

$$\leq Z^{(n)}+A(K^{2}+Kh^{2})$$

where A depends on the upper bounds for $$\frac{\partial^{ 2}u}{\partial t^{2}}$$ and $$\frac{\partial^{4}u}{\partial x^{4}}$$ and $Z^{(n)}$is the maximum modulus value of $Z^{n}_{m}$ over the required range of $m$. Thus

$$Z^{(n+1)}\leq Z^{(n)}+A(K^{2}+Kh^{2})$$

and so if $Z^{(0)}=0$(the same initial data for differential and difference equations),

$$Z^{(n)}\leq nA(K^{2}+Kh^{2})$$

$$=TA(K+h^{2})$$

$\rightarrow 0$ as $h, k\rightarrow 0$ for fixed $X,T$. This establishes convergence if the expression (21) is satisfied.

Stability:
The problem of stability of a finite difference scheme for solving equation (1) consists of finding conditions under which

$$U^{n}_{m}-U^{rn}_{m}\left(\equiv Z^{n}_{m}\right)$$

the difference between the theoretical and numerical solutions of the difference equation, remains bounded as $n$ increases, K remaining fixed for all $m$. There are two methods which are commonly used for examining stability of a finite difference scheme.

The Von Neumann Method:
In this method, a harmonic decomposition is made of the error Z at grid points, at a given time level, leading to the error function.

$$E(x)=\sum_{j} A _{j} e ^{i \beta_{j} x}$$

where in general the frequencies $\vert\beta_{ j}\vert$ and $j$ are arbitrary. It is necessary to consider only the single term $e^{i\beta x}$ where $\beta$ is any real number. For convenience, suppose that the time level being considered corresponds to t=0. To investigate the error propagation as t increases, it is necessary to find a solution of the finite difference equation which reduces to $e^{i\beta x}$ when t=0. Let such a solution be

$$e^{\alpha t} e^{i\beta x}$$

where $\alpha=\alpha(\beta)$ is, in general, complex. The original error component $e^{i\beta x}$ will not grow with time if
$\vert e^{\alpha k}\vert\leq 1$ for all $\alpha$. This is Von Neumann criterion for stability. As an example, let us examine the stability of finite difference scheme (10). Since $Z^{n}_{m}$ satisfies the original difference equation, we get

$$Z^{n+1}_{m}=(1-2r)Z^{n}_{m}+r\left(Z^{n}_{m+1}+Z^{n}_{m-1}\right)$$ (22)

Let $Z^{n}_{m}=e^{\alpha n k} e ^{i \beta m h }= \xi^{n}e ^{i \beta m h}$ , where $\xi={e^\alpha k}$. Then equation (22) gives

$$\xi^{n+1}e^{i\beta m h}=\xi^{n}\left[(1-2r)e^{i\beta m h}+r\left(e^{i\beta(m+1)}e^{h}e^{i\beta(m-1)h}\right)\right]$$

Cancelling $\xi^{ n } e^{i \beta m h}$ on both sides leads to

$$\xi = (1-2r) + r (e ^{i \beta h }+e ^{-i \beta h })$$

$$= 1-2r(1 - cos \beta h)$$

$$= 1-4 r sin^{2} \frac{\beta h}{2}$$

The quantity $\xi$ is called the amplification factor. For stability, $\vert\xi\vert\leq 1$, for all values of $\beta h$, and so

$$-1\leq 1-4 r sin^{2} \frac{\beta h}{2}\leq 1\quad (\forall\beta h)$$

The right hand side of the inequality is satisfied if $r>0$ and the left hand side gives

$$r\leq \frac{1}{2 sin^{ 2 }\frac{\beta h}{2}}$$

leading to the stability condition $0<r\leq\frac{1}{2}$.

The Matrix Method:
If $Mh=1$, the totality of difference equation connecting values of $V$ at two neighboring time levels can be written in the matrix form

$$AU^{n+1}=BU^{n}$$ (23)

where $U^{s}\quad(s=n,n+1)$ denotes the column vector

$$\left[U^{s}_{1},U^{s}_{2},......U^{s}_{M-1}\right]^{T}$$

and A,B are square matrices of order $(M-1)$. If the difference formula is explicit A=I. Now equation (23) can be written in the explicit form

$$U^{n+1}=C U^{n}$$

where $C=A^{-1}B$ provided $\vert A\vert\neq 0$. The error vector

$$U^{n}-U^{rn}(\equiv Z^{n})$$

satisfies

$$Z^{n+1}=CZ^{n}$$

from which it follow that

$$\vert\vert Z^{n+1}\vert\vert\leq \vert\vert C \vert\vert\qquad \vert\vert Z^{n}\vert\vert$$

Where $\vert\vert.\vert\vert$ denotes a suitable norm. The necessary and sufficient condition for the stability of a finite difference scheme based on a constant time step and proceeding indefinitely in time is $\vert\vert C\vert\vert\leq 1$. When $C$ is symmetric, $\vert\vert C\vert\vert _{2}=\max \limits_{s} \vert\lambda_{s}\vert$where $\lambda _{s} (s=1,2,.....,M-1)$ are the eigen values of $C$, and $\vert\vert.\vert\vert _{2 }$ denotes the $L_{2}$ norm. As an example of the matrix method for examining stability , we consider the finite difference scheme (10). Here we have,

$$B=C=\left(% \begin{array}{ccccccccc} (1-2r) & r & 0 & 0 & 0 ... ...0 & 0 & 0 & 0 & 0 & 0 & 0 & r & (1-2r) \\ \end{array}% \right)$$


The eigen values of this matrix are $$\lambda _{s}=1-4r sin^{2}\frac{s \pi}{2N},\quad s=1,2....,\quad N-1$$ and thus the method is stable if

$$-1\leq 1-4r sin^{2}\frac{s\pi}{2N }\leq 1\qquad s=1,2,.....,N-1$$

which leads to

$$0<r\leq\frac{1}{2},$$

an identical condition obtained by the method of Von Neumann.
A difference approximation to a parabolic equation is consistent, if truncation error $k\rightarrow 0$ as $h, k\rightarrow 0$.

Hyperbolic Equation in one space variable: The simplest hyperbolic problem is that of the vibrating string

$$\frac{\partial^{2}u}{\partial t^{2}}=\frac{\partial^{2}u}{\partial x^{2}}$$ (24)

in the domain $R=[0\leq x\leq 1]\times[t>0],$ satisfying the following initial conditions

$$\left. \begin{aligned}u(x,0)=f_{1}(x) \\ \frac{\partial u}{\partial t}(x,0)=f_{2}(x) \end{aligned} \right\}\qquad 0\leq x\leq 1$$

and boundary conditions

$$\left. \begin{aligned}u(0,t)=g_{1}(t)\\ u(1,t)=g_{2}(t) \end{aligned} \right\}\qquad t>0$$

we place a mesh of points $(x_{m},t_{n})$ on R, where $$x_{m}=m h\qquad m=0,1,2,......,M,\quad Mh=1$$

$$t_{n} = nk\qquad n=0,1,2,......$$

The exact difference replacement of (24) at the nodal points $(x_{m},t_{n})$ is given by

$$\left(sinh ^{-1}\frac{\delta t}{2}\right)^{2} u(x_{m},t_{n})= p^{2} \left(sinh ^{-1} \frac{\delta x}{2}\right)^{2}u(x_{m},t_{ n})$$ (27)

where $p=k/h$ is mesh ratio and

$$4\left(sinh^{-1} \frac{\delta}{2}\right)^{2}=\delta^{ 2}-\frac{1}{12}\delta^{4}+\frac{1}{90}\delta^{ 6}-.......$$

The explicit and implicit difference scheme for (24) will be obtained by approximating equation (27).

Explicit Difference Schemes:
An explicit difference scheme for (24) is given by

$$\delta^{2}_{t}U^{n}_{m}=p^{2}\delta^{ 2}_{x}U^{n}_{m}$$

which may be written in the form

$$U^{n+1}_{m}=2(1-p^{2})U^{n}_{m}+ p^{2}(U^{n}_{m-1}+U^{n}_{m+1})-U^{n-1}_{m}$$ (28)

where $U^{n}_{m}$ is the approximation to $u^{n}_{m} \equiv u(x_{m},t_{n})$.

If each term in (28) is expanded in Taylor's series about the nodal point $(x_{m},t_{n})$ and the function $u(x_{m},t_{ n})$ satisfies (24), then we fit the function error

$$T^{n}_{m}=K^{2}h^{2}\left[\frac{1}{12}(p^{2}-1)\frac{\partial^{4}... ...^{2}(p^{4}-1)+ \frac{\partial^{6}u(x_{m},t_{n})}{\partial x^{6}}+.......\right]$$

For p=1, the truncation error vanishes and so the difference representation of (24) is obtained as

$$U^{n+1}_{m}=U^{n}_{m-1}+U^{n}_{m+1}-U^{n-1}_{m}$$

In order to start computation, we require data on the two lines $t=0$ and $t=k$. The first condition in (25) gives $U^{0}_{m}$ on the initial line as

$$U^{0}_{m}=f_{1}(mh)\qquad 0\leq m\leq M$$

We can use the second condition in (25) to find values on the line $t=k$. Using the central difference approximation for the derivative, i.e

$$\frac{\partial u^{0}_{m}}{\partial t} \simeq \frac{U^{1}_{m}-U^{-1}_{m}}{2k}$$

In the second condition in (25) and eliminating $U^{-1}_{m}$ from (28) for n=0, we get the formula to give the values on the first time level.

$$U^{1}_{m}=(1-p^{2})f_{1}(mh)+Kf_{2}(mh)+\frac{1}{2}p^{2}[f_{1}(m-1)h+f_{1}(m+1)h]\quad 1\leq m\leq M-1$$

The boundary conditions (26) become $U^{n}_{0}=g^{n}_{1} \quad$ and $U^{n}_{M}=g^{n}_{2};\quad n=1,2,.......$

Formula (27) may now be used to advance computation for $n\geq1$.

To examine the finite difference formula (27) for stability, we replace $U^{n}_{m}$ by $\xi^{ n } e^{i \beta m h}$ and get

$$\xi+\frac{1}{\xi}=2-4p^{2}sin^{2}\left(\frac{\beta h}{2}\right)$$

or

$$\xi^{2}-2A\xi+1=0$$ (29)

where $A=1-2p^{2}sin^{2}\left(\frac{\beta h}{2}\right)$
The solution of (29) is given by

$$\xi_{1}=A+\sqrt{A^{2}-1}$$   and$$\quad \xi_{2}=A-\sqrt{A^{2}-1}$$

This gives that $$\xi_{1}=\frac{1}{\xi_{2}}$$.

$$\left . \begin{array}{ccc} If & A>1 & \vert\xi_{1}\vert>1 \\... ...1 & \vert\xi_{1}\vert=\vert\xi_{2}\vert=1 \end{array}% \right .$$

Thus for stability $-1\leq A\leq 1$ or $-1\leq 1-2 p^{2}sin^{2}(\frac{\beta h}{2})\leq1$, which gives $p\leq 1$.
This is the condition of stability.

Implicit Difference Scheme:
In order to improve stability, we now consider an implicit difference replacement of equation (25). This takes the form

$$p^{2}\left\{ a \left[U^{n+1}_{m+1}-2U^{n+1}_{m}+U^{n+1}_{m-1}\rig... ...{n}_{m-1}\right]+a\left[U^{n-1}_{m+1}-2U^{n-1}_{m}+U^{n-1}_{m-1}\right]\right\}$$

$$=\left[U^{n+1}_{m} -2U^{n}_{m}+U^{n-1}_{m}\right]\qquad \qquad \qquad a>0$$

When $$p=\frac{K}{h}$$, one can examine this scheme for stability and find that the scheme is stable for all $p>0$ if $a\geq\frac{1}{4}$. Thus, if $a=\frac{1}{4}$, the implicit difference formula

$$\left[U^{n+1}_{m+1}-2U^{n+1}_{m}+U^{n+1}_{m-1}\right]+2\left[U^{n... ...}_{m}+U^{n}_{m-1}\right]+ \left[U^{n-1}_{m+1}-2U^{n-1}_{m}+U^{n-1}_{m-1}\right]$$

$$=\frac{4}{p^{2}}\left[U^{n+1}_{m}-2U^{n}_{m}+U^{n-1}_{m}\right]$$

The consistency of a finite difference approximation to a hyperbolic equation can be defined briefly. A difference scheme to a hyperbolic equation is consistent; if $$\frac{Truncation}{k^{2}}\rightarrow 0$$ as $h, k\rightarrow 0$.


Elliptic equations in two dimensions:
Suppose that R is a bounded region in the $(x_{1},x_{2})$ plane with boundary $\partial R$. The equation

$$a(x_{1},x_{2})\frac{\partial^{2}x}{\partial x_{1}^{2}}+2b(x_{1},x... ...2},u,\frac{\partial u}{\partial x_{1}},\frac{\partial u}{\partial x_{2}}\right)$$ (30)

is said to be elliptic in R if $b^{2}-ac < 0$ for all points $(x_{1},x_{2})$ in R. Three distinct problems involving equation (30) arise depending on the boundary conditions prescribed on $\partial R:$

1. The first boundary value problem or Dirichlet problem, requires a solution u of equation (30) which takes on prescribed values

   $$u=f(x_{1},x_{2})$$ (31)

   
   
   on the boundary $\partial R.$
2. The second boundary value problem, or Neumann problem, where

   $$\frac{\partial u}{\partial x}=g(x_{1},x_{2}))$$ (32)

   
   
   on the boundary $\partial R$. Here $$\frac{\partial}{\partial x}$$ refers to derivative along the normal to $\partial R$ directed array from the interior of $R$.
3. The third boundary value problem, or Robbins problem, with

   $$\alpha(x_{1},x_{2})u+\beta(x_{1},x_{2})\frac{\partial u}{\partial x}=\gamma(x_{1},x_{2})$$ (33)

   
   
   on $\partial R,$ where $\alpha (x_{1},x_{2})>0,\beta(x_{1},x_{2})>0$ for $(x_{1},x_{2})\epsilon \partial R.$

Before developing finite difference methods of solving elliptic equations, a most useful analytical tool in the study of elliptic partial differential equation will be introduced. This is the Maximum Principle which will be stated for the linear elliptic equation.

$$a\frac{\partial^{2}u}{\partial x_{1}^{2}}+ 2b\frac{\partial^{2}u}... ...}_{2}}+ d\frac{\partial u}{\partial x_{1}}+e\frac{\partial u}{\partial x_{2}}=0$$

where a,b,c,d and e are functions of the independent variables $x_{1},x_{2}$. It is clear that in this case any constant represents a solution of the equation. The maximum principle states that the constants are the only solutions which can assume a maximum or minimum value in the interior of the bounded region R. Alternatively, it states that every solution of the elliptic equation achieves its maximum and minimum values on the boundary $\partial R$ of R.

Laplace's equation in a square:
We consider Laplace's equation

$$\frac{\partial^{2}u}{\partial x^{2}_{1}}+\frac{\partial^{2}u}{\partial x^{2}_{2}}=0$$ (34)

subject to $u=f(x_{1},x_{2})$ on the boundary of the unit square $0\leq x_{1},x_{2}\leq1$. The square region is covered by a grid with sides parallel to the coordinate axes and the grid spacing is $h$. If $Mh=1$, the number of internal grid points or nodes is $(M-1)^{2}$. The coordinates of a typical internal grid point are $x_{1}=lh$, $X_{2}=mh$ ( l and m integers ) and the value of $u$ at this grid point is denoted by $u_{l,m}$. Using Taylor's Theorem, we obtain

$$u_{l+1,m}=(u+h\frac{\partial u}{\partial x_{1}}+ \frac{1}{2}h^{2}... ...{1}^{3}}+\frac{1}{24}h^{4}\frac{\partial^{4}u}{\partial x_{1}^{4}}+.....)_{l,m}$$

and

$$u_{l-1,m}=(u-h\frac{\partial u}{\partial x_{1}}+ \frac{1}{2}h^{2}... ...}^{3}}+ \frac{1}{24}h^{4}\frac{\partial^{4}u}{\partial x_{1}^{4}}+......)_{l,m}$$

and after addition

$$u_{l+1,m}+u_{l-1,m}-2u_{l,m}=(h^{2}\frac{\partial^{2}u}{\partial ... ..._{1}}+ \frac{1}{12}h^{4}\frac{\partial^{4}u}{\partial x_{1}^{4}}+.......)_{l,m}$$

Similarly.

$$u_{l,m+1}+u_{l,m-1}-2u_{l,m}=(h^{2}\frac{\partial^{2}u}{\partial ... ..._{2}}+ \frac{1}{12}h^{4}\frac{\partial^{4}u}{\partial x_{2}^{4}}+.......)_{l,m}$$

and so

$$u_{l+1,m}+u_{l-1,m}+u_{l,m+1}+u_{l,m-1}-4u_{l,m}= \left[ h^{2}\le... ..._{1}^{4}}+\frac{\partial^{ 4}u}{\partial x_{2}^{ 4}}\right)+......\right]_{l,m}$$

leading to the five point finite difference scheme

$$U_{l+1,m}+U_{l-1,m}+U_{l,m+1}+U_{l,m-1}-4U_{l,m}=0$$ (35)

for the Laplace's equation, with a local truncation error

$$\frac{1}{12}h^{4}\left(\frac{\partial^{4}u}{\partial x^{4}_{1}}+ \frac{\partial^{4}u}{\partial x^{4}_{2}}\right)_{l,m}+.......$$

where $U_{l,m}$ denotes the function satisfying the difference equation at the grid point

$$X_{1}=lh,X_{2}=mh$$

The totality of equations at the $(M-1)^{2}$ interval grid points of the unit square leads to the matrix equation

$$AU=K$$ (36)

where $A$ is a matrix of order $(M-1)^{2}$ given by

$$\left(% \begin{array}{ccccccccc} B & -J & 0 & 0 & 0 & 0 & 0 ... ...\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & -J & B \\ \end{array}% \right)$$

with $J$ the unit matrix of order $(M-1)$ and $B$ a matrix of order $(M-1)$ given by

$$\left(% \begin{array}{ccccccccc} 4 & -1 & 0 & 0 & 0 & 0 & 0 ... ...\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & -1 & 4 \\ \end{array}% \right)$$

The vector $U$ and $K$ are given by

$$\{U_{1,1},....,U_{1,M-1};U_{2,1},....,U_{2,M-1};........;U_{M-1,1},.....,U_{M-1,M-1}\}^{T}$$

and

$$K=\{K_{1,1},....,k_{1,M-1};k_{2,1},....,k_{2,M-1};......;k_{M-1,1},.....k_{M-1,M-1}\}^{T}$$

respectively, where $\{\}^{T}$ denotes the transpose. The elements of the vector U constitute the $(M-1)^{2}$ unknowns $U_{l,m}(1\leq l,m,M-1)$ and the elements of the vector K depend on the boundary values $f(x_{1},x_{2})$ at the grid points on the perimeter of the unit square. Because of the large number of zero element in the matrix A, iterative methods are often used to solve the system (36).