\documentclass{article}
\textwidth 6.5in
\oddsidemargin 0in
\evensidemargin 0in

\def\dx2i {\frac{1}{\left(\Delta x\right)^2}}
\def\dy2i {\frac{1}{\left(\Delta y\right)^2}}
\def\dr2i {\frac{1}{\left(\Delta r\right)^2}}
\def\dz2i {\frac{1}{\left(\Delta z\right)^2}}

\def\nx {\mathrm{nx}}
\def\ny {\mathrm{ny}}

\def\cc {\mathrm{cc}_{ij}}
\def\ce {\mathrm{ce}_{ij}}
\def\cw {\mathrm{cw}_{ij}}
\def\cn {\mathrm{cn}_{ij}}
\def\cs {\mathrm{cs}_{ij}}
\def\rhs {\mathrm{rhs}_{ij}}

\begin{document}

\title{Laplace2D}
\author{Karen Camarda}
\date{February 21, 2001}
\maketitle

\abstract{}

\section{Purpose}

\section{The equation being solved}
This thorn solves Laplace's equation
\begin{equation}
\nabla^2 u = 0.
\end{equation}
In the two-dimensions. Hopefully I can use this thorn to solve it in both
Cartesian and (axisymmetric) cylindrical coordinates.

We solve this equation numerically using CG...

\section{Derivation of Finite-differenced equations}

\subsection{Cartesian coordinates}

\subsection{Cylindrical coordinates}
Laplace's equation in cylindrical coordinates takes the form
\begin{equation}
\frac{1}{r} \frac{\partial}{\partial r} \left( r \frac{\partial u}{\partial r}
\right) + \frac{1}{r^2} \frac{\partial^2 u}{\partial \theta^2} +
\frac{\partial^2 u}{\partial z^2} = 0,
\end{equation}
which in axisymmetry simplifies to
\begin{equation}
\frac{1}{r} \frac{\partial}{\partial r} \left( r \frac{\partial u}{\partial r}
\right) + \frac{\partial^2 u}{\partial z^2} = 0.
\end{equation}
A finite-differenced form of this equation is
\begin{equation}
\dr2i (u_{i+1,j}-2 u_{ij} + u_{i-1,j}) +
\frac{1}{2 r_i \Delta r} (u_{i+1,j}-u_{i-1,j}) +
\dz2i (u_{i,j+1}-2 u_{ij} + u_{i,j-1}) = 0.
\end{equation}
This equation can be rewritten in the form
\begin{equation}
\cc u_{ij} + \ce u_{i+1,j} + \cw u_{i-1,j} + \cn u_{i,j+1} + \cs u_{i,j-1} =
\rhs,
\end{equation}
where
\begin{eqnarray}
\cc &=& -2 \left(\dr2i + \dz2i \right) \\
\ce &=& \dr2i + \frac{1}{2 r_i \Delta r} \\
\cw &=& \dr2i - \frac{1}{2 r_i \Delta r} \\
\cn &=& \dz2i \\
\cs &=& \dz2i \\
\rhs &=& 0
\end{eqnarray}

\subsubsection{Boundary conditions}
\label{sec:cylbound}
In axisymmetric cylindrical coordinates, it must be the case that
$\left.\frac{\partial u}{\partial r}\right|_{r=0} = 0$. To enforce this
condition, we take $u$ to be parabolic in $r$ at $r=0$. {\em i.e.}, 
ignoring the $z$-dependence, we take $u$ to have the form $u = A + B r^2$.
To convert this requirement into something we can use, we solve the
system of equations
\begin{equation}
\left\{
\begin{array}{l}
u(0) = A \\
u(\Delta r) = A + B (\Delta r)^2 \\
u(2 \Delta r) = A + 4 B (\Delta r)^2
\end{array}
\right.
\end{equation}
giving the numerical boundary condition
\begin{equation}
u(0) = \frac{4 u(\Delta r) - u(2 \Delta r)}{3}
\end{equation}

\section{Numerical Methods}

\section{Appendix: Coding of boundary conditions for CG}

The equation
\begin{equation}
\cc u_{ij} +
\ce u_{i+1,j} +
\cw u_{i-1,j} +
\cn u_{i,j+1} +
\cs u_{i,j-1} = \rhs
\end{equation}
is valid for the interior points of the computational domain. In C these points
are given by the indices
\begin{eqnarray}
1 \le i < \nx-1 \\
1 \le j < \ny-1
\end{eqnarray}
We incorporate the boundary conditions by adjusting the stencil coefficients
next to the boundary in such a way as to eliminate the boundary points from
the computation.

\subsection{Dirichlet}
This section gives the coding for the case of all boundaries having Dirichlet
boundary conditions.

% Case 1
$i=1$, $j=1$:
\begin{eqnarray}
\rhs &=& \rhs - \cw u_{0j} - \cs u_{i0} \\
\cw &=& 0 \\
\cs &=& 0
\end{eqnarray}

% Case 2
\noindent $2 \le i < \nx-2$, $j=1$:
\begin{eqnarray}
\rhs &=& \rhs - \cs u_{i0} \\
\cs &=& 0
\end{eqnarray}

% Case 3
\noindent $i=\nx-2$, $j=1$:
\begin{eqnarray}
\rhs &=& \rhs - \ce u_{\nx-1,j} - \cs u_{i0} \\
\ce &=& 0 \\
\cs &=& 0
\end{eqnarray}

% Case 4
\noindent $i=1$, $2 \le j < \ny-2$:
\begin{eqnarray}
\rhs &=& \rhs - \cw u_{0j} \\
\cw &=& 0
\end{eqnarray}

% Case 5
% All points interior

% Case 6
\noindent $i=\nx-2$, $2 \le j < \ny-2$:
\begin{eqnarray}
\rhs &=& \rhs - \ce u_{\nx-1,j} \\
\ce &=& 0
\end{eqnarray}

% Case 7
\noindent $i=1$, $j=\ny-2$:
\begin{eqnarray}
\rhs &=& \rhs - \cw u_{0j} - \cn u_{i,\ny-1} \\
\cw &=& 0 \\
\cn &=& 0
\end{eqnarray}

% Case 8
\noindent $2 \le i < \nx-2$, $j=\ny-2$:
\begin{eqnarray}
\rhs &=& \rhs - \cn u_{i,\ny-1} \\
\cn &=& 0
\end{eqnarray}

% Case 9
\noindent $i=\nx-2$, $j=\ny-2$:
\begin{eqnarray}
\rhs &=& \rhs - \ce u_{\nx-1,j} - \cn u_{i,\ny-1} \\
\ce &=& 0 \\
\cn &=& 0
\end{eqnarray}

\subsection{Dirichlet for Cylindrical}
This section shows the coding for the boundary conditions when all boundaries
have Dirichlet boundary conditions, but we are using cylindrical coordinates,
so the $r=0$ boundary actually has a Neumann boundary condition, as
described in section \ref{sec:cylbound}.

% Case 1
$i=1$, $j=1$:
\begin{eqnarray}
\rhs &=& \rhs - \cs u_{i0} \\
\cc &=& \cc + \frac{4}{3} \cw \\
\ce &=& \ce - \frac{1}{3} \cw \\
\cw &=& 0 \\
\cs &=& 0
\end{eqnarray}

% Case 2
\noindent $2 \le i < \nx-2$, $j=1$:
\begin{eqnarray}
\rhs &=& \rhs - \cs u_{i0} \\
\cs &=& 0
\end{eqnarray}

% Case 3
\noindent $i=\nx-2$, $j=1$:
\begin{eqnarray}
\rhs &=& \rhs - \ce u_{\nx-1,j} - \cs u_{i0} \\
\ce &=& 0 \\
\cs &=& 0
\end{eqnarray}

% Case 4
\noindent $i=1$, $2 \le j < \ny-2$:
\begin{eqnarray}
\cc &=& \cc + \frac{4}{3} \cw \\
\ce &=& \ce - \frac{1}{3} \cw \\
\cw &=& 0
\end{eqnarray}

% Case 5
% All points interior

% Case 6
\noindent $i=\nx-2$, $2 \le j < \ny-2$:
\begin{eqnarray}
\rhs &=& \rhs - \ce u_{\nx-1,j} \\
\ce &=& 0
\end{eqnarray}

% Case 7
\noindent $i=1$, $j=\ny-2$:
\begin{eqnarray}
\rhs &=& \rhs - \cn u_{i,\ny-1} \\
\cc &=& \cc + \frac{4}{3} \cw \\
\ce &=& \ce - \frac{1}{3} \cw \\
\cw &=& 0 \\
\cn &=& 0
\end{eqnarray}

% Case 8
\noindent $2 \le i < \nx-2$, $j=\ny-2$:
\begin{eqnarray}
\rhs &=& \rhs - \cn u_{i,\ny-1} \\
\cn &=& 0
\end{eqnarray}

% Case 9
\noindent $i=\nx-2$, $j=\ny-2$:
\begin{eqnarray}
\rhs &=& \rhs - \ce u_{\nx-1,j} - \cn u_{i,\ny-1} \\
\ce &=& 0 \\
\cn &=& 0
\end{eqnarray}

\section{Appendix: Bessel Functions}

\section{Appendix: Comparison with Analytic Solutions}

\subsection{Cartesian coordinates}

\subsection{Cylindrical coordinates}
The Laplace equation in axisymmetric cylindrical coordinates is
\begin{equation}
\frac{1}{r} \frac{\partial}{\partial r} \left( r \frac{\partial u}{\partial r}
\right) + \frac{\partial^2 u}{\partial z^2} = 0.
\end{equation}
This equation can be solved analytically using separation of variables. Assume
the solution $u(r,z)$ can be written $u(r,z)=R(r)Z(z)$. Plugging this form
into the origin differential equation gives
\begin{equation}
\frac{1}{rR} \frac{d}{dr} \left( r \frac{dR}{dr} \right) + \frac{1}{Z}
\frac{d^2 Z}{dz^2} = 0.
\end{equation}
In order for this equation to be true for all $(r,z)$, each term must be
equal to a constant. This gives the set of equations
\begin{eqnarray}
&\frac{1}{Z} \frac{d^2Z}{dz^2} = \alpha^2& \label{eqn:zeqn} \\
&\frac{1}{rR} \frac{d}{dr} \left( r \frac{dR}{dr} \right) + \alpha^2 = 0&
\label{eqn:reqn}
\end{eqnarray}
Taking $\alpha$ to be real, equation (\ref{eqn:zeqn}) has the general solution
\begin{equation}
Z(z) = A \cosh(\alpha z) + B \sinh(\alpha z).
\end{equation}
To solve equation (\ref{eqn:reqn}) we first perform the substitution
$r=\rho/\alpha$, valid for $\alpha \ne 0$. This results in the equation
\begin{equation}
\frac{d^2R}{d\rho^2} + \frac{1}{\rho} \frac{dR}{d\rho} + R = 0.
\end{equation}
This equation is a special case of Bessel's equation and has the general
solution
\begin{eqnarray}
R &=& C J_0(\rho) + D N_0(\rho) \\
  &=& C J_0(\alpha r) + D N_0(\alpha r)
\end{eqnarray}
where $J_0$ and $N_0$ are the zero-th order Bessel and Neumann functions,
respectively. Since we will only consider domains that contain the line
$r=0$, and since the Neumann functions are not regular on this line, we will
set $D=0$. So, the general solution to our original differential equation is
\begin{equation}
u(r,z) = J_0(\alpha r) \left[ A \cosh(\alpha z) + B \sinh(\alpha z) \right].
\end{equation}
Actually that solution is not complete general, since it neglects the
possibility of $\alpha=0$. The solution in this case is
\begin{equation}
u(r,z) = \frac{1}{r} \left( A + B z \right),
\end{equation}
but we will neglect this solution since it blows up at $r=0$.

\begin{thebibliography}{10}

\bibitem{Press86}
W. H. Press, B. P. Flannery, S. A. Teukolsky, and W. T. Vetterling.
{\em Numerical Recipes}. Cambridge University Press, Cambridge, England, 1986.

\bibitem{Wyld76}
H. W. Wyld. {\em Mathematical Methods for Physics}. Addison-Wesley Publishing
Company, Inc., Advanced Book Program, Reading, Massachusetts, 1976.

\end{thebibliography}

\end{document}
