% SIMIODE-TeacherStudent-Template.tex and simiode.cls for LaTeX2e
% Version 22 November 2016.
%
% This template is to be used with the  class file "simiode.cls".
%
% The latest version of this class and template files
% can be found at
%
%       <http://www.SIMIODE.org> in the Author section.
%
% A "Quick Start Guide" appears below.
%
% This file must be processed with version LaTeX2e or higher.
% Appreciation to Mark Frantz, IUPUI; John Thoo, Yuba College; and Eric Sullivan, Carroll College.
%
\documentclass{simiode}
%
% With this version of the tex file and the associated SIMIODE.cls file you can produce BOTH a Student Version and a Teacher Version of your modeling scenario or technique narrative.
%
% Student Version contains the title STUDENT VERSION along with the NAME OF SCENARIO, with no abstract, no keywords, no tags, no author identifiers, and no COMMENTS section. The Student Version contains only the STATEMENT of the modeling scenario only with appropriate figures, data, instructions, and perhaps a bibliography appropriate for student resounding, but nothing in the bibliography which would "give away" the nature of the modeling scenario, etc.
%
%The Teacher Version contains the material in the Student Version, namely the title NAME OF SCENARIO, now headed with   TEACHER VERSION, but also full author  identification, abstract, key words, tag words, and, following the STATEMENT, a section of COMMENTS which can be discussion of presentation  ideas, "solutions," techniques, etc.
%
% Authors can prepare a full Teacher Version using the notions immediately below  to create a pdf for review, but in review please eliminate all information about the author as we use a double-blind referee system. Such information will be restored when accepted and published in SIMIODE at www.simiode.org.
%
%***********   Toggle the following line to turn the teacher version on (1) or off (0). ***********
%
\def\TeacherEdition{0} % Toggle this line's value.
\newcommand{\teacher}[1]{\ifnum\TeacherEdition=1 #1 \fi}
%
% The \teacher{ } command wraps any text that will appear only in the teacher version of
% the document. See this file for examples. This has been done for you in the document below.
% Rename the pdf after you compile each final version.
%
\usepackage{amsmath,amssymb}
\usepackage{color}
\usepackage{graphicx}
\usepackage{hyperref}
\usepackage{url}
\newtheorem{examp}{Example}
\newenvironment{example}{\begin{examp}\rm}{\end{examp}}
\newtheorem{prob}{Problem}
\newenvironment{problem}{\begin{prob}\rm}{\end{prob}}
%
\renewcommand{\baselinestretch}{1.3}
%
% NOTE regarding macros %%
%
% Any macros defined for your paper should be contained in
% the top matter. Likewise, any environment definitions such
% as \newtheorem or \newenvironment should also go in the
% top matter.
%
% Do not \input your macros as separate files within your
% final version because more files make it harder for the
% editors and users to keep track of and/or modify and
% append your paper.
%
% Top matter is above and the beginning of the document below.
%

\begin{document}

\newpage
%
% You only need to change the NAME OF SCENARIO in the next line and delete the footnote.
% See above for instructions on toggling from Teacher Version to Student Version.
%
\title{\ifnum\TeacherEdition=1 TEACHER VERSION
%\footnote{DELETE THIS FOOTNOTE BEFORE SUBMISSION.  In the \TeX \, source code there is a toggle/option to produce a Teacher Version or a Student Version by ``carving'' out with  $\backslash$teacher\{\} that material which should only appear in the Teacher version. See the source code.}
\else STUDENT VERSION \fi\\Unearthing the Truth}

\markboth{Unearthing the Truth}{Unearthing the Truth}

\author{Kurt Bryan\\
Rose-Hulman Institute of Technology\\
Terre Haute IN 47840 USA\\
kurt.bryan@rose-hulman.edu
}

\makeStitlePDFLaTex


\begin{abstract}
This project introduces \textit{electrical resistivity tomography}, a technique of interest for geophysical imaging, used to produce images of underground features or structures by using electrical current. Specifically, a known electrical current is injected into an object (for example, the earth) and the resulting induced voltages are measured on the surface. From this we attempt to determine the interior electrical resistivity of the object. This is an example of an \textit{inverse problem}.

In this project we first model how electrical current flows through a conductive medium, which results in a classic boundary value problem in two or three dimensions. In the simplest case of interest this is just Laplace's equation. With this model in hand we then consider the inverse problem of recovering a nonconductive internal circular ``void'' in an otherwise uniformly resistive object, as a simple example of how this imaging methodology works.

To tackle this project students would benefit from exposure to basic Vector Calculus and the Divergence Theorem, at a level typical in a third semester calculus course. Some exposure to boundary value problems or Laplace's equation is also helpful, with allied topics like separation of variables and Fourier series solutions, as well as familiarity with the notions of an electric field and electrical potential that one might see in an introductory physics or vector calculus course.
\end{abstract}


\section*{SCENARIO DESCRIPTION}

Prior to the Second World War, the city of Vilnius in Lithuania had over 200,000 residents, about 80,000 of whom were Jewish.
Vilnius was major center for European and Jewish culture and had been described as the ``Jerusalem of the North.'' When war broke out in 1939, Vilnius fell under the control of the
Soviet Union, then briefly Lithuania, and again the Soviets. Nazi Germany captured the city in 1941 and by July of that year began the systematic extermination of the Jewish population, months before they
began their genocide against the Jews elsewhere in Europe. By the end of 1941 it is estimated that about 50,000 Jews had been executed, and
more followed in later years. In total it is estimated that 60,000-70,000 Jews and 30,000 others were executed by the Nazis or Lithuanian collaborators in or near Vilnius during the war,
in what became known as the ``Ponary Massacre'' \cite{sakowicz,shaer}.
%\footnote{''Ponary — the Golgotha of Wilno"='' https://web.archive.org/web/20081224124407/http://www.wsp.krakow.pl/konspekt/konspekt5/ponary.html}

Most of the victims were marched in large groups to the forest outside of Vilnius, stripped, and shot. Their bodies were then dumped into one of six large cylindrical pits,
12 to 30 meters in diameter and 5 to 8 meters deep, that had previously been dug for oil storage tanks by the Soviets, but never used.
By 1943 the Soviet Union was threatening to recapture Vilnius, so the Nazis began an attempt to hide their atrocities. Approximately
80 Jewish prisoners were forced to enter and reside in the pits, exhuming and burning the remains of the bodies, sometimes finding the bodies of their own
wives and children. The ashes were then mixed with sand and scattered in the forest,
eliminating evidence of what the Nazis had done. The prisoners realized that as soon as their work was done, their ashes would also be added to the forest, and so they came up with
an escape plan. Each night under the cover of darkness, over a period of two and a half months, they dug an escape tunnel from one of the pits, using only bare hands
and implements recovered from the bodies. The tunnel
ranged as deep as several meters, just wide enough for one person at a time, and led 30 meters to an opening away from the pit. On a moonless
night in April of 1944 the prisoners quietly exited the pit through this tunnel, but not quietly enough. The guards were alerted and pursued the escapees
with dogs and guns. Only 12 prisoners escaped into the forest. Of these, 11 survived the war to bear witness to what had happened in the forests outside of Vilnius \cite{shaer}.

The stories of these survivors and the escape tunnel they helped dig became a part of the history of the Holocaust. The pits and executions outside of Vilnius have been well-documented, but until recently no evidence or remains of the tunnel had been found. As this is considered hallowed ground, any type of excavation that would disturb victims' remains was out the question.
But in 2016 a team of researchers were able to use new imaging methods to locate and map the remains of the tunnel below ground, without disturbing the area \cite{McClymont},
and so contribute modern forensic evidence to the story of the survivors.
The Public Broadcasting Service NOVA special \cite{nova} has a
fascinating account of the work done to discover this tunnel, or see \cite{bauman} for an interesting day-to-day blog of one of the geophysicists involved in the work.

One of the main methods used to image the remains of the tunnel was \textit{electrical resistivity tomography} (ERT),
a technique used in geophysical imaging. The goal of ERT is to produce an image of the subsurface structure of the earth, a
sort of x-ray or CT scan of what lies below the surface in a localized area. However, ERT doesn't use x-rays, but rather, electrical measurements. It is
closely related to \textit{electrical impedance tomography} (EIT), a method that has been explored
for medical diagnosis and nondestructive testing.
In either technique we want to image the inaccessible interior of an object, for example, subsurface soil or a human body, using electrical current.

Refer to the left panel in Figure \ref{fig0} for an illustration of how ERT works. One injects an electrical
current into an object or region of interest. In Figure \ref{fig0} a current is being injected into the earth at the leftmost electrode,
and withdrawn at the rightmost electrode; the manner in which the current flows through
the earth is affected by the subsurface electrical resistivity, and
this in turn affects the voltage throughout the region, including the surface.  By measuring this
voltage using the other electrodes we can gather information about the subsurface electrical resistivity of the soil. This
might be repeated by using other input/output electrodes, and with this information and appropriate
mathematics we can build a  ``resistivity image'' of the subsurface soil. A tunnel, even collapsed, might manifest its presence in altered resistivity of the disturbed soil.
\begin{figure}[ht]
\centerline{\includegraphics[width=0.5\textwidth]{geoprospect1.png}\ \ \ \includegraphics[width=0.5\textwidth]{eit.png}}
\caption{Left panel: Illustration of ERT for imaging subsurface structure. Right panel: Analogous problem, imaging a hole
of unknown radius inside a circular region.}
\label{fig0}
\end{figure}

In this project we'll look at a special case of ERT in a simplified setting, illustrated in the right panel of Figure \ref{fig0}, in which we seek to determine the radius of
a circular nonconductive hole or void inside a bounded two-dimensional region, instead of an unbounded half-plane. In our setting the region will be the unit disk centered at the
origin (the outer circle in the right panel) and the nonconductive void will be a concentric disk of unknown radius $R<1$ (the black disk in the right panel). We have
access only to the outer surface (the unit circle). We attach electrodes to the outer surface and pass electrical current through the region, illustrated by the red curves,
which are flow lines for the current and must flow around the void. This alters the voltage throughout the conductive region, and in particular on the outer surface. By measuring this voltage on
the outer surface we attempt to determine the radius of the void.

To understand how this type of imaging is done, in this simplified setting or in general, we need to know a bit about how electrical current flows through
a conductive material.

\section*{MODELING ELECTRICAL CONDUCTION}

\subsection*{Current Flux}

First, all of our models will be steady-state, so nothing will depend on time.
There will be no magnetic fields, just electrostatic fields.
We'll also stick to two dimensions; the extension to three dimensions is straightforward. Suppose
that $\Omega$ denotes an electrically conductive object or region in two dimensions; $\Omega$ could be a half-space, if imaging
subsurface features of the earth, or a bounded region, for example, a cross-section of a human body.

To describe the flow
of electrical current through $\Omega$ we will use a vector field $\mathbf{J}(x,y)$, and denote the
rectangular coordinate components of this vector field by $J_1(x,y)$ and $J_2(x,y)$, so $\mathbf{J}(x,y) = \langle J_1(x,y),J_2(x,y)\rangle$.
At each point in $\Omega$ the current flux $\mathbf{J}$ has both a direction and a magnitude.
The direction of $\mathbf{J}$ indicates the direction in which charge is flowing, and we'll take the conventional
approach in which charge is positive, though it doesn't really matter. The magnitude $\|\mathbf{J}\| = \sqrt{J_1^2+J_2^2}$
at a point $\mathbf{p} = (x,y)$ is defined to be the ``intensity'' of the current flowing past a neighborhood of $\mathbf{p}$. This ``intensity'' is quantified
precisely in the following manner: imagine a hypothetical short line segment
``$L$'' centered at $\mathbf{p}$, with length $|L|$ and unit normal vector $\mathbf{n}$, oriented so that $\mathbf{n}$
is parallel to $\mathbf{J}$ at $\mathbf{p}$ (so $L$ is at a right angle to $\mathbf{J}$). Refer to the left panel in Figure \ref{fig1}, with
representative points $\mathbf{p}=\mathbf{p}_1$ or $\mathbf{p}=\mathbf{p}_2$. The red line segment in each case is a short segment $L$, with normal $\mathbf{n}$, oriented
to be parallel to the local flow of $\mathbf{J}$.
The magnitude $\|\mathbf{J}\|$ at $\mathbf{p}$ is defined so that the current
flux over $L$ (charge per unit time) is given approximately as
$$
\textrm{current flux over $L$} \approx |L| \|\mathbf{J}\|,
$$
at least if $L$ is very short. More precisely, the magnitude $\|\mathbf{J}\|$ is defined as
\begin{equation}
\label{flux0}
\|\mathbf{J}\| = \lim_{|L|\rightarrow 0} \left (\frac{\textrm{current flux over $L$}}{|L|}\right ).
\end{equation}
This means $\mathbf{J}$ has units of charge per time per length.
\begin{figure}[ht]
\centerline{\includegraphics[width=0.4\textwidth]{conduct1.png}\ \ \ \includegraphics[width=0.4\textwidth]{conduct2.png}}
\caption{Left panel: Current flux field with points $\mathbf{p}_1,\mathbf{p}_2$ and short reference segments (red) with unit normal vectors in the same
direction as the current flux.
Right panel: point $\mathbf{p}$ with reference segment (red) and unit normal vector at arbitrary orientation.}
\label{fig1}
\end{figure}

On the right in Figure \ref{fig1} we show the more general case with a point $\mathbf{p}$ and short line segment $L$ that is not oriented at a right angle
to the local current field $\mathbf{J}$. In this case the flux of $\mathbf{J}$ over $L$ is equal to $\|\mathbf{J}\| |L| \cos(\theta)$, where
$\theta$ is the angle between the unit normal $\mathbf{n}$ and the current flux $\mathbf{J}$ at $\mathbf{p}$. If $\theta=0$ this is just (\ref{flux0}).
See Problem \ref{prob1} for some intuition about cases where $\theta>0$.
Since $\|\mathbf{n}\|=1$ we can write the flux in the general case as $|L|\|\mathbf{J}\|\|\mathbf{n}\|\cos(\theta)$ or, using
the basic properties of the dot product, as
\begin{equation}
\label{flux1}
\textrm{current flux over $L$} \approx |L|\|\mathbf{J}\|\|\mathbf{n}\|\cos(\theta) = |L|(\mathbf{J}\cdot \mathbf{n}).
\end{equation}
\begin{problem}
\label{prob1}
Sketch a situation in which the current flux $\mathbf{J}$ at a point $\mathbf{p}$ cuts orthogonally across a short line segment $L$
centered on $\mathbf{p}$ and $\mathbf{n}$ points in the same direction at $\mathbf{J}$. In this case $\theta=0$. Repeat for the case in which
$\mathbf{n}$ is orthogonal to $\mathbf{J}$ at $\mathbf{p}$. What is $\theta$ here, and why does (\ref{flux1}) make physical sense in this case? Suppose that $\mathbf{n}$
is exactly opposed to $\mathbf{J}$---what is $\theta$ now? What does (\ref{flux1}) give for the current flux over $L$, and how does it compare to the case $\theta=0$?
\end{problem}

\subsection*{Conservation of Charge}

An important fact is that if electric charge is conserved then
\begin{equation}
\label{chargecon}
\textrm{div}(\mathbf{J}) = 0
\end{equation}
inside $\Omega$, where $\textrm{div}(\mathbf{J})$ is the \textit{divergence of $\mathbf{J}$}. In rectangular coordinates
$\textrm{div}(\mathbf{J}) = \frac{\partial J_1}{\partial x} + \frac{\partial J_2}{\partial y}$.

To see why (\ref{chargecon}) is true, suppose that $C$ is a simple smooth curve inside $\Omega$, and let $\mathbf{n}$
be a unit normal vector on $C$, consistently oriented (so as we move along $C$ the vector
$\mathbf{n}$ changes continuously, and always points to the same ``side'' of $C$.) The total current flux of $\mathbf{J}$
over $C$ in the direction of $\mathbf{n}$ is given by the integral
\begin{equation}
\label{fluxgen}
\textrm{flux of $\mathbf{J}$ over $C$} = \int_C \mathbf{J}\cdot\mathbf{n}\,ds
\end{equation}
where $ds$ is arc length on $C$. Equation (\ref{fluxgen}) follows by chopping the curve $C$ up into many short pieces of length ``$ds$'', applying (\ref{flux1}) to each piece
(yielding net flux $\mathbf{J}\cdot\mathbf{n}\,ds$ over each piece) and summing (integrating) the results.

If $C$ is also a closed curve that bounds a region $D$ and charge is conserved then, since
we're in the steady-state case, it must be that the current flux over $C$ is zero, for otherwise charge
would be building up or being depleted in $D$, which is not steady-state. Thus
\begin{equation}
\label{conserve1}
\int_C \mathbf{J}\cdot\mathbf{n}\,ds = 0.
\end{equation}

If we apply the Divergence Theorem and use (\ref{conserve1}) we find
\begin{equation}
\label{conserve2}
\int_D \textrm{div}(\mathbf{J})\,dA = \int_C \mathbf{J}\cdot\mathbf{n}\,ds = 0
\end{equation}
where $dA$ is the area element.
But the curve $C$ is an arbitrary simple closed curve, so $D$ on the left in (\ref{conserve2}) could be any region with such a boundary,
for example, any disk. It is a fact that if $f(x,y)$ is a continuous function defined in some two-dimensional
region and
\begin{equation}
\label{con3}
\int_D f(x,y)\,dA = 0
\end{equation}
for \textit{every} disk $D$ then $f$ must be identically zero. (To see this, suppose $f(x_0,y_0)>0$ at some point
$(x_0,y_0)$. Then $f(x,y)>0$ for all $(x,y)$ close enough to $(x_0,y_0)$, say within a distance $\epsilon$, since $f$ is continuous. Take $D$
to be a disk of radius $\epsilon$ centered at $(x_0,y_0)$ and (\ref{con3}) would be violated since the integral of a positive function must be positive.)
Using this reasoning we see that the integrand in the left integral of (\ref{conserve2}) is forced to be identically zero in $\Omega$, and so (\ref{chargecon}) must hold
inside the conductive region.\\\ \\
\begin{problem}
\label{prob2}
Use (\ref{chargecon}) to decide which of the following vector fields $\mathbf{J}$ could represent the flux for a steady-state current flow
in some region $\Omega$.
\begin{enumerate}
\item[a.] $\mathbf{J}(x,y) = <1,3>$.
\item[b.] $\mathbf{J}(x,y) = <x,y>$.
\item[c.] $\mathbf{J}(x,y) = <x,-y>$.
\item[d.] $\mathbf{J}(x,y) = <x^2y,-y>$.
\item[e.] $\mathbf{J}(x,y) = <xy,-y^2/2>$.
\end{enumerate}
\end{problem}
\begin{problem}
\label{prob3}
Suppose that $\mathbf{J}(x,y)=<x+x^2y, f(x,y)>$ represents the current flux of a steady-state flow. What can you deduce about $f(x,y)$?
\end{problem}

\subsection*{Ohm's Law in Two (or Three) Dimensions}

Recall Ohm's Law for a resistor: $V=IR$, where $V$ is the voltage (potential)
drop across the resistor, $I$ is the current through the resistor in the direction of the voltage drop, and $R$ is the \textit{resistance} of the resistor. It is the drop from high potential at one end of the resistor to a lower potential at the other end that indicates the presence of an
electric field in the resistor, which impels charge to move and current to flow.
The rate at which current flows is (often) proportional to the voltage drop, that is, $I = V/R$, hence Ohm's Law. In this context $1/R$ is called the \textit{conductivity}
of the resistor.

We will generalize Ohm's Law to two (or three) dimensions.
We now make a modeling assumption.  Electrical current is just moving charge.  Why
would charge move?  Charges move because they are being pushed by an electrical field.  In the simplest case a
positively charged particle will move in the direction of the electric field.  Based on
this reasoning, we are going to assume that inside whatever electrical conductor
we're dealing with we have
\begin{equation}
\label{ohm2}
\mathbf{J} = \gamma \mathbf{E}
\end{equation}
where $\mathbf{E}$ is the electric field inside $\Omega$ and $\gamma \geq 0$ is called the
\textit{electrical conductivity} of $\Omega$.  The scalar quantity $\gamma$ could be constant
throughout the object, or it could be a function of position.  It doesn't matter for now. Equation (\ref{ohm2}) will
be the two-dimensional analogue of Ohm's Law. It says that charge flows in the direction it's being pushed by $\mathbf{E}$, in proportion
to the magnitude of $\mathbf{E}$.

You should think about why $\gamma$ is called the ``conductivity.''  For a fixed electrical
field $\mathbf{E}$, large values of $\gamma$ (high conductivity) result in larger $\mathbf{J}$,
while small $\gamma$ results in little current flow.  The case $\gamma\equiv 0$
corresponds to an insulator, in which no current flows at all, regardless of the electric
field.  On the other hand if ``$\gamma=\infty$'' (a perfect
conductor) then we must have $\mathbf{E}=0$, or else unlimited current would flow.  This is the
(possibly familiar) fact that the electric field inside a perfect conductor is zero.
The quantity $1/\gamma$ is a called the \textit{resistivity}, and of course may
be a function of position.

%The model in equation (\ref{ohm2}) is just that---a model. It's really a higher dimensional
%version of Ohm's Law, that current is proportional to the electric field strength.  There
%are materials that don't obey equation (\ref{ohm2}), so our model wouldn't work in
%those cases.

%So based on equations (\ref{chargecon}) and (\ref{ohm2}) we have $\nabla \cdot (\gamma \mathbf{E}) = 0$.

\subsection*{The Equation for the Electric Potential}

We will use $\phi(x,y)$ for the electric potential (voltage) inside $\Omega$. Recall that
potential is always ``relative'' to some fixed zero point, which we'll worry about later.
Our ultimate goal is to derive a differential equation that $\phi$ must satisfy.

It is a fact that for an electrostatic field $\mathbf{E}$ we have $\mathbf{E} = -\nabla \phi$
where $\phi$ is the potential or voltage function; that is, $\mathbf{E}$ is a \textit{conservative} vector field. From (\ref{ohm2}) this means that
\begin{equation}
\label{currvolt}
\mathbf{J} = -\gamma \nabla \phi.
\end{equation}
Combining this with (\ref{chargecon}) shows that $\phi$ must satisfy
\begin{equation}
\label{condeqn}
\textrm{div}(\gamma \nabla \phi) = 0
\end{equation}
inside $\Omega$. This is a partial differential equation (PDE) for $\phi$.
Keep in mind that in two dimensions $\phi$ and $\gamma$ are functions of $x$ and $y$, while in
three dimensions they're functions of $x,y$, and $z$.  You can write
out (\ref{condeqn}) explicitly in two dimensional rectangular coordinates, as
\begin{equation}
\label{condeqn2}
\gamma \left ( \frac{\partial^2 \phi}{\partial x^2} + \frac{\partial^2 \phi}{\partial y^2}\right ) + \frac{\partial\phi}{\partial x}
\frac{\partial\gamma}{\partial x} + \frac{\partial\phi}{\partial y}
\frac{\partial\gamma}{\partial y} = 0.
\end{equation}

In the special case in which $\gamma>0$ is constant, (\ref{condeqn}) becomes $\textrm{div} (\nabla \phi) = 0$, or
in two dimensions with rectangular coordinates,
\begin{equation}
\label{condeqn2b}
\frac{\partial^2 \phi}{\partial x^2} + \frac{\partial^2 \phi}{\partial y^2} = 0.
\end{equation}
Equation (\ref{condeqn2b}) is \textit{Laplace's equation}, one of the most common and important equations of mathematical physics.
A function $\phi(x,y)$ that satisfies (\ref{condeqn2b}) is said to be \textit{harmonic}.

Given the modeling assumptions we have made, (\ref{condeqn}) must be satisfied
at all points $(x,y)\in \Omega$.  Thus if we knew the conductivity $\gamma$
inside the object, we could try to solve the PDE (\ref{condeqn}) to find $\phi$.  But there still
isn't enough information to nail down a unique solution for the potential $\phi(x,y)$.  We must also specify boundary
conditions, either the voltage $\phi$ or the input current, on
the boundary of $\Omega$.

\subsection*{Boundary Conditions---Injecting a Current Flux}

We will use $\partial\Omega$ to denote the boundary of $\Omega$. For example, if $\Omega$ is a two-dimensional half-space $\{(x,y); y<0\}$ then $\partial \Omega$
consists of the $x$ axis; if $\Omega$ is a disk, $\partial\Omega$ is the outer boundary, a circle.
In our application we inject a specified current into $\partial\Omega$ and then measure the resulting potential (voltage) on $\partial\Omega$.
Refer to Figure \ref{fig2}, in which a hypothetical current field is shown in a portion of $\Omega$ and a corresponding portion of
$\partial\Omega$ is shown as a quarter circle. At the point $\mathbf{p}$ on $\partial\Omega$ the current flux is generally outward from $\Omega$.
A short portion of $\partial\Omega$ near $\mathbf{p}$ is shown in red; this short portion is, to good approximation, a line segment $L$, with indicated outward unit normal vector $\mathbf{n}$. As per equation (\ref{flux1}) above, the rate at which current is leaving $\Omega$ over $L$ is $|L|(\mathbf{J}\cdot\mathbf{n})$, or equivalently, the rate at which current is entering $\Omega$ over $L$ is minus this, so
\begin{equation}
\label{rateentry}
\textrm{rate current enters $\Omega$ across $L$} = -|L|(\mathbf{J}\cdot\mathbf{n}).
\end{equation}

When we inject current into a region $\Omega$ at the boundary it is conventional to specify the amount of charge injected per unit time (current)
\textit{per length of the boundary} at any point on $\partial\Omega$. That is to say, we specify $-\mathbf{J}\cdot\mathbf{n}$
on $\partial\Omega$ (so $|L|$ is divided out in (\ref{rateentry})). Let $g$ be such a specified input current density, so that our boundary condition becomes
$-\mathbf{J}\cdot\mathbf{n} = g$ on $\partial\Omega$. If we make use of (\ref{currvolt}) we obtain a boundary condition
\begin{equation}
\label{bcinput}
\gamma \nabla \phi \cdot \mathbf{n} = g
\end{equation}
for $\phi$ on $\partial\Omega$. A boundary condition of the form (\ref{bcinput}), in which we specify information related to $\nabla\phi\cdot\mathbf{n}$, the normal derivative of $\phi$, is called a \textit{Neumann boundary condition}.
\begin{figure}[ht]
\centerline{\includegraphics[width=0.5\textwidth]{conduct3.png}}
\caption{Current flux field on a portion of $\Omega$, with a short ``linear'' piece of $\partial\Omega$ near $\mathbf{p}$, unit outward normal vector $\mathbf{n}$
Charge exits $\Omega$ at a rate of $|L|(\mathbf{J}(\mathbf{p})\cdot\mathbf{n})$ Coulombs per second.}
\label{fig2}
\end{figure}

\subsection*{The Forward Problem}

For a given region $\Omega$ with boundary $\partial\Omega$, a specified conductivity $\gamma(x,y) > 0$ defined on $\Omega$, and a specified input current flux function
$g$ defined on $\partial\Omega$, equation (\ref{condeqn}) and the boundary condition (\ref{bcinput}) constitute a \textit{boundary value problem}, to be solved for the
function $\phi(x,y)$. Physically, we inject a current flux $g$ into a region $\Omega$ with internal conductivity $\gamma$ and seek to find $\phi(x,y)$, the electrical
potential or voltage at all points in $\Omega$.
This is the traditional thing to do with a differential equation and boundary data---find the solution. We refer to this as the \textit{forward} or \textit{direct} problem.

In the case that $\Omega$ is a bounded region it turns out that equation (\ref{condeqn}) with the boundary condition (\ref{bcinput}) always has a solution if $g$
satisfies a certain condition (see Problem \ref{prob7} below), but the solution is not quite unique; you can easily check that if $\phi(x,y)$ is any solution then so is $\phi(x,y)+c$ for any constant $c$. However, this is the only way solutions can fail to be unique: if $\phi_1$ and $\phi_2$ both
satisfy (\ref{condeqn}) with the boundary condition (\ref{bcinput}) then $\phi_2=\phi_1+c$ for some constant $c$. Physically, this means that potentials are only
unique up to an arbitrary additive constant (determined by the ``ground'' reference point.) If $\Omega$ is unbounded (e.g., a half plane) then we also need to specify that the solution
$\phi$ is bounded as $(x,y)\rightarrow\infty$ to obtain a unique solution.
\begin{problem}
\label{prob4}
Suppose that $\Omega$ is the square $0<x,y<1$ in the plane. We input a current flux $g$ equal to $-1$ on the left side of $\Omega$ (where $x=0$), $0$
on the top and bottom of $\Omega$, and $g=1$ on the right ($x=1$) side of $\Omega$. Physically, we are injecting a current of $1$ ampere per meter of boundary on the right edge
of $\Omega$ and withdrawing $1$ ampere per meter of boundary on the left edge; the top and bottom have no current input or output, i.e., they are insulated.  Verify that the function $\phi(x,y) = x+c$ satisfies equation (\ref{condeqn}) with $\gamma\equiv 1$, and that the boundary condition (\ref{bcinput}) holds, for any constant $c$. Hint: the unit outward normal vector $\mathbf{n}$ equals $<-1,0>$ on the left side of $\Omega$ and is equally easy to find on the other sides.
\end{problem}
\begin{problem}
\label{prob5}
Suppose that $\Omega$ is the square $0<x,y<1$ in the plane, with conductivity $\gamma(x,y) = 1+xy$. We input a current flux $g$ given by
$$
g = \left \{\begin{array}{ll}
0, & \textrm{  on  }x=0\textrm{  (the left side of $\Omega$)}\\
  -2x-2, & \textrm{  on  }y=1\textrm{  (the top of $\Omega$)}\\
    2y+2, & \textrm{  on  }x=1\textrm{  (the right side of $\Omega$)}\\
    0, & \textrm{  on  }y=0\textrm{  (the bottom of $\Omega$)}\\
\end{array}
\right .
$$
Verify that the function $\phi(x,y) = x^2-y^2+c$ satisfies equation (\ref{condeqn}) with $\gamma(x,y) = 1+xy$, and that the boundary condition (\ref{bcinput}) holds,
for any constant $c$.
\end{problem}
\begin{problem}
\label{prob6}
Show that if $\phi(x,y)$ is a solution to (\ref{condeqn}) with boundary condition (\ref{bcinput}) then so is $\phi(x,y)+c$ for any constant $c$.
\end{problem}
\begin{problem}
\label{prob7}
Show that if $\textrm{div}(\mathbf{J})=0$ holds in $\Omega$ then it must be the case that
$$
\int_{\partial\Omega} \mathbf{J}\cdot\mathbf{n}\,ds = 0
$$
where $\mathbf{n}$ is an outward unit normal vector on $\partial\Omega$. Hint: Look at (\ref{conserve1}). Use this to deduce that if (\ref{bcinput}) holds then it must be the
case that
$$
\int_{\partial\Omega} g\,ds = 0.
$$
Thus equation (\ref{condeqn}) with the boundary condition (\ref{bcinput}) is solvable only if $g$ integrates to zero around $\partial\Omega$. What is the
physical interpretation of this condition? Hint: in a steady-state problem, what goes in must come out.
\end{problem}

\subsection*{Solution to the Forward Problem on an Annulus}

In general the solution to (\ref{condeqn}) can't be written out explicitly, except in special cases. One such important special case is that of constant conductivity and
a simple geometry for $\Omega$, for example, a rectangle, disk, or annulus. It is the last geometry that we'll consider, with conductivity $\gamma \equiv 1$. This
simplified setting will serve to illustrate the essential idea of ERT: it is possible to determine something about the interior electrical properties or geometry of a region from
input current/output voltage data.

So to narrow our focus, we suppose that $\Omega$ is an annulus.
Given the circular symmetry of the annulus, it's most convenient to work in polar coordinates $(r,\theta)$. Let $\Omega$ have inner radius $R$
and outer radius $1$, so in polar coordinates with $r=0$ as the center of the annulus we have
\begin{equation}
\label{Omegadef}
\Omega = \{(r,\theta); R<r<1,-\pi<\theta\leq \pi\}.
\end{equation}
Think of $\Omega$ as an object to which we have access only at $r=1$, the outer surface. Our goal is to determine the radius $R$ of the inside ``hole'' from
data taken on the outside at $r=1$. We will assume that $\Omega$ has constant conductivity $1$. The hole $r<R$ is to be interpreted as a nonconductive void, through which
electrical current cannot pass. One might think of it as a cross-section of a cylindrical tunnel.

To image the inside of $\Omega$ we apply an input current flux $g$ at $r=1$; the function $g$ thus varies only with the angular variable $\theta$, that is, $g=g(\theta)$.
This causes current to flow through $\Omega$, and an electric potential is induced throughout $\Omega$.
Let's write this potential as $\phi(r,\theta)$, a function of polar radius $r$ and angle $\theta$. In polar equations Laplace's equation,
which $\phi$ must satisfy, becomes
\begin{equation}
\label{lappolar}
\frac{\partial^2\phi}{\partial r^2} + \frac{1}{r}\frac{\partial \phi}{\partial r} + \frac{1}{r^2}\frac{\partial^2 \phi}{\partial\theta^2} = 0.
\end{equation}
The usual process of separating variables, as $\phi(r,\theta) = p(r)q(\theta)$ leads us to separable solutions to Laplace's equation of the form
$$
1, \ln(r), r^k\sin(k\theta), r^k\cos(k\theta), r^{-k}\sin(k\theta), r^{-k}\cos(k\theta)
$$
where $k\geq 1$ is an integer. Moreover, any linear combination of these basic functions is also a solution to Laplace's equation. Note that since $r>R$ in $\Omega$,
the function $\ln(r)$ and $r^{-k}$ pieces are valid possibilities, since we cannot approach $r=0$.  We thus consider solutions that consist of a linear combination
\begin{equation}
\label{sol1}
\phi(r,\theta) = c_0 + d_0\ln(r) + \sum_{k=1}^{\infty} \left ((c_k r^k + c_{-k}r^{-k})\cos(k\theta) + (d_k r^k + d_{-k}r^{-k})\sin(k\theta)\right )
\end{equation}
of the above functions, for appropriate constants $c_k,d_k$, $-\infty<k<\infty$.
By adjusting the $c_k$ and $d_k$ for $-\infty<k<\infty$ we will be able to obtain the required boundary conditions.

On the outer boundary $r=1$ it's easy to see that the unit normal vector $\mathbf{n}$ points radially outward; the quantity $\nabla \phi\cdot \mathbf{n}$ is thus
just the directional derivative of $\phi$ in the radial direction, that is, $\nabla \phi\cdot \mathbf{n}=\frac{\partial \phi}{\partial r}$; we want this quantity
to equal $g(\theta)$ on the outer boundary $r=1$. By differentiating the right side of (\ref{sol1})
term-by-term with respect to $r$ we find that, from the power rule,
\begin{equation}
\label{sol1dr}
\frac{\partial\phi}{\partial r}(r,\theta) = d_0/r + \sum_{k=1}^{\infty} \left ((kc_k r^{k-1} - kc_{-k}r^{-k-1})\cos(k\theta) + (kd_k r^{k-1} - kd_{-k}r^{-k-1})\sin(k\theta)\right ).
\end{equation}
(The constant $c_0$ drops out.)  At $r=1$ we find
\begin{equation}
\label{sol1dra}
\frac{\partial\phi}{\partial r}(1,\theta) = d_0 + \sum_{k=1}^{\infty} \left ((kc_k - kc_{-k})\cos(k\theta) + (kd_k - kd_{-k})\sin(k\theta)\right ).
\end{equation}
We must choose the $c_k$ and $d_k$ so that the expression on the right equals $g(\theta)$. Now $g(\theta)$ itself has a Fourier sine/cosine expansion of the form
$$
g(\theta) = a_0/2 + \sum_{k=1}^{\infty} (a_k\cos(k\theta) + b_k\sin(k\theta))
$$
where
\begin{equation}
\label{gcoefs}
a_k = \frac{1}{\pi}\int_{-\pi}^{\pi} g(\theta)\cos(k\theta)\,d\theta,\;\;\;b_k = \frac{1}{\pi}\int_{-\pi}^{\pi} g(\theta)\sin(k\theta)\,d\theta.
\end{equation}
In the present case $a_0=0$ (recall Problem \ref{prob7}).
In order for the right side of (\ref{sol1dra}) to match $g(\theta)$ we must match the remaining Fourier sine/cosine coefficients, so that we need $d_0 = a_0/2=0$ and
\begin{align}
k(c_k-c_{-k}) & = a_k\label{bceqn1}\\
k(d_k-d_{-k}) & = b_k\label{bceqn2}
\end{align}
for $k\geq 1$. For each integer $k\geq 1$ equation (\ref{bceqn1}) is a single equation for two unknowns $c_k$ and $c_{-k}$. The same holds true for
equation (\ref{bceqn2}) and unknowns $d_k,d_{-k}$. We need another equation for each pair of unknowns.

The additional equations come from a boundary condition on $r=R$.
On this inside radius $r=R$ no electrical current can enter or exit $\Omega$, since the region $r<R$ is nonconductive. On $r=R$ the outward unit normal
vector $\mathbf{n}$ to $\partial \Omega$ points directly toward the origin, and in fact $\nabla \phi\cdot\mathbf{n} = -\frac{\partial\phi}{\partial r}(R,\theta)$.
Since no current can cross the boundary $r=R$ we must have $-\frac{\partial\phi}{\partial r}(R,\theta)=0$. Making use of (\ref{sol1dr}) with $r=R$ leads us to
\begin{equation}
\label{sol1drb}
\frac{\partial\phi}{\partial r}(r,\theta) = \sum_{k=1}^{\infty} \left (k(c_k R^{k-1} - c_{-k}R^{-k-1})\cos(k\theta) + k(d_k R^{k-1} - d_{-k}R^{-k-1})\sin(k\theta)\right ) = 0
\end{equation}
since we deduced that $d_0=0$ above. In order for the sum in (\ref{sol1drb}) to be identically zero the Fourier sine/cosine coefficients must all be zero, so that
\begin{align}
k(c_k R^{k-1} - c_{-k}R^{-k-1}) & = 0\label{bceqn1b}\\
k(d_k R^{k-1} - d_{-k}R^{-k-1}) & = 0\label{bceqn2b}
\end{align}
for each integer $k\geq 1$. Note that $k$ can be divided out of each equation above.

Equations (\ref{bceqn1}) and (\ref{bceqn1b}) are a pair of equations for $c_k$ and $c_{-k}$, for each $k\geq 1$, with solution
\begin{align}
c_k & = \frac{R^{-k}}{k(R^{-k}-R^k)}a_k\label{cksol}\\
c_{-k} & = \frac{R^{k}}{k(R^{-k}-R^k)}a_k.\label{ck1sol}
\end{align}
Similar reasoning shows that we need
\begin{align}
d_k & = \frac{R^{-k}}{k(R^{-k}-R^k)}b_k\label{dksol}\\
d_{-k} & = \frac{R^{k}}{k(R^{-k}-R^k)}b_k.\label{dk1sol}
\end{align}

In summary, we can solve Laplace's equation (\ref{lappolar}) on the annulus $\Omega$ defined by (\ref{Omegadef}) with boundary data $\frac{\partial u}{\partial r}(1,\theta) = g(\theta)$ (physically, a current flux injected into the outer boundary) and $\frac{\partial u}{\partial r}(R,\theta) = 0$ (a nonconductive void of radius $R$ in the middle)
with the following procedure:
\begin{enumerate}
\item Compute all Fourier sine/cosine coefficients $a_k,b_k$ for $g(\theta)$ via (\ref{gcoefs}) for $k\geq 1$ (note $a_0=0$).
\item Set $d_0=0$ and compute $c_k,d_k$ for all integers $k\neq 0$ according to (\ref{cksol})-(\ref{dk1sol}).
\item Form solution $\phi(r,\theta)$ via (\ref{sol1}). Here $c_0$ is an arbitrary constant that is dictated by our reference or ``ground'' value for $\phi$.
It is common to chose $c_0$ so that ${\displaystyle \int_{-\pi}^{\pi} \phi(1,\theta)\,d\theta = 0}$ (the solution is normalized to have zero integral over the
outer boundary $r=1$). This dictates the value $c_0=0$.
\end{enumerate}

Let's use the procedure above to work an example.
\begin{example}
\label{solexamp1}
Suppose that the input current flux is given by $g(\theta) = 2\cos(\theta)$ on $r=1$ and that the inside hole has radius $R=1/2$. Note that $g$ integrates to zero over this curve, that is, $\int_{-\pi}^{\pi} g(\theta)\,d\theta = 0$, so the boundary value problem is solvable.
In Step 1 above we find Fourier sine/cosine coefficients $a_1=2$ and all other $a_k=0$, as well as all $b_k=0$ (that is, happily, $g$ has a finite Fourier sine/cosine
expansion). In Step 2 using (\ref{cksol}) shows that $c_1=8/3,c_{-1}=2/3$, and
all other $c_k=0$, as well as $d_k=0$ for all $k$. In Step 3 we make use of(\ref{sol1}) to see that the solution to Laplace's equation with this outer boundary data is
$$
\phi(r,\theta) = \frac{8}{3} r\cos(\theta) + \frac{2}{3r}\cos(\theta).
$$
It's straightforward to check that $\frac{\partial\phi}{\partial r}(1,\theta)=2\cos(\theta)$ and $\frac{\partial\phi}{\partial r}(1/2,\theta)=0$.
\end{example}
\begin{problem}
\label{prob8}
Suppose $R=1/3$, $\gamma\equiv 1$, and the input flux is $g(\theta) = 2\sin(4\theta)$. Use the procedure illustrated in Example \ref{solexamp1}) to compute $\phi(r,\theta)$.
\end{problem}
\section*{THE INVERSE PROBLEM AND ERT}

The recipe above and its application in Example \ref{solexamp1} illustrates the classic task in differential equations: Given a domain, a differential equation, and
boundary (and/or initial) data, find the solution. In ERT we're interested in a slight twist on this, the \textit{inverse problem}. We are actually given the solution,
or at least some information about the solution, and asked to find coefficients in the differential equation, or deduce something about the geometry of the domain
that we can't see from the outside. In the present case we will be given the input
current flux $g(\theta)$, the ``measured'' potential $\phi(1,\theta)$ on the outer boundary of $\Omega$, and then asked to deduce the inner radius $R$, under the
assumption that $\Omega$ has a constant conductivity $\gamma=1$.
\begin{example}
\label{invexamp1}
Suppose that, as in Example \ref{solexamp1}, we are told that a current flux $g(\theta)=2\cos(\theta)$ was injected into $\Omega$ and the resulting
boundary potential was measured to be $\phi(1,\theta)=\frac{10}{3} \cos(\theta)$ (the solution from Example \ref{solexamp1} at $r=1$). The radius $R$ of the void (tunnel?) is unknown; possibly even
$R=0$, so the void isn't present. Can we deduce the value of $R$?

The solution (\ref{sol1}) with the associated recipe for the $c_k$ and $d_k$ is still perfectly valid. We simply don't know $R$. But we do know from
$g(\theta)=2\cos(\theta)$ that $a_1=2$ and all other $a_k=0$, as well as all $b_k=0$, just as in Example \ref{solexamp1}. The formulae (\ref{cksol})-(\ref{dk1sol}) for the $c_k$ and $d_k$ then dictate that
\begin{align}
c_1 & = \frac{2R^{-1}}{(R^{-1}-R)},\label{cksolaa}\\
c_{-1} & = \frac{2R}{(R^{-1}-R)}\label{ck1solaa}
\end{align}
with all other $c_k=0$ and all $d_k=0$. We might be able to use (\ref{cksolaa}) or (\ref{ck1solaa}) and determine $R$ if we know something about $c_1$ and/or $c_{-1}$,
which we do! Equation (\ref{sol1}) indicates that
$$
\phi(r,\theta) = (c_1 r + c_{-1}/r)\cos(\theta).
$$
When $r=1$ as on the outer boundary of $\Omega$ this becomes
$$
\phi(1,\theta) = (c_1+c_{-1}) \cos(\theta).
$$
Since we ``measured'' $\phi(1,\theta)=\frac{10}{3} \cos(\theta)$ we conclude that
$$
c_1+c_{-1} = \frac{10}{3}
$$
after cancelling the cosine terms. With (\ref{cksolaa}) and (\ref{ck1solaa}) this becomes
\begin{equation}
\label{Reqn}
\frac{2(R^{-1}+R)}{(R^{-1}-R)} = \frac{10}{3},
\end{equation}
an equation for $R$ with unique nonnegative solution $R=1/2$ ($R=-1/2$ also works, but our context demands $R\geq 0$).
If it's there, we can find the radius of the void!
\end{example}

The natural questions to ask about this inverse problem are those one asks for any inverse problem. Among these are\\\ \\
{\bf Uniqueness:} Does the given data uniquely determine the unknown quantity? If not, what else is needed?\\
{\bf Reconstruction:} How can we efficiently find the unknown from the given data?\\
{\bf Stability:} If the data are noisy or erroneous, how will this affect our estimate of the unknown?\\\ \\
In Example \ref{invexamp1} above we saw that the value of $R$ was uniquely determined by equation (\ref{Reqn}), and easy to reconstruct (it results in a quadratic equation for $R$).
In Problems \ref{prob10}-\ref{prob11} below you can examine a more general case.\\\ \\
\begin{problem}
\label{prob9}
Suppose that the input current flux $g(\theta) = \cos(3\theta)$; the resulting boundary voltage is $\phi(1,\theta) = \frac{4097}{12285}\cos(3\theta)$ (no noise!)
Determine $R$ by emulating Example \ref{invexamp1}. A computer algebra system might help. Is $R\geq 0$ uniquely determined?\\\\ \\
\end{problem}
\begin{problem}
\label{prob10}
Suppose that the input current flux $g(\theta) = \cos(m\theta)$ for some integer $m\geq 1$; we measure the boundary voltage $\phi(1,\theta)$ for $-\pi<\theta\leq \pi$.
\begin{enumerate}
\item[a.] Suppose the true inner radius is $R_0$ (unknown, but $0<R_0<1$). Solve the forward problem to show that the resulting boundary voltage is $\phi(1,\theta) = q_m \cos(m\theta)$ where
$$
q_m = \frac{1+R_0^{2m}}{m(1-R_0^{2m})}.
$$
Show that $q_m>1/m$.

\item[b.] Show that to solve the inverse problem we must solve the equation
$$
\frac{1+R^{2m}}{m(1-R^{2m})} = q_m
$$
for unknown $R$, and that there is only one nonnegative solution for $R$ if $q_m>1/m$. Hint: Try substituting $T=R^{2m}$ so that the equation above becomes
\begin{equation}
\label{Teqn}
\frac{1+T}{m(1-T)} = q_m.
\end{equation}
When $T=0$ the left side of (\ref{Teqn}) equals $1/m$; as $T\rightarrow 1^-$ the left side approaches $\infty$. Differentiate the left side of (\ref{Teqn}) with
respect to $T$ and show that $(1+T)/(1-T)$ is strictly increasing for $0<T<1$, and conclude that (\ref{Teqn}) has a unique solution for $T$. Since $R^{2m}=T$, what does this say about $R$?
\end{enumerate}
\end{problem}
\begin{problem}
\label{prob11}
Suppose the true inner radius is $R_0$.
Show that if $g$ is any input flux that is not identically zero (but that integrates to zero around $\partial\Omega$, of course) then the data $\phi(1,\theta)$ allows us to uniquely determine $R_0$. Hint: If $g$ is not identically zero then at least one Fourier coefficient $a_m$ or $b_m$ is not zero.
\end{problem}

\section*{EXTENSIONS AND CONCLUSION}

We've considered a rather restricted case of ERT, but it does illustrate the essential idea---one can use input current/output voltage measurements on the boundary of a region
to deduce information about the interior. In the present case we considered the unknown to be a nonconductive circular void centered at $r=0$ in the unit disk, but this can be relaxed to allow the
void to be centered elsewhere, or not a disk, or with some nonzero electrical conductivity. Of course $\Omega$ does not have to be the unit disk. Depending on how many
quantities are considered unknown, it may take multiple input current/output boundary voltage measurements to find the unknowns, if it is even possible.
Also, the solutions to the boundary value problem will not have a simple closed form, so more sophisticated or approximate methods may be needed.
The most general extension allows the unknown interior conductivity to be a function of position, so $\gamma = \gamma(x,y)$ (or $\gamma(x,y,z)$ in three dimensions).
In this case multiple input/output data sets will definitely be required to find any significant information about $\gamma$ and form a resistivity image. See \cite{loke} for a survey of methods
used for ERT, or \cite{cheney} for an overview of the closely related EIT problem.

\teacher{
\section*{NOTES FOR TEACHERS}

This is definitely more sophisticated material than is seen in an introductory course on differential equations.
I've used this material in both my Boundary Value Problems course and a Vector Calculus course (mostly third year students), after we've discussed Laplace's equation, the types of boundary conditions, and
how to solve Laplace's equation on a disk or annulus. Most of our students have been exposed to ``Freshman Electromagnetics'' (e.g., electric fields, potentials) so that helps. I usually take them through the basic model leading up to equation (\ref{condeqn}) in class, then assign the problems as homework.

Generalizing this material to recovering a circular region of constant (nonzero) conductivity $\gamma\neq 1$, still centered at the origin, isn't too hard. Generalizing it
to a circular region centered at an arbitrary point $(a,b)$ seems more challenging, but would be a nice research project for a dedicated undergraduate. In two-dimensions
it would probably help to use a conformal map of the unit disk to itself to reduce the problem to that of a region centered at the origin (but the author hasn't tried it.)



}



\teacher{

\section*{SOME ANSWERS}
\begin{enumerate}
\item The first case is illustrated in Figure \ref{fig1}, left panel, and we obtain flux $|L|\|\mathbf{J}\|$ as in the displayed equation before (\ref{flux0}). When $\mathbf{n}$ is
    orthogonal to $\mathbf{J}$ at $\mathbf{p}$ we have $\theta=\pi/2$, and (\ref{flux1}) yields a flux of zero over $L$. This makes sense since here $L$ is parallel to $\mathbf{J}$, so the current just skims over $L$, but doesn't cross. In the last case in which $\mathbf{n}$ is opposed to $\mathbf{J}$ we have $\theta=\pi$. Here the flux over $L$ is
    $|L|\mathbf{J}\cdot\mathbf{n} = -|L|\|\mathbf{J}\|$. This indicates that the flux is exactly the same as when $\theta=0$, but negated. The charge is crossing $L$ in the direction opposite $\mathbf{n}$.

\item The vectors fields that might represents a steady-state current flux are those of (a), (c), and (e), for these satisfy $\textrm{div}(\mathbf{J}) = 0$. The rest do not.

\item If we want $\textrm{div}(\mathbf{J}) = 0$ then we need
$\partial J_1/\partial x + \partial J_2/\partial y = 0$, or
$$
1+2xy+\frac{\partial f}{\partial y} = 0.
$$
Write this as $\frac{\partial f}{\partial y} = -1-2xy$ and integrate in $y$ to find $f(x,y) = -y-xy^2+c(x)$, where $c(x)$ is
any function of $x$.

\item First, $\nabla \phi = <1,0>$ here. If we dot this with $\mathbf{n}=<-1,0>$ (left side of the domain) we find
    $\nabla\phi\cdot\mathbf{n}=-1$, as required. On the top, right side, and bottom ($\mathbf{n}=<0,1>,<1,0>,<0,-1>$, respectively) we have $\nabla\phi\cdot\mathbf{n}=0$, $\nabla\phi\cdot\mathbf{n}=1$, and $\nabla\phi\cdot\mathbf{n}=0$, as required.

\item With $\phi(x,y)$ as given we find $\nabla \phi = <2x,-2y>$, so
$\gamma \nabla \phi = <2x+2x^2y, -2y-2xy^2>$ and it follows that $\textrm{div}(\gamma\nabla\phi) = 0$.

On the left side $x=0$ we find $\nabla \phi = <0,-2y>$ and $\gamma = 1$, so $\gamma\nabla \phi = <0,-2y>$ and $\gamma\nabla\phi\cdot\mathbf{n} = 0$ since $\mathbf{n}=<-1,0>$.
On the top $y=1$ we find $\nabla \phi = <x,-2>$ and $\gamma = 1+x$, so $\gamma\nabla \phi = <x+x^2,-2-2x>$ and $\gamma\nabla\phi\cdot\mathbf{n} = -2x-2$ since $\mathbf{n}=<0,1>$.
Similar computations on the right side and bottom of $\Omega$ show that this solution satisfies the Neumann boundary data.

\item This is immediate---in the PDE $\textrm{div}(\gamma\nabla \phi) = 0$ the constant $c$ appears under a derivative, as $c$ also does
    in the gradient $\nabla \phi$ in the boundary condition.

\item This is just the Divergence Theorem applied to $\mathbf{J}$ on
$\Omega$:
$$
\int_{\partial \Omega} \mathbf{J}\cdot\mathbf{n}\,ds = \int_{\Omega} \textrm{div}(\mathbf{J})\,dA.
$$
where $dA$ is the area element. But since $\mathbf{J} = -\gamma\nabla \phi$ we find
$$
\int_{\partial \Omega} \gamma\nabla \phi \cdot\mathbf{n}\,ds = \int_{\Omega} \textrm{div}(\gamma\nabla \phi)\,dA.
$$
The integral on the right above is zero (since $\textrm{div}(\gamma\nabla \phi)=0$) and the integral on the left above is ${\displaystyle \int_{\partial \Omega} g\,ds}$.

\item Here $b_4=2$ with all other $b_k$ and all $a_k=0$. From (\ref{cksol})-(\ref{dk1sol}) we find $d_4=\frac{6561}{13120},d_{-4}=\frac{1}{13120}$, and all
other $d_k$ and all $c_k$ equal zero. From (\ref{sol1}) the solution is
$$
\phi(r,\theta) = \frac{6561}{13120}r^4\sin(4\theta) + \frac{1}{13120r^4}\sin(4\theta).
$$

\item In this case we know $g$ has Fourier coefficients $a_3=1$, all other coefficients equal to zero. We compute $c_3=\frac{R^{-3}}{3(R^{-3}-R^3)}$
and $c_{-3}=\frac{R^{3}}{3(R^{-3}-R^3)}$.  The solution is therefore of the form
$$
\phi(r,\theta) = \frac{R^{-3}}{3(R^{-3}-R^3)}r^3\cos(3\theta) + \frac{R^{3}}{3(R^{-3}-R^3)}r^{-3}\cos(3\theta).
$$
When $r=1$ this is, after simplifying
$$
\phi(r,\theta) = \frac{1+R^6}{3(1-R^6)}\cos(3\theta).
$$
Matching to the boundary data $4097\cos(3\theta)/12285$ yields $\frac{1+R^6}{3(1-R^6)}=\frac{4097}{12285}$. The only nonnegative solution is $R=1/4$ (all other solutions are
either real negative or complex).

\item
\begin{enumerate}
\item[a.] In this case we know $g$ has Fourier coefficients $a_m=1$ for the fixed $m$, all other coefficients equal to zero. We compute $c_m=\frac{R_0^{-m}}{m(R_0^{-m}-R_0^m)}$
and $c_{-m}=\frac{R_0^{m}}{m(R_0^{-m}-R_0^m)}$.  The solution is therefore of the form
$$
\phi(r,\theta) = \frac{R_0^{-m}}{m(R_0^{-m}-R_0^m)} r^m\cos(m\theta) + \frac{R_0^{m}}{m(R_0^{-m}-R_0^m)} r^{-m}\cos(m\theta).
$$
When $r=1$ this is, after simplifying
$$
\phi(r,\theta) = \frac{1+R_0^{2m}}{m(1-R_0^{2m})}\cos(m\theta),
$$
or $q_m\cos(m\theta)$ where $q_m$ is as asserted.
Since $0<R_0^{2m}<1$ for $0<R_0<1$ we see that $1+R_0^{2m}>1$ and $1-R_0^{2m}<1$ in this range, so that $\frac{1+R_0^{2m}}{m(1-R_0^{2m})}>1$ and so
$q_m>1/m$.

\item[b.] To solve the inverse problem we seek a value of $R$ so that the solution $\phi(1,\theta)$ agrees with the ``measured'' $q_m\cos(m\theta)$ at the boundary $r=1$, leading to
    $$
\frac{1+R^{2m}}{m(1-R^{2m})} = q_m
    $$
    for unknown $R$. If we make the substitution $T=R^{2m}$ this above displayed equation becomes
$$
\frac{1+T}{m(1-T)} = q_m.
$$
When $T\rightarrow 0^+$ the left side above approaches $1/m$, and so dips below $q_m$. As $T\rightarrow 1^-$ the left side approaches $\infty$. Moreover,
$\frac{d}{dT}((1+T)/(1-T)) = 2/(1-T)^2>0$, so the left side above is strictly increasing. We conclude from the Intermediate Value Theorem that $\frac{1+T}{m(1-T)} = q_m$
has a solution in $(0,1)$, and that the solution is unique. By similar reasoning, the equation $R^{2m}=T$ has a unique solution $R\in (0,1)$ (since $R^{2m}$ ranges from $0$ to $1$, is continuous, and strictly increasing). Of course, $R=R_0$ clearly works, and so must be the unique (and correct) solution.
\end{enumerate}

\item If $g$ is not zero them simply choose a nonzero Fourier coefficient $a_m$ or $b_m$. For simplicity, suppose $a_m\neq 0$, the coefficient of the $\cos(m\theta)$ term
in $g$ (which we can compute). From (\ref{sol1}) we know that (with $c_0=0$)
$$
\phi(1,\theta) =\sum_{k=1}^{\infty} \left ((c_k + c_{-k})\cos(k\theta) + (d_k + d_{-k})\sin(k\theta)\right ).
$$
We can pick out the value of $c_m+c_{-m}$ (the $\cos(m\theta)$ Fourier cosine coefficient of $\phi(1,\theta)$) by integrating $\phi(1,\theta)\cos(m\theta)$ for
$\theta=-\pi$ to $\pi$. Thus we know $a_m$ from the input flux $g$, and we know $c_m+c_{-m}$ from the boundary potential data, so by using (\ref{cksol})-(\ref{ck1sol})
we arrive right back at
$$
\frac{1+R^{2m}}{m(1-R^{2m})} = q_m
$$
with $q_m = \frac{1+R_0^{2m}}{m(1-R_0^{2m})}$ as in Problem \ref{prob10}, and so we can uniquely determined $R_0$.

\end{enumerate}

}

%\section*{ACKNOWLEDGEMENTS}


\begin{thebibliography}{99}

\bibitem{bauman} Bauman, P. Geophysics Blog. \href{http://www.paulbaumangeophysics.com/blog/previous/6}{http://www.paulbaumangeophysics.com/blog/previous/6}. Accessed 15 August 2019.

\bibitem{cheney} Cheney M., D.~Isaacson,  and J.~C.~Newell. 1999. Electrical impedance tomography. \textit{SIAM Review}  41(1):  85-101.

\bibitem{loke} Loke, M. and R.~Barker. 1996. Practical techniques for 3D resistivity surveys and data inversion. \textit{Geophysical Prospecting}. 44: 499-523.

\bibitem{McClymont} McClymont, A., R.~Bauman, R.~Freund, J.~Seligman, J., H.~Jol, K.~Bensimon, K. and P.~Reeder. 2017. The Forgotten Holocaust in Lithuania - Geophysical Investigations at the Ponary Extermination Site. In the \textit{Proceedings of the 23$^{\rm rd}$ European Meeting of Environmental and Engineering Geophysics}, Malm\"{o}. Sweden. 3-7 September.

\bibitem{nova} Public Broadcasting Service NOVA Special ``Holocaust Escape Tunnel.'' 19 April  2017.  \href{https://www.pbs.org/wgbh/nova/video/holocaust-escape-tunnel}{https://www.pbs.org/wgbh/nova/video/holocaust-escape-tunnel},
also freely available at \href{https://www.dailymotion.com/video/x5j57yg}{https://www.dailymotion.com/video/x5j57yg}.

\bibitem{sakowicz} Sakowicz, K. and Y.~Arad.  2005. \textit{Ponary Diary, 1941-1943: A Bystander's Account of a Mass Murder}. (Translated by Laurence Weinbaum).New Haven: Yale University Press.

\bibitem{shaer} Shaer, M. 2017. The Holocaust's Great Escape,  \textit{Smithsonian Magazine}. 47(10):  42-53.

\end{thebibliography}

\end{document}
