< >
 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}$

The potential equation then reads
$\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:

Figure: ADI method

Both $A_{1}$ and A_{2}$are tridiagonal, and not pentadiagonal as the original matrix$A$.$\Longrightarrow$Recursion method  ADI method:$ \begin{eqnarray} (A_{1}+ \omega I) \cdot v^{n+1/2} &=& b-(A_{2} \cdot w^{n} - \omega v^{n}) \\ && \\ w^{n+1/2} & = & U \cdot v^{n+1/2} \\ && \\ (A_{2}+ \omega I) \cdot w^{n+1} & = & b-(A_{1} \cdot v^{n+1/2} - \omega w^{n+1/2}) \end{eqnarray} $The optimal value of the relaxation parameter is given by$\omega = \sqrt{\lambda_{1} \lambda_{2}} ,$where$\lambda_{1}$and$\lambda_{2}$are the smallest and largest eigenvalue, respectively, of the matrix$A$. In the specific case of the potential equation, assuming a lattice with$M=N$, we have$\omega \approx \pi/N$. EXERCISE: Apply the ADI method to the Laplace problem with$M=N=5$. ### 5.3.3 Fourier Transform Method (FT) Let$u_{0,l}=u_{M,l}$and$u_{k,0}=u_{k,N}$(periodic boundary conditions.) Then we may write$ \begin{eqnarray} u_{k,l} &=& \frac{1}{MN} \sum \limits_{m=0}^{M-1}\sum \limits_{n=0}^{N-1} U_{m,n} e^{\textstyle -2\pi i km/M} e^{\textstyle -2\pi i nl/N} \end{eqnarray} $with$ \begin{eqnarray} U_{m,n} &=& \sum \limits_{k=0}^{M-1}\sum \limits_{l=0}^{N-1} u_{k,l} e^{\textstyle 2\pi i km/M} e^{\textstyle 2\pi i nl/N} \end{eqnarray} $A similar expansion is used for the charge density$\rho_{k,l}$:$ R_{m,n} = \sum \limits_{k=0}^{M-1}\sum \limits_{l=0}^{N-1} \rho_{k,l} e^{\textstyle 2\pi i km/M} e^{\textstyle 2\pi i nl/N} $Inserting these expressions in$ \begin{eqnarray} u_{k+1,l} -2u_{k,l}+ u_{k-1,l}+u_{k,l+1}-2u_{k,l}+u_{k,l-1} &=& - (\Delta l)^{2} \rho_{k,l} \end{eqnarray} $we find$ U_{m,n} = \frac{\textstyle -R_{m,n} (\Delta l)^{2}}{\textstyle 2 [\cos 2\pi m/M + \cos 2\pi n/N -2]} $FT method for periodic boundary conditions: • Determine$R_{m,n}$by Fourier analysis of$\rho_{k,l}$• Compute$U_{m,n}$from$R_{m,n}$according to the above equation • Insert$U_{m,n}$in the Fourier series for$u_{k,l}$Use Fast Fourier Transform! ($N \ln N$operations instead of$N^{2}$.) Variants of the method cover other than periodic boundary conditions. ### 5.3.4 Cyclic Reduction (CR) Again, consider$ u_{k+1,l}-2u_{k,l}+ u_{k-1,l} +u_{k,l+1}-2u_{k,l}+u_{k,l-1} = - \rho_{k,l} (\Delta l)^{2} $Let the number of columns in the lattice be$M=2^{p}$, and define the column vectors$ u_{k} \equiv \{ u_{k,l} ; \; l=0, \dots N-1 \}^{T} ;\;\;\;k=0,\dots,M-1 $Then,$ \begin{eqnarray} u_{k-1} & + & T \cdot u_{k} + u_{k+1} = - \rho_{k} (\Delta l)^{2} \end{eqnarray} $where$ \begin{eqnarray} T& \equiv & \left( \begin{array}{ccccc} -2 & 1 & 0 & \dots & \\ 1 &-2 & 1 & & \\ &\ddots &\ddots &\ddots & \end{array} \right) -\left( \begin{array}{ccccc} 2 & 0 & 0 & \dots & \\ 0 &2 & 0 & & \\ &\ddots &\ddots &\ddots & \end{array} \right) \\ & = & \;\;\;\;\;\; B \;\;\;\;\;\; - \;\;\;\;\;\; 2 I \end{eqnarray} B$and$T$are tridiagonal. Now form linear combinations$[k-1]-T \cdot [k] + [k+1]$:$ u_{k-2} + T^{(1)} \cdot u_{k} + u_{k+2} = - \rho_{k}^{(1)} (\Delta l)^{2} $with$ \begin{eqnarray} T^{(1)} & \equiv & 2 I - T^{2} \\ \rho_{k}^{(1)} & \equiv & \rho_{k-1} - T \cdot \rho_{k} + \rho_{k+1} \end{eqnarray} $The "reduced" equation has the same form as before, except that only every other vector$u_{k}$appears in it.$\Longrightarrow$Iterate until$ u_{0} + T^{(p)} \cdot u_{M/2} + u_{M} = - \rho_{M/2}^{(p)} (\Delta l)^{2} $Now the vectors$u_{0}$and$u_{M}$are just the given boundary values$u_{0,l}$and$u_{M,l}(l=0,1, \dots N-1)$. The matrix$T^{(p)}$is known since it arose from the$p$-fold iteration of the above combination rule. While$T^{(p)}$is not tridiagonal any more, it may at least be represented by a$2^{p}$-fold product of tridiagonal matrices:$T^{(p)}=-\prod \limits_{\textstyle l=1}^{\textstyle 2^{\textstyle p}} \left[ T-\beta_{l}I \right] $with$\beta_{l}= 2 \cos\left[ \frac{\textstyle 2(l-1)\pi}{\textstyle 2^{\textstyle p+1}}\right]\Longrightarrow$Solve for the vector$u_{M/2}$by inverting$2^{p}$tridiagonal systems of equations. Now go back to ever more intermediate columns: compute$u_{M/4}$and$u_{3M/4}$follow from$ \begin{eqnarray} u_{0} + T^{(p-1)} \cdot u_{M/4} + u_{M/2} & = & - \rho_{M/4}^{(p-1)} (\Delta l)^{2} \\ && \\ u_{M/2} + T^{(p-1)} \cdot u_{3M/4} + u_{M} & = & - \rho_{3M/4}^{p-1} (\Delta l)^{2} \end{eqnarray} $and so on. Lit. Hockney 81: Combination of the CR technique and the Fourier transform method is superior to most other techniques for solving the potential equation.$\Longrightarrow$"FACR" method (Fourier analysis and cyclic reduction): In place of the column vector$u_{k} \equiv \{ u_{k,l};l=0, \dots N-1\}^{T}$use its$N$Fourier components,$ U_{k}(n) \equiv \sum \limits_{l=0}^{N-1} u_{k,l}\, e^{\textstyle 2\pi i nl/N} ; \;\;(n=0, \dots N-1) \$

Details: see Hockney 81, or Vesely 01.

vesely 2005-10-10

 < >