next up previous
Next: An example tridiagonal matrix Up: Poisson's equation Previous: Introduction


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
\begin{displaymath}
\frac{d^2 u(x)}{dx^2} = v(x),
\end{displaymath} (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

\begin{displaymath}
x_i = x_l + \frac{i\,(x_h-x_l)}{N+1},
\end{displaymath} (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

\begin{displaymath}
\frac{d^2u(x_i)}{dx^2} = \frac{u_{i-1} - 2\,u_i+u_{i+1}}{(\delta x)^2} + O(\delta x)^2.
\end{displaymath} (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
\begin{displaymath}
u_{i-1} -2\,u_i + u_{i+1} = v_i\,(\delta x)^2,
\end{displaymath} (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

\begin{displaymath}
{\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]
\end{displaymath} (117)

be the vector of the source terms. The discretized equations can be written as:
\begin{displaymath}
{\bf M}\,{\bf u} = {\bf w}.
\end{displaymath} (118)

The matrix ${\bf M}$ takes the form
\begin{displaymath}
{\bf M} = \left(\begin{array}{rrrrrr}
-2 & 1 & 0 & 0 & 0 & 0...
...1 & -2 & 1 \\
0 & 0 & 0 & 0 & 1 & -2 \\
\end{array} \right)
\end{displaymath} (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

\begin{displaymath}
{\bf u} = {\bf M}^{-1}\,{\bf w},
\end{displaymath} (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

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

for $i=1,N$. Let us search for a solution of the form
\begin{displaymath}
u_{i+1} = x_i\,u_i + y_i.
\end{displaymath} (122)

Substitution into Eq. (121) yields
\begin{displaymath}
a_i\,u_{i-1} + b_i\,u_i+c_i\,(x_i\,u_i+y_i) = w_i,
\end{displaymath} (123)

which can be rearranged to give
\begin{displaymath}
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}.
\end{displaymath} (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
\begin{displaymath}
x_{i-1} = - \frac{a_i}{b_i+c_i\,x_i},
\end{displaymath} (125)

and
\begin{displaymath}
y_{i-1} = \frac{w_i-c_i\,y_i}{b_i+c_i\,x_i}.
\end{displaymath} (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,
\begin{displaymath}
x_{N-1} = -\frac{a_N}{b_N},\mbox{\hspace{1cm}}y_{N-1} = \frac{w_N}{b_N},
\end{displaymath} (127)

since $c_N=0$. Furthermore,
\begin{displaymath}
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}}
\end{displaymath} (128)

for $i=N-2, 1$. Finally,
\begin{displaymath}
x_0 = 0,\mbox{\hspace{1cm}}y_0 = \frac{w_1-c_1\,y_1}{b_1+c_1\,x_1},
\end{displaymath} (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,
\begin{displaymath}
u_1 = y_0,
\end{displaymath} (130)

since $x_0=0$, and
\begin{displaymath}
u_{i+1} = x_i\,u_i + y_i
\end{displaymath} (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$.


next up previous
Next: An example tridiagonal matrix Up: Poisson's equation Previous: Introduction
Richard Fitzpatrick 2006-03-29