% begin preamble
\documentclass[11pt,titlepage]{article}
\usepackage[margin=1in,noheadfoot]{geometry}
\usepackage[parfill]{parskip}
\usepackage{graphicx}
\usepackage[colorlinks,urlcolor=blue]{hyperref}
\usepackage{amssymb}
\usepackage{amsbsy}
\usepackage{wrapfig}

\newcommand{\e}[1]{\times 10^{#1}}
\newcommand{\q}{\mathrm}
\newcommand{\where}[1][where]{\qquad\mbox{#1}\qquad}
\newcommand{\bs}[1]{\boldsymbol{#1}}
\newcommand{\bh}[1]{\hat{\boldsymbol{#1}}}
\newcommand{\bv}[1]{\vec{\boldsymbol{#1}}}
\newcommand{\bt}[1]{\tilde{\boldsymbol{#1}}}
\newcommand{\ad}{^{\dagger}}
\newcommand{\flow}[1]{\mathcal{e}_{#1}}
\newcommand{\flux}[1]{\phi_{#1}}
\newcommand{\dslash}{d \hspace{-0.5ex}\rule[1.2ex]{0.6ex}{.1ex}}
\newcommand{\ts}{\textsuperscript}
\def\rcurs{{\mbox{$\resizebox{.16in}{.08in}{\includegraphics{scriptr}}$}\!\!}}
\def\brcurs{{\mbox{$\resizebox{.16in}{.08in}{\includegraphics{boldr}}$}\!\!}}
\def\hrcurs{{\mbox{$\hat \brcurs$}\!}}
%end preamble

\begin{document}
\begin{center}
  {\bf University of Kentucky, Physics 306 \\
    Homework \#11, Rev. A, Bonus }
\end{center}

\vspace{.3 in}

{\bf 1.~The Relaxation method,} largely due to R.~V.~Southwell, is a
straightforward method of solving PDEs based on the mean value property of the
Laplacian.  It is applicable to a wide variety of problems, not only Laplace's
equation.  We will use it to calculate the potential distribution within a
square with two adjacent edges at 100~V, one at 0~V and one at 50~V.

{\bf \quad a)} Draw a square grid of 36 points: 20 on the edges of the square,
and 16 in the interior.

{\bf \quad b)} Fill in the boundary potentials on all edges, and initial
estimates of the potentials at the interior points on your grid.  This is your
first approximation.

{\bf \quad c)} Correct your first approximation as follows: replace the value
of each interior point with the average of its four adjacent squares.  Explain
why this procedure `relaxes' the potential closer to $\nabla^2 V = 0$ (see
\#2b).

{\bf \quad d)} Repeat part c) until none of the changes are greater than 0.5
V, and sketch the equipotentials for 20~V, 40~V, 60~V, and 80~V.  You may want
to solve this problem on a computer, using Excel, Matlab, or
\href{http://www.gnu.org/software/octave/}{Octave} (open source program
similar to Matlab).

\quad How would one extend this method to solve Neumann (flux) boundary
conditions?


{\bf 2.~Finite Differences Method}.  Re-solve problem \#1 noniteratively as
outlined below:

{\bf \quad a)} Given the discrete approximation of a function $f_i = f(x_i)$
evaluated at regular intervals $x_i=x_0+i\, \Delta x$, show that $f'(x_i)
\approx (f_{i+1}-f_{i})/\Delta x \approx (f_{i}-f_{i-1})/\Delta x$ are two
approximations of the first derivative.  What is their average?

{\bf \quad b)} Derive a formula for the discrete second derivative using the
two first derivatives of part a), and extend it to a discrete 2-dimensional
Laplacian.

{\bf \quad c)} Turn Laplace's equation into a matrix equation by treating the
function $V(x)$ on interior points of the grid in problem \#1 as a
16-component vector $\bs V = [V_{11},V_{12},V_{13},V_{14},V_{21}, \ldots,
V_{44}]$.  Use part b) to represent $\nabla^2V$ as a $16\times 16$ matrix $\bs
K$ multiplied by $\bs V$.  Note that the Laplacian can be approximated by a
matrix because it is a linear operator.  Each row should represent Laplace's
equation evaluated at one interior point.  There is no place in $\bs K \bs V$
for the fixed boundary potentials $V_{0n}$, $V_{5n}$, $V_{m0}$, and $V_{m5}$
for $m,n=1,2,3,4$, so put these values in a fixed vector $\bs B$, such that
the matrix equation $\bs K \bs V + \bs B=\bs 0$ represents the full boundary
value problem.

{\bf \quad d)} Solve $\bs K \bs V = - \bs B$ for the potential $V$ at each
point and compare your result with \#1.  DO NOT do this by hand!  Try using
the online applet
\href{http://www.bluebit.gr/matrix-calculator}{http://www.bluebit.gr/matrix-calculator}
or one of the software packages
\href{http://www.gnu.org/software/octave/}{Octave}, Matlab, Maple, or
Mathematica.

\newpage

{\bf 3.~Finite Element Method}.  In this problem we will investigate one of
the most common methods of solving partial differential equations numerically
(especially involving the Laplacian).  This method is very flexible and can be
used to solve PDEs on irregularly shaped domains (like cars, or electrodes for
your awesome new experiment).  This problem is self-contained, but for extra
details or hints, refer to the article
\href{http://en.wikipedia.org/wiki/Finite_element}
{http://wikipedia.org/wiki/Finite\_element}.

{\bf \quad a)} Show that 
\begin{equation}
\int_R (\nabla^2 u) v d\tau = \oint_{\partial R} (\nabla u) v
\cdot {\bs{da}} - \int_R \nabla u \cdot \nabla v  d\tau.
\end{equation}
Assuming that
$v=0$ on the boundary, this means that the equation $\nabla^2 u=f$ can be
written 
\begin{equation}
  - \int_R \nabla u \cdot \nabla v  d\tau = \int_R f v d\tau,
  \label{integral_eq}
\end{equation} 
which is now a first order integral equation which must hold true for any {\it
  test function} $v({\bs r})$.  This is called the {\it weak form} of
Laplace's equation.

{\bf \quad b)} In a one-dimensional space, Eq.~\ref{integral_eq} can be
discretized by choosing appropriate ``basis functions'' for $v({\bs r})$.  The
$i^{\rm{th}}$ {\it tent function} is defined as $v_i(x) = (1-|x-i|)\;
\theta(1-|x-i|), $ where $\theta(x)= \{ 1$ if $x>0$, and $0$ if $x<0 \}$ is
the Heaviside step function and $i=1,2,3,4$.  Plot each of these functions on
the same graph.

{\bf \quad c)} Sketch the function $f(x)=2 v_1(x) + 4 v_2(x) + 3 v_3(x) + 1
v_4(x)$.  In the same way, any function defined on $0<x<5$ can be approximated
by a linear combination of these basis functions $f(x)\approx \sum_{i=1}^4 f_i
v_i(x)$, where $f_i = f(i)$.

{\bf \quad d)} Convert Eq.~\ref{integral_eq} into four linear equations
($i=1,2,3,4$) by substituting $v(x) \to v_i(x)$ and approximating $u(x)=\sum_j
u_j v_j(x)$ and $f(x)=\sum_j f_j v_j(x)$.  Combine these equations into a
single matrix equation $-Lu=Mf$, where $u=[u_1,u_2,u_3,u_4]^T$,
$v=[v_1,v_2,v_3,v_4]^T$, $M_{ij}=\int_{-\infty}^{\infty} v_i(x) v_j(x) dx$ and
$L_{ij}=\int_{-\infty}^{\infty} \nabla v_i(x) \cdot \nabla v_j(x) dx$.  Note
that $\nabla = d/dx$ in one dimension.  Perform the integrals of each
component of $L$ and $M$ to obtain numerical matrices.
 
{\bf \quad e)} Use part d) to solve the boundary value problem $\nabla^2 u =
3$ on the region $0<x<5$, with boundary conditions $u(0)=0$ and $u(5)=0$, by
solving the above matrix equation for $u_i$.  Solve the boundary value problem
analytically and compare with the finite element result.
 
\quad Describe how this method could apply to higher dimensions.  Hint:
Define the 2-dimensional tent function $v_{ij}(x,y) = v_i(x) \; v_j(y)$, and
show how it can be used as above.  The matrices get very large for 2- or
3-dimensional problems, but there are packages like FlexPDE
(\href{http://www.pdesolutions.com/index.html} {http://www.pdesolutions.com},
free student version) and COMSOL (\href{http://comsol.com}{http://comsol.com})
which do all the bookkeeping, and implement this method for arbitrary
geometries and complex partial differential equations from a GUI interface.

\end{document}

%%% Local Variables:
%%% mode: latex
%%% TeX-master: t
%%% End:
