< >
 Ch. 5, Sec. 3

## 5.3 Boundary Value Problems: Elliptic DE

Standard problem: two-dimensional potential equation,

$\frac{\textstyle \partial^{2} u}{\textstyle \partial x^{2}} + \frac{\textstyle \partial^{2} u}{\textstyle \partial y^{2}} = -\rho(x,y)$

For general $\rho(x,y)$ this is Poisson's equation; if $\rho \equiv 0$ it is called Laplace's equation.

Assuming $\Delta y = \Delta x \equiv \Delta l$ we have

$\frac{\textstyle 1}{\textstyle (\Delta l)^{2}} \left[ \delta_{i}^{2} u_{i,j} +\delta_{j}^{2} u_{i,j} \right] = -\rho_{i,j}$

or

$\begin{eqnarray} \frac{\textstyle 1}{\textstyle (\Delta l)^{2}} \left[ u_{i+1,j} \right. & - & \left. 2u_{i,j}+ u_{i-1,j} +u_{i,j+1}-2u_{i,j}+u_{i,j-1} \right] = - \rho_{i,j} \\ && \;\;\;\;\;\;(i=1,2,\dots N;\;\; j=1,2,\dots M) \end{eqnarray}$

Construct a vector $v$ of length $N.M$ by linking together the rows of the matrix $\{ u_{i,j}\}$:

$\begin{eqnarray} v_{r} &=& u_{i,j},\;\;\;\;{\rm with}\;\; r=(i-1)M+j \end{eqnarray}$

$\begin{eqnarray} v_{r-M}+v_{r-1}-4v_{r}+v_{r+1}+v_{r+M}= - (\Delta l)^{2} \rho_{r} \end{eqnarray}$

or
$\begin{eqnarray} A \cdot v & = & b \end{eqnarray}$

with $b \equiv - (\Delta l )^{2} \{ \rho_{1}, \dots \rho_{N.M}\}^{T}$ and

$\begin{eqnarray} A & \equiv & \left( \begin{array}{cccccc} -4 & 1 & \dots & 1 & & \\ 1 &-4 & 1 & & \ddots & \\ \vdots &\ddots &\ddots &\ddots & & \\ 1 & & & & & \\ &\ddots & & & & \end{array} \right) \end{eqnarray}$

Treating the boundaries:

Assume $u_{i,j}=u_{i,j}^{0}$ to be given along the sides: $\Longrightarrow$ $v_{1}=u_{1,1}^{0}$ etc.

At the interior points we have

$\begin{eqnarray} -4v_{7}+v_{8}+v_{12}&=&- (\Delta l)^{2} \rho_{2,2}-u_{2,1}^{0}-u_{1,2}^{0} \end{eqnarray}$

etc., yielding a modified system matrix $A$ and vector $b$.
More specifically, the matrix $A$ has the form

$\left( \begin{array}{rrrrrrrrrrr} -4& 1& & |& 1& & & |& & & \\ 1&-4& 1& |& & 1& & |& & & \\ & 1&-4& |& & & 1& |& & & \\ -& -& -& |& -& -& -& |& -& -& -\\ 1& & & |&-4& 1& & |& 1& & \\ & 1& & |& 1&-4& 1& |& & 1& \\ & & 1& |& & 1&-4& |& & & 1\\ -& -& -& |& -& -& -& |& -& -& -\\ & & & |& 1& & & |&-4& 1& \\ & & & |& & 1& & |& 1&-4& 1\\ & & & |& & & 1& |& & 1&-4\\ -& -& -& |& -& -& -& |& -& -& - \end{array} \right)$

Figure: Treatment of Dirichlet-type boundary conditions $u_{i,j}=u_{i,j}^{0}$ in the case of a $5 \times 5$ lattice

The vector $v$ consists of the nine elements $v_{7},v_{8},v_{9},$ $v_{12},v_{13},v_{14},$ $v_{17},$ $v_{18},$ $v_{19}$, and the vector $b$ has components

$\begin{eqnarray} b_{7} &=& - (\Delta l)^{2} \rho_{7}-u_{1,2}^{0} - u_{2,1}^{0} \\ b_{8} &=& - (\Delta l)^{2} \rho_{8}-u_{1,3}^{0} \\ b_{9} &=& - (\Delta l)^{2} \rho_{9}-u_{1,4}^{0} - u_{2,5}^{0} \\ b_{12} &=& - (\Delta l)^{2} \rho_{12}-u_{3,1}^{0} \\ b_{13} &=& - (\Delta l)^{2} \rho_{13} \\ b_{14} &=& - (\Delta l)^{2} \rho_{14}-u_{3,5}^{0} \\ b_{17} &=& - (\Delta l)^{2} \rho_{17}-u_{4,1}^{0} - u_{5,2}^{0} \\ b_{18} &=& - (\Delta l)^{2} \rho_{18}-u_{5,3}^{0} \\ b_{19} &=& - (\Delta l)^{2} \rho_{19}-u_{4,5}^{0} - u_{5,4}^{0} \end{eqnarray}$

Figure: Potential equation on a $5 \times 5$ lattice: at the points $\circ$ the values of the potential $u(x,y)$ are given (Dirichlet boundary conditions)

Neumann boundary conditions of the form

$\left( \frac{\textstyle \partial u}{\textstyle \partial x} \right)_{i,j} = \alpha_{i,j} \; ; \;\; \left( \frac{\textstyle \partial u}{\textstyle \partial y} \right)_{i,j} = \beta_{i,j}$

