Franz J. Vesely > CompPhys Tutorial > Partial Differential Equations  
 
 





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

< >