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