1-d problem with Dirichlet boundary conditions

As a simple test case, let us consider the solution of Poisson's equation in one dimension. Suppose that

$$\frac{d^2 u(x)}{dx^2} = v(x),$$ (113)

for $x_l\leq x\leq x_h$, subject to the Dirichlet boundary conditions $u(x_l) = u_l$ and $u(x_h) = u_h$.

As a first step, we divide the domain $x_l\leq x\leq x_h$ into equal segments whose vertices are located at the grid-points

$$x_i = x_l + \frac{i\,(x_h-x_l)}{N+1},$$ (114)

for $i=1,N$. The boundaries, $x_l$ and $x_h$, correspond to $i=0$ and $i=N+1$, respectively.

Next, we discretize the differential term $d^2u/dx^2$ on the grid-points. The most straightforward discretization is

$$\frac{d^2u(x_i)}{dx^2} = \frac{u_{i-1} - 2\,u_i+u_{i+1}}{(\delta x)^2} + O(\delta x)^2.$$ (115)

Here, $\delta x = (x_h-x_l)/(N+1)$, and $u_i \equiv u(x_i)$. This type of discretization is termed a second-order, central difference scheme. It is ``second-order'' because the truncation error is $O(\delta x)^2$, as can easily be demonstrated via Taylor expansion. Of course, an $n$th order scheme would have a truncation error which is $O(\delta x)^n$. It is a ``central difference'' scheme because it is symmetric about the central grid-point, $x_i$. Our discretized version of Poisson's equation takes the form

$$u_{i-1} -2\,u_i + u_{i+1} = v_i\,(\delta x)^2,$$ (116)

for $i=1,N$, where $v_i\equiv v(x_i)$. Furthermore, $u_0=u_l$ and $u_{N+1} = u_h$.

It is helpful to regard the above set of discretized equations as a matrix equation. Let ${\bf u} = (u_1,\, u_2,\, \cdots,\, u_N)$ be the vector of the $u$-values, and let

$${\bf w} = [v_1\,(\delta x)^2 - u_l,\, v_2\,(\delta x)^2,\, v... ...\,\cdots,\, v_{N-1}\, (\delta x)^2,\, v_N\,(\delta x)^2 - u_h]$$ (117)

be the vector of the source terms. The discretized equations can be written as:

$${\bf M}\,{\bf u} = {\bf w}.$$ (118)

The matrix ${\bf M}$ takes the form

$${\bf M} = \left(\begin{array}{rrrrrr} -2 & 1 & 0 & 0 & 0 & 0... ...1 & -2 & 1 \\ 0 & 0 & 0 & 0 & 1 & -2 \\ \end{array} \right)$$ (119)

for $N=6$. The generalization to other $N$ values is fairly obvious. Matrix ${\bf M}$ is termed a tridiagonal matrix, since only those elements which lie on the three leading diagonals are non-zero.

The formal solution to Eq. (118) is

$${\bf u} = {\bf M}^{-1}\,{\bf w},$$ (120)

where ${\bf M}^{-1}$ is the inverse matrix to ${\bf M}$. Unfortunately, the most efficient general purpose algorithm for inverting an $N\times N$ matrix--namely, Gauss-Jordan elimination with partial pivoting--requires $O(N^3)$ arithmetic operations. It is fairly clear that this is a disastrous scaling for finite-difference solutions of Poisson's equation. Every time we doubled the resolution (i.e., doubled the number of grid-points) the required cpu time would increase by a factor of about eight. Consequently, adding a second dimension (which effectively requires the number of grid-points to be squared) would be prohibitively expensive in terms of cpu time. Fortunately, there is a well-known trick for inverting an $N\times N$ tridiagonal matrix which only requires $O(N)$ arithmetic operations.

Consider a general $N\times N$ tridiagonal matrix equation ${\bf M}\,{\bf u} = {\bf w}$. Let ${\bf a}$, ${\bf b}$, and ${\bf c}$ be the vectors of the left, center and right diagonal elements of the matrix, respectively Note that $a_1$ and $c_N$ are undefined, and can be conveniently set to zero. Our matrix equation can now be written

$$a_i\,u_{i-1}+ b_i\,u_i + c_i\,u_{i+1} = w_i,$$ (121)

for $i=1,N$. Let us search for a solution of the form

$$u_{i+1} = x_i\,u_i + y_i.$$ (122)

Substitution into Eq. (121) yields

$$a_i\,u_{i-1} + b_i\,u_i+c_i\,(x_i\,u_i+y_i) = w_i,$$ (123)

which can be rearranged to give

$$u_i = - \frac{a_i\,u_{i-1}}{b_i+c_i\,x_i} + \frac{w_i-c_i\,y_i}{b_i+c_i\,x_i}.$$ (124)

However, if Eq. (122) is general then we can write $u_i = x_{i-1}\,u_{i-1} + y_{i-1}$. Comparison with the previous equation yields

$$x_{i-1} = - \frac{a_i}{b_i+c_i\,x_i},$$ (125)

and

$$y_{i-1} = \frac{w_i-c_i\,y_i}{b_i+c_i\,x_i}.$$ (126)

We can now solve our tridiagonal matrix equation in two stages. In the first stage, we scan up the leading diagonal from $i=N$ to $1$ using Eqs. (125) and (126). Thus,

$$x_{N-1} = -\frac{a_N}{b_N},\mbox{\hspace{1cm}}y_{N-1} = \frac{w_N}{b_N},$$ (127)

since $c_N=0$. Furthermore,

$$x_i = - \frac{a_{i+1}}{b_{i+1}+c_{i+1}\,x_{i+1}},\mbox{\hspa... ...i = \frac{w_{i+1}-c_{i+1}\,y_{i+1}}{b_{i+1}+c_{i+1}\,x_{i+1}}$$ (128)

for $i=N-2, 1$. Finally,

$$x_0 = 0,\mbox{\hspace{1cm}}y_0 = \frac{w_1-c_1\,y_1}{b_1+c_1\,x_1},$$ (129)

since $a_1=0$. We have now defined all of the $x_i$ and $y_i$. In the second stage, we scan down the leading diagonal from $i=0$ to $N-1$ using Eq. (122). Thus,

$$u_1 = y_0,$$ (130)

since $x_0=0$, and

$$u_{i+1} = x_i\,u_i + y_i$$ (131)

for $i=1, N-1$. We have now inverted our tridiagonal matrix equation using $O(N)$ arithmetic operations.

Clearly, we can use the above algorithm to invert Eq. (118), with the source terms specified in Eq. (117), and the diagonals of matrix ${\bf M}$ given by $a_i=1$ for $i=2,N$, plus $b_i = -2$ for $i=1,N$, and $c_i = 1$ for $i=1, N-1$.