may be treated using a linear approximation for $u(x,y)$ at the rim. Again considering the previous example, we proceed as follows:
• Add an additional "skin layer"of lattice points points with linearly extrapolated vales of $u$:

$\begin{eqnarray} u_{0,1} &=& u_{2,1}- 2\alpha_{1,1} \Delta l \\ u_{0,2} &=& u_{2,2}- 2\alpha_{1,2} \Delta l \\ \vdots && \\ u_{1,0} &=& u_{1,2}- 2\beta_{1,1} \Delta l \\ \vdots && \end{eqnarray}$

• At the original boundary points, such as $(1,1)$, we have

$u_{2,1}-2u_{1,1}+u_{0,1}+u_{1,2}-2u_{1,1}+u_{1,0}=-\rho_{1,1} (\Delta l)^{2}$

Expressing the external values by the internal ones we have

$u_{2,1}-2u_{1,1}+u_{2,1}+u_{1,2}-2u_{1,1}+u_{1,2}=-\rho_{1,1} (\Delta l)^{2} +2\alpha_{1,1} \Delta l + 2\beta_{1,1} \Delta l$

$\Longrightarrow$ Same form of the discretized Poisson equation as in the interior region but with "effective" charge densities.

Putting things together, we find

$\left( \begin{array}{ccccccccccccccc} -4& 2& & & & |& 2& & & & & |& & &\\ 1& -4& 1& & & |& & 2& & & & |& & &\\ & 1& -4& 1& & |& & & 2& & & |& & &\\ & & 1& -4& 1& |& & & & 2& & |& & &\\ & & & 2& -4& |& & & & & 2& |& & &\\ - & - & - & - & - & |& - & - & - & - & - & |& - & - & -\\ 1& & & & & |& -4& 2& & & & |& 2& &\\ & 1& & & & |& 1& -4& 1& & & |& &2 &\\ & & & & & |& & & & & & |& & &\\ \end{array} \right)$

Treatment of Neumann-type boundary conditions in the case of a $5 \times 5$ lattice

Subsections

### 5.3.1 Relaxation and Multigrid Techniques

Having transformed the given elliptic PDE into a set of linear equations $A \cdot v = b$, we can apply any iterative technique to find $v$.
For example, the Jacobi scheme reads
$v^{n+1} = \left[ I + \frac{1}{4} A \right] \cdot v^{n} + \frac{(\textstyle \Delta l)^{2}}{\textstyle 4} \rho$

The somewhat faster Gauss-Seidel and SOR methods are also easy to apply.

Problem:

The estimated rate of convergence is linked to the largest (by absolute value) eigenvalue of the iteration matrix constructed from $A$.

In the case of the Poisson equation, the Jacobi and Gauss-Seidel iteration matrices have spectral radii that are close to $1$.
$\Longrightarrow$ Slow convergence!

Solution: Let $e \equiv v - v^{n}$ be the residual error after the $n$-th relaxation step. Convergence depends not only on $A$, but also on the properties of the iterated vector $e$. By applying multigrid schemes we may manipulate these properties so as to improve convergence by orders of magnitude.

Let $K \equiv N \cdot M$ be the total number of grid points; the iterated vectors $v^{n}$ and $e$ are then both of dimension $K$. $e \equiv (e_{0}, \dots e_{K-1})$ may be written as a sum of Fourier components, or modes,

$e_{j} = \frac{\textstyle 1}{\textstyle K} \sum \limits_{k=0}^{K-1} E_{k} \exp^{\textstyle -2 \pi ijk/K}$

where

$E_{k} = \sum \limits_{j=0}^{K-1} e_{j} \exp^{\textstyle 2 \pi ijk/K}$

Relaxation is faster for the "oscillatory" high wave number modes of $e$ than for the "smooth", low wave number components.

Multigrid: Double the lattice width: $\Longrightarrow$ Each wave number is automatically increased by a factor of two. $\Longrightarrow$ Faster convergence of the modes!

Basic procedure:

• Do several iterations on the fine grid, to get rid of the oscillatory modes
• Double the grid spacing by considering only every other lattice point (or by a more refined recipe) and relax the long-wavelength modes through a number of steps.
• If necessary, go through a cascade of such coarsening stages
• Reconstruct the fine grid gradually, using some interpolation scheme.

### 5.3.2 ADI Method for the Potential Equation

Alternating Direction Implicit technique: In addition to $v$, construct another long vector $w$ by linking together the columns of the matrix $\{ u_{i,j}\}$:

$\begin{eqnarray} w_{s} &=& u_{i,j} ,\;\;\;{\rm with}\;\;s=(j-1)N+i \end{eqnarray}$

Conversely,
$j=\mbox{int} \left( \frac{\textstyle s-1}{\textstyle N}\right)+1 ; \;\; i=\left[ (s-1) \; \mbox{mod} \; N \right] + 1$

The discretized potential equation is then

$\begin{eqnarray} w_{s+1}-2w_{s}+w_{s-1}+v_{r+1}-2v_{r}+v_{r-1} &=& - (\Delta l)^{2} \rho_{i,j} \end{eqnarray}$

or

$\begin{eqnarray} A_{1} \cdot v + A_{2} \cdot w & = & b \end{eqnarray}$

$A_{1}$ acts on the "rows" of the $u_{i,j}$ lattice, $A_{2}$ affects the "columns" only: