\documentclass[reqno]{amsart}
\usepackage{hyperref}

\AtBeginDocument{{\noindent\small
{\em Electronic Journal of Differential Equations},
Vol. 2005(2005), No. 132, pp. 1--12.\newline
ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu
\newline ftp ejde.math.txstate.edu  (login: ftp)}
\thanks{\copyright 2005 Texas State University - San Marcos.}
\vspace{9mm}}

\begin{document}
\title[\hfilneg EJDE-2005/132\hfil Level set method]
{Level set method for solving Poisson's equation with 
discontinuous nonlinearities}
\author[J. Kolibal\hfil EJDE-2005/132\hfilneg]
{Joseph Kolibal}

\address{Joseph Kolibal\hfill\break
 Department of Mathematics \\
University of Southern Mississippi \\
Hattiesburg, MS 39406, USA}
 \email{Joseph.Kolibal@usm.edu}

\date{}
\thanks{Submitted October 14, 2005. Published November 25, 2005.}
\subjclass[2000]{35J05, 35J60}
\keywords{Laplace equation; reduced wave equation (Helmholtz);
\hfill\break\indent Poisson equation;  nonlinear elliptic PDE}

\begin{abstract}
 We study semi-linear elliptic free boundary problems 
 in which the forcing term is discontinuous;
 i.e., a Poisson's equation where the  forcing term is
 the Heaviside function applied to the unknown minus a constant.
 This approach uses level sets to construct a monotonic
 sequence of iterates which converge to a class of solutions to the
 free boundary problem. The monotonicity of the construction
 based on nested sets provides insight into the connectivity of
 the free boundary sets associated with the solutions.
\end{abstract}

\maketitle

\numberwithin{equation}{section}
\newtheorem{theorem}{Theorem}[section]
\newtheorem{lemma}[theorem]{Lemma}
\newtheorem{corollary}[theorem]{Corollary}


\section{Introduction}
\label{sec:introduction}

We are interested in computing solutions to the
semi-linear elliptic boundary-value problem
\begin{equation}
\begin{gathered}
   - \Delta \psi  =  \lambda H(\psi-1) \quad \mbox{in }\Omega, \\
         \psi  =  0  \quad\mbox{on }\partial \Omega,
\end{gathered}  \label{eq:basic}
\end{equation}
on bounded domains $\Omega \subset \mathbb{R}^n$,  $1 \leq n \leq 3$,
where $H$ is the Heaviside function,
$H(t) = 0$ for $t  \leq  0$ and $H(t) = 1$ for $t  >  0$,
and $\lambda > 0$ is a scalar parameter.
The semi-linear problem may be reformulated as
the equivalent free boundary problem, find $\psi$ such that
\begin{equation}  \label{eq:char}
\begin{gathered}
 -\Delta \psi = \begin{cases}
    \lambda & \mbox{on }A,    \\
         0 & \mbox{on }\Omega \setminus A,
\end{cases}\\
\psi = 0 \quad \mbox{on } \partial \Omega
\end{gathered}
\end{equation}
where  $A=\{ x \in \Omega : \psi(x)  \geq  1 \}$ is
the free boundary set
and $\partial A$ is the free boundary which is to be determined.
Since we may interpret the forcing term of (\ref{eq:char}) as $\lambda \chi_A$,
where $\chi_A$ is the characteristic function of the set $A$,
we need only insist that the solution be as regular as the solution
of Poisson's equation, $-\Delta \psi = \lambda \chi_A$ in $\Omega$
with $\psi|_{\partial \Omega} = 0$
in order to completely specify $\psi$.

Equations of this type
arise in the consideration of vortex flow,
porous medium combustion problems and
plasma dynamics \cite{friedman,kan,turkington,elcrat}.
In particular, problems of the form $-\Delta \psi = f(x,\psi) \geq 0$ where
the forcing
term $f$ is H\"{o}lder continuous with exponent $\alpha$, $0 < \alpha
< 1$, with suitable restrictions on the growth of $f$
have been shown to have solution $\psi \in C^{2+\alpha}$
using isoperimetric variational methods \cite{baf}.
In a related approach,
Fraenkel and Berger \cite{fab} in $\mathbb{R}^3$,
and Norbury \cite{norbury} in $\mathbb{R}^2$,
have proven the existence of a solution of the semi-linear elliptic
differential equation by maximizing certain functionals on the surface of
a sphere in a Sobolev space.
By taking the limit of H\"{o}lder continuous functions and
using Steiner symmetrization \cite{steiner,berger} and the generalized
maximum principle of Littman \cite{littman}, Fraenkel and Berger \cite{fab}
are able to construct a solution $\psi \in C^{1 + \alpha}$ when
the forcing term $f$ is discontinuous.
Keady \cite{keady,keadykloeden} has dealt explicitly with (\ref{eq:basic});
 however, the interest was in obtaining estimates for doubly symmetrized
domains in $\mathbb{R}^2$.

The primary reason for considering the semi-linear form $(\ref{eq:basic})$ of
the problem
is the ease with which variational methods \cite{berger,friedman,turkington}
yield results on the existence of solutions
with the computational aspects of the problem being reduced to
considering an isoperimetric problem in the calculus of variations.
The approach we take
to finding nontrivial solutions of (\ref{eq:basic})
is based on considering the
equivalent free boundary problem (\ref{eq:char}).
Although a direct attack on the semi-linear problem is appealing,
the alternative free boundary formulation
can be utilized to construct an efficient solution algorithm,
and we consider our numerical scheme to be an improvement on that of Alexander
and Fleishman \cite{aaf}.
Each stage of the construction is based on solving Poisson's equation
\begin{equation}
\begin{gathered}
 -\Delta u_{i+1}  =  \lambda  \chi_{A_i}  \quad\mbox{in }\Omega, \\
         u_{i+1}  =  0  \quad\mbox{ on }\partial \Omega,
\end{gathered}  \label{eq:pois}
\end{equation}
where the sets $A_i$ are iteratively corrected so that
$A_i  \to A$, without searching for the
boundary of $A$, that is,  where  $u(x) = 1$ in $\Omega$.
The approach exploits the linearity of
(\ref{eq:pois})
along with the maximum principle
and the scalability of (\ref{eq:basic})
to construct an algorithm which converges
monotonically to the desired free boundary set.
By scalability we mean that if $k > 0$
and if $\upsilon = k \psi$ and $\mu = \lambda k$,
then $- \Delta \upsilon = \mu H(\upsilon - k)$
whenever $-\Delta \psi = \lambda H(\psi -1)$ in $\Omega$
with $\psi|_{\partial \Omega}     = \upsilon|_{\partial \Omega} = 0$.
Hence the rescaled solution satisfies a closely related problem.
The value of $\upsilon$ on the boundary of the free boundary set
which is implicit in the problem is now $k$ instead of $1$.
This algebraic result will be important
in formulating the solution algorithm.
To avoid confusion
with solving the related equation $-\Delta  \psi = \lambda H(\psi - k)$,
we consider only the canonical case when $k=1$.


To illustrate some key aspects of this problem,
consider an example on the unit disk $\subset \mathbb{R}^2$.
By symmetry the solution of (\ref{eq:char})
satisfies an ordinary differential equation in the radial variable $r$,
\begin{equation}
\begin{gathered}
- \frac{1}{r} \frac{d}{dr}  \big(r  \frac{d \psi }{dr}   \big)
=  \lambda H(\psi-1) \quad \mbox{in }r \leq 1,    \\
       \psi  =  0 \quad \mbox{at }r=1\,.
\end{gathered} \label{eq:line}
\end{equation}
The solution to (\ref{eq:line}) is completely determined by imposing
$C^1$ continuity of $\psi$ across the free boundary at $r=a$, thus
\begin{equation}
\psi(r) = \begin{cases}
 1 + \lambda(a^2 - r^2) / 4,     & 0 \leq r \leq a, \\
 \log (r) / \log (a), & a \leq r \leq 1,
\end{cases}  \label{eq:sdsk}
\end{equation}
and it can be verified easily that $\lambda=-2/a^2 \log (a)$.
Since for all $a > 0$, $\psi(0) = 1 + \lambda a^2/4 > 0$,
it might appear that radial solutions to (\ref{eq:sdsk}) exist
for  all $\lambda > 0$; however, since
$\lambda$ and $a$ are not free parameters,
a unique nontrivial solution pair $(\psi; \lambda)$ may not exist
for all choices of $\lambda$.
For solutions of (\ref{eq:line})
there are two values of $a$ which
will yield exactly the same solution,
corresponding to each branch
of the solution curve.
Only the trivial solution exists if $\lambda < \lambda_{min} =  4 \mbox{e}$.
Moreover, any algorithm which relies on solving (\ref{eq:char})
by successively solving a sequence of related linear problem in which the
free boundary set is fixed at each iteration
and then searching for the free boundary on the next iteration,
must avoid having
the estimate for the free boundary set collapse into the trivial solution.
This is precluded from happening in the algorithm we present
because the construction is monotone increasing and solutions
do not depend on the choice of $\lambda$.
We never attempt to compute solutions for values of $\lambda$
which are below the minimum for which solution
will collapse into the trivial solution.
In solving this example problem in (\ref{eq:line})
the approach relied on choosing the value of $a$,
or more generally the free boundary set as the the free parameter,
and then computing $\lambda$.
In much the same way, in developing a solution algorithm
on arbitrary domains in $\mathbb{R}^2$, or $\mathbb{R}^3$, it will be easier
working with the free boundary set,
computing the values of  $\lambda$ after having obtained a solution.
Thus, while $\lambda$ may be taken as a free parameter in (\ref{eq:basic}),
we prefer to work with (\ref{eq:char}) and the
free boundary set, requiring instead that the value
of $\lambda$ be determined after having found a free boundary set $A$.

The numerical algorithm we describe
for solving this problem is an easy and direct implementation
of the constructive
proof of the existence of solutions to (\ref{eq:char}).
In Section \ref{sec:sequence} we describe in detail the construction
of the sequences of sets $A_i$
which in Section \ref{sec:existence} are shown to converge
to the set free boundary $A$ on which the solution of
(\ref{eq:basic}) depends.
Since construction of the free boundary
set $A$ relies on solving an elementary linear equation (\ref{eq:pois}) at
each step of the construction,
numerical results are easily obtained using a standard linear
finite element formulation.

\section{A monotone algorithm for constructing the sets $A_i$}
\label{sec:sequence}

We construct a sequence of nested sets,
$A_0 \subseteq A_1 \subseteq \dots \subseteq A_n$, which
in Section \ref{sec:existence} will be
shown to converge to the free boundary set $A$.
Presently our concern is only with the specifics of the
construction, which is done recursively, and some immediate
consequences, notably the monotonicity of the
sequence $\{c_n\}$ of numbers associated with these subsets.

Let $u \in C^0(\Omega)$, $u > 0$ in $int(\Omega)$,
and consider a  subset of the graph of a function $u$,
defined by selecting the level set
\begin{equation}
\Sigma(c) = \{(u,x) : u(x) \geq  c, x \in \Omega \}  ,
\end{equation}
where $0 < c < \max u$.
We are interested in the projection of $\Sigma(c)$ onto $\Omega$,
defined by $A = \{ x \in \Omega \mid (u,x) \in \Sigma(c) \}$.
The recursive construction of the nested sets $\{A_i\}$ is based
on beginning with an initial solution to the boundary-value problem
\begin{equation}
\begin{gathered}
 -\Delta u_0   =  \lambda  \chi_B     \quad\mbox{in }\Omega,  \\
         u_0   =  0      \quad\mbox{on }\partial \Omega,
\end{gathered}  \label{eq:de}
\end{equation}
where, as before, $\lambda >  0$ and $B  \subseteq    \Omega$.
The special case when $B  =  \Omega$ is examined initially
(and subsequently extended to more general solutions in which $B$ is any open
subset of $\Omega$).
We shall refer to $u_0$ obtained by considering $B  =   \Omega$ as the basic
solution to our problem, partly because on some elementary domains the
free boundary set $A$ is easily generated
or derived from this fundamental initial choice.

Thus the basic solution is initiated by solving
(the elastic torsion problem),
\begin{equation}
\begin{gathered}
-\Delta u_0  = \lambda  \quad \mbox{in }\Omega,\\
u_0\big|_{\partial \Omega}  = 0.
\end{gathered}\label{eq:it1}
\end{equation}
Fix an arbitrary $c_0 \in (0,\max u_0)$ to be the starting level and
denote the projection of $\Sigma(c_0)$ onto $\Omega$ as
\begin{equation}
A_0 =  \{x \in \Omega : u_0(x)  \geq  c_0 \}.
 \label{eq:adef}
\end{equation}
Continuing recursively, for $i \geq 1$ we denote by $u_i$ the solution of
\begin{equation}
 -\Delta u_i = \lambda \chi_{A_{i-1}} \quad \mbox{in }\Omega
\quad\mbox{and}\quad
 u_i\big|_{\partial \Omega}  = 0,
  \label{eq:uset}
\end{equation}
and determine
\begin{equation}
c_i  =  \min_{x \in A_{i-1}}  u_i(x).
 \label{eq:cdef}
\end{equation}
By projecting $\Sigma  (c_i)$ onto $\Omega$, the set
\begin{equation}
A_i  =  \{x \in \Omega : u_i(x)  \geq  c_i \}
 \label{eq:it2}
\end{equation}
is obtained, and the motivation for the algorithm is now apparent.

Let $A_n$ be the estimate for the free boundary set
on the $n$-th iteration, then solve (\ref{eq:pois}) assuming
that $A_n$ solves (\ref{eq:char}).
If it does, then the solution on the boundary of $A_n$
equals some constant, say $k$.
By rescaling the solution, we can make it equal to one on the boundary of
$A_n$ and we are done.
Otherwise, find the minimum of this new solution over $A_n$
and extend $A_n$ to $A_{n+1}$, constructed
so that it includes all those point where this new solution is
greater than $k$.
Then test whether $A_{n+1}$ is a solution, continuing in this
fashion until a sufficiently converged estimate of the
free boundary set $A$ is attained.

The choice of $B = \Omega$ is the largest set which
solves the linear problem (\ref{eq:pois}),
and the extent to which (\ref{eq:basic}) is satisfied,
or nearly satisfied by
any particular set $A_0$ generated by projecting $\Sigma(c_0)$
on the first iterative step
depends on the extent to which the boundary of the free boundary
set is similar to the boundary of $\Omega$.
For example, for symmetric solutions on the unit disk $D$,
solving $-\Delta u_0 = \lambda$ in $D$ gives
the set $A_0 =\{x \in D \mid x \cdot x \leq a \}$ where
$0 < a < 1$ is determined only by the choice of $0  <  c_0 < \|u_0\|_{\infty}$.
In this example, each solution to (\ref{eq:basic}) is a disk contained in $D$
 which is obtained immediately on the first projection to $A_0$.
On more complex domains,
we will find
that $\partial A$ corresponds, in a sense,
to a homotopy between the initial
guess $\partial A_0$ and $\partial \Omega$.
Indeed, the number of steps taken to arrive at a solution will in part
depend on the number of projections needed to achieve this
deformation by iterative projection.

If the $u_i$ are continuous on the domain and the domain is sufficiently
regular, this procedure produces a monotonic sequence of numbers $\{c_i\}$
and an associated sequence of nested subsets $\{A_i\} \subseteq \Omega$.


\section{Preliminaries}
\label{sec:existence}

In proving the existence of solutions to the semi-linear equation,
we take $\Omega$ sufficiently smooth and regular
so as to apply standard results
from the theory of partial differential
equations \cite{adams,showalter,schecter}.
We take $\Omega$ bounded, obeying the cone condition \cite{showalter}.
We also require that $\partial \Omega$ is either $C^2$ or that $\Omega$
is convex polyhedral \cite{grisvard}.

As usual, we introduce the Sobolev spaces $H^m$ and $H^m_0$, defined as the
closure of $C^{\infty}$ and $C^{\infty}_0$, respectively,  in the norm
\begin{equation}
\| u \|_{H^m} = \Big\{ \sum_{| k | \leq m}
 \| D^k u \|^2_{L_2} \Big\} ^{1/2},
\end{equation}
where $k=(k_1,k_2,\dots,k_n)$ is an $n$-tuple of non-negative
integers, and
where $D^k = D^{k_1}_1 D^{k_2}_2 \dots D^{k_n}_n$ is the partial
derivative of order $\mid k \mid =k_1 + k_2 + \dots + k_n$.
We define the bilinear form $(u,v)$ corresponding to the Laplace operator as
\begin{equation}
( u,v ) = \int_{\Omega} \nabla u \cdot \nabla v  \,  d \Omega,
\end{equation}
It will also be convenient to introduce the inner
product of two functions in $L_2$,
\begin{equation}
 \langle u,v \rangle  =  \int_{\Omega} u v \, d\Omega  .
\end{equation}
Throughout, we take $x \in \mathbb{R}^n$, and denote the Lebesgue
measure of the set $A$ by $m(A)$.
By a solution $u$ of the linear partial differential equation we mean $u$
satisfies
\begin{equation}
     (u,v)  =  \lambda \langle \chi_A, v \rangle,
   \qquad \forall v \in H_0^1(\Omega);
\end{equation}
in other words, $u$ is a weak solution of the
differential equation, although we will continue to
pose the problem classically using $-\Delta u = \lambda_A$.
We remark that $u \equiv 0$ for all $\lambda$ is the trivial solution
and henceforth
solution means nontrivial, weak solution.

\begin{lemma} \label{thm:slb}
Consider the construction
using \eqref{eq:uset}--\eqref{eq:it2}.
Let $-\Delta u_i=\lambda \chi_{A_{i-1}}$
and $A_i \subset \Omega$, $u_i\big|_{\partial \Omega}=0$.
Then $u_i \in H^2(\Omega) \cap H^1_0(\Omega)$. Also, $u_i \in C^0 (\Omega)$ if
$\Omega \subset \mathbb{R}^n$, $n \leq 3$.
\end{lemma}

\begin{proof}
This is a straightforward consequence of
elliptic regularity and the Sobolev Imbedding Theorem since
$\chi_{A_i} \in L_2(\Omega)$ \cite{adams,showalter,hutsonpym}.
\end{proof}


See \cite[Chapter~11.6, Theorem 11.6.6]{hutsonpym} for a  discussion and in
particular \cite[Chapter 7]{gilbargandtrudinger} which shows that we could
improve the results of this lemma (for more general $\Omega$,
or for more general elliptic operators, or with regard to the
smoothness of $u_i$), however we have chosen
to exploit only the continuity of each $u_i$.
Although the solution is smoother than this on either side
of the boundary of $A_i$ and for that matter across $A_i$,
continuity of the solutions is entirely sufficient
to demonstrate the construction.
Note that for $A_i  \subset  \Omega$, $\chi_{A_i}$ will vanish in a
neighborhood of $\partial \Omega$ so that we are only considering harmonic
functions outside of $A_i$ near $\partial \Omega$.
We remark that in the construction defined in Section \ref{sec:sequence},
since we consider weak solutions of the
differential equation, the inequality defining the subset
of the graph of $u_i$ given by $\Sigma(c_i)$ could have excluded the equality
without effect, with $A_i$ being defined using open sets and the infimum
replacing the minimum.
In particular, this construction and the continuity of $u$
results in sets $A_i$ which
are nice in the sense that $int A_i$ is open, being the projection
of $int \Sigma(c_i)$ onto $\Omega$.

We are now ready to formally demonstrate
that the construction described in Section \ref{sec:sequence}
produces an algorithm which is convergent to a solution of (\ref{eq:basic}).
We begin by showing that on each iteration,
when we examine the graph of the solution
$u_i$ restricted to the set $A_i$, that no minima can exist
in the interior of $A_i$.
Otherwise the algorithm would fail since the level sets $\Sigma(c_i)$
used in constructing the next support set $A_{i+1}$
are based on finding $c_i$, the minimum
of $u_i$ over the set $\partial A_i$.
Algorithmically, this avoids the examination
of $u_i$ over all of $A_i$; equally,
it shows that the free boundary set cannot develop
on any open set which is interior to the
current estimator $A_i$ of $A$.

\begin{lemma} \label{thm:awmax}
Consider $A_i  \subset  \Omega$ constructed
using \eqref{eq:uset}--\eqref{eq:it2}
with $-\Delta u_i= \chi_{A_{i-1}}$ in $\Omega$,
$u_i\big|_{\partial \Omega}=0$.
Then the minimum of $u_i$ over $A_{i-1}$ occurs on the boundary of $A_{i-1}$.
\end{lemma}

\begin{proof}
 From Lemma~\ref{thm:slb}, $u_i \in  H^2(\Omega)  \cap  H^1_0(\Omega)$
and is continuous on $\Omega$.
Define $w_i$ to be the restriction of $u_i$ to $A_{i-1}$.
The solution $w_i$ of $-\Delta w_i =1$
on $A_{i-1}$, $w\big|_{\partial A_{i-1}} = u_i\big|_{\partial A_{i-1}}$
We apply the weak maximum principle to $w_i$
along with $w_i(x) = u_i(x)|_{A_{i-1}}, \forall x \in A_{i-1}$, to
conclude that
\[
 \min_{\partial A} u_i  \leq  \min_{A_{i-1}}        u_i .
\]
 From the continuity of $u_i$, equality can be achieved at most on
a set of measure zero.
Otherwise there is an open ball, inside a set of positive measure
in $A_{i-1}$, on which $\Delta  u_i  =  0$ by direct differentiation.
This contradicts $-\Delta  u_i  =  1$ on $A_{i-1}$.
\end{proof}

Lemma \ref{thm:awmax} demonstrates
the equivalence between the semi-linear and free boundary formulation.
Since the minimum occurs on the boundary of $A_{i-1}$, if $u_i = 1$ on $\partial
 A_{i-1}$,
then $u_i$ solves (\ref{eq:char}) with $u_i > 1$ in the interior of $A_{i-1}$,
consequently $u_i$ solves (\ref{eq:basic}).
Observe, that since we construct $A_i$ from the level set of $u_i$, the
existence of a minimum (on a set of measure 0) in the interior equal to that
which occurs on the
boundary of $A_i$ also has no effect on the construction.


\begin{lemma} \label{thm:mos}
Consider the construction
\eqref{eq:it1}--\eqref{eq:it2}
for $i \geq 1$, if $u_{i-1} \in H^2(\Omega) \cap C^0(\Omega)$,
then $A_{i-1} \subseteq A_i$, $u_i  \leq  u_{i+1}  \leq  u_0$, and
$c_i \leq c_{i+1}  \leq  c_0$.
\end{lemma}

\begin{proof}
The set inclusion follows since on the set $A_{i-1}$, $u_i$ has minimum
value
$c_i$ while the set $A_i$ includes all those values of $x$
for which $u_i \ge c_i$.

We show next that $c_0 \geq c_i$.  For functions in $H^2(\Omega)$
the necessary generalization of the usual maximum principle to
precisely the weak solutions that we require is discussed in Gilbarg and
Trudinger \cite[Chapter 8]{gilbargandtrudinger}.
We have constructed $u_0$ such that the support of the forcing
term is $\Omega$, $A_0$ is the projection of $\Sigma(c_0)$
into $\Omega$, and $u_i$ is the function
which solves the differential equation where the support of the forcing
term is $A_{i-1}$.
Hence the maximum principle
and the inclusion of $A_0$ in $A_{i-1}$
may be used to find
\begin{align*}
A_{i-1} \subseteq \Omega
   & \Rightarrow  \chi_{A_{i-1}} \leq \chi_{\Omega}  \\
   &  \Rightarrow    u_i(x) \leq u_0(x)   \\
   &  \Rightarrow    \min_{x\in A_{i-1}} u_i(x)
                       \leq \min_{x\in A_{i-1}} u_0(x)
                       \leq \min_{x\in A_0} u_0(x)  .
\end{align*}
By the definition \eqref{eq:cdef} of $c_i$ and \eqref{eq:adef}
 of $A_0$ we have
\[
 c_i \leq c_0 \,.
\]
For $i \geq 1$, we recall constructing
$A_i$ as the set of all $x$ such that $u_i(x)  >  c_i$,
and $u_{i+1}$ as the solution of the
differential equation where the support of the forcing term is $A_i$.
Hence, again using the maximum principle
\begin{align*}
A_{i-1} \subseteq  A_i & \Rightarrow \chi_{A_i} \geq \chi_{A_{i-1}} \\
    &  \Rightarrow    u_{i+1}(x) \geq u_i(x)  \\
    &  \Rightarrow     \min_{x\in A_i} u_{i+1}(x)
                         \geq \min_{x\in A_i} u_i(x)  .
\end{align*}
Since $c_{i+1}$ is the minimum value taken by the function $u_{i+1}$ on $A_i$,
while $A_i$ is the set where $u_i  >  c_i$, we have $c_{i+1} \geq c_i$.
\end{proof}

An obvious special case occurs if $A_i =  A_{i-1}$ for some $i$.
Then an easy argument
shows that $c_i = c_{i+1}$ (and thus all subsequent iterations $A_{i+n}$
 coincide with $A_i$, with $c_{i+n}  =  c_i$ for all $n$)
and we have found a set $A$ which solves the free boundary problem.
Indeed, we have two cases to consider:
either we trivially find a solution for some fixed $i$ and there is no growth
in the sets $A_i$, or there is growth and $A_{i-1} \subset  A_i$ strictly,
for every $i$.


The next result affirms that if the sets $A_i$ are not all identically
the free boundary set,
then the sequence, $\{c_i\}$, $i = 1, 2, \dots$ induced by the iterative
construction is strictly monotone increasing, and consequently
the functions $u_{i}$ and $u_{i+1}$ solving \eqref{eq:uset}
are strictly monotone increasing, i.e., $u_{i+1} > u_i$ in $int(\Omega)$.
Specifically, we wish to show that in the non-trivial case,
if each solution $u_{i}$ is continuous, then
the constructive algorithm
discussed in \eqref{eq:uset}--\eqref{eq:it2}
provides the strict monotonicity required for convergence.

\begin{lemma} \label{thm:maxprin}
Let $u_i$ and $u_{i+1}$ be continuous solutions
of $-\Delta u_{i}=\chi_{A_{i-1}}$
in $\Omega$ based on the construction in
\eqref{eq:it1}--\eqref{eq:it2} with
$A_{i-1} \subset A_{i} \subset \Omega$, and $u_i$,
$u_{i+1}\big|_{\partial\Omega}=0$,
then $u_{i+1} > u_{i}$, for all $x$ in  interior or $(\Omega)$.
\end{lemma}

\begin{proof}
Let $v = u_{i+1} - u_{i}$ with $v\big|_{\partial \Omega} = 0$.
Since $v$ is continuous, there is an open disk $D$ contained in
$A_{i} \setminus A_{i-1}$.
Construct
$f(x) \in C^{\infty}(\Omega)$
with $f = 1$ in $int(D)$ and $f = 0$
on $\Omega  \setminus  D$.
Then $\chi_{A_{i} \setminus A_{i-1}} \geq f$
on $\Omega$. Since $f \in L^2(\Omega)$, by comparing the solution of
$-\Delta v = \chi_{A={i} \setminus A_{i-1}}$, ${v}\big|_{\partial \Omega}$,
with $-\Delta \hat{v} = f$, $\hat{v}\big|_{\partial \Omega} = 0$, using the
weak maximum principle for weak solutions and the continuity of $v$ and
$\hat{v}$, from Lemma~\ref{thm:mos} we obtain that $v \geq \hat{v}$ on $\Omega$.
Now since $\hat{v} \in C^2(\Omega) \cap C(\Omega)$,
we obtain by the strong maximum principle of classical solutions of Poisson's
equation that $\hat{v} > 0$ in $\Omega$.
Hence $u_{i+1} > u_{i}$ in $int(\Omega)$.
\end{proof}

It remains to show that
the sets $\{ A_i \}$ will always grow boundedly to
the set $A$ which solves the free boundary problem.

In the construction of the sets $A_i$
in Section \ref{sec:sequence}, $A_{i-1}\subseteq A_i$,
however the sets may have the same measure, i.e., $m(A_i \setminus A_{i-1})=0$.
This would be insufficient to assure convergence,
except in the special case when $A_i$ corresponds to a free boundary set,
when $A_i = A_{i-1}$ for all $i$.
However,
in the nontrivial case the construction described in Section~\ref{sec:sequence},
along with the continuity of the $u_i$, is sufficient
to assure that $m(A_i \setminus A_{i-1}) > 0$. Then we have

\begin{lemma} \label{thm:conw}
Consider the construction in \eqref{eq:it1}--\eqref{eq:it2}
with $A_{i-1}  \subseteq  A_i$,
$c_i  = \min \{u_i(x) \mid x \in \partial A_{ i-1} \}$,
and $-\Delta u_i = \chi_{A_{i-1}} \geq 0$ in $\Omega$.
If $d_i = \max \{u_i(x) \mid x \in \partial A_{i-1} \} >  c_i  $,
then $m(   A_i         \setminus  A_{i-1}) > 0$.
\end{lemma}

\begin{proof}
Since $u_i \in C^0  (\Omega)$, we have that
$int\/A_{i} = \{x \in \Omega | u_i(x) > c_i \}$,
and similarly $int A_{i-1}$, are open.
Since $u_i(x) > c_i$  for all $x \in intA_{i-1}$
and since
$c_i  <  u_i(x) < d_i$ for some $x \in A_{i}$,
$\exists y \in   A_{i}$ such that $y \notin  A_{i-1}$.
Thus for some $\epsilon > 0$, $\exists D(y,\epsilon)$,
an open ball around $y$,
such that $D(y,\epsilon) \cap  A_{i-1} = \emptyset
\Rightarrow m(  A_i     \setminus  A_{i-1}) \geq m(D(y,\epsilon)) > 0$.
 \end{proof}

We introduce formally now the re-scaling property of the forcing term,
$\lambda H(u-k)$ which is central to the development of this algorithm.

\begin{lemma} \label{thm:cwmax}
Let $A  \subset  \Omega$ and $u$ be a solution
of $-\Delta u= \chi_A$ in $\Omega$, $u\big|_{\partial \Omega}=0$.
If $u$ is constant on $\partial A$, then $u$ solves (\ref{eq:basic}).
\end{lemma}

\begin{proof}
If $u$ is constant on $\partial A$, say $k >  0$, then
$\partial A=\{ x \in \Omega \mid u(x)  =   k   \}$,
consequently $-\Delta u  =  \lambda H(u - k)$.
If we write $k \upsilon =  u$ and $k \mu = \lambda$,
then $-k \Delta \upsilon = k \mu H(k \upsilon - k)$ in $\Omega$
and $k \upsilon\big|_{\partial \Omega} = 0$.  Thus
\[
-\Delta \upsilon   = \mu  H(k(\upsilon -1))  = \mu H(\upsilon - 1)
\quad   \mbox{and}  \quad  \upsilon|_{\partial \Omega}  =  0.
\]
Hence the re-scaled solution $\upsilon$ satisfies (\ref{eq:basic})
as desired.
\end{proof}

If $(u; \lambda)$ is a solution such that $u\big|_{\partial A} = k$,
then $(\upsilon; \mu)$ is a
solution such that $\upsilon\big|_{\partial A} = 1$.
With these ideas in place, we prove that the sequence
of nested subsets
converges to give a solution of the free boundary problem.

\begin{theorem} \label{thm:slw}
Given $\Omega \subset  \mathbb{R}^n$, $1 \leq n \leq 3$, $\lambda > 0$, and
the sequence of nested subsets, $A_0 \subseteq A_1 \subseteq A_2 \dots
\subseteq A_i \dots \subseteq \Omega$
constructed in Section \ref{sec:sequence},
the sequence of solutions of $-\Delta u_i = \lambda \chi_{A_{i-1}}$
converges to a solution u of $-\Delta u = \lambda \chi_A$, where
$A = \cup A_i$, and where $u$, $u_i\big|_{\partial \Omega} = 0$.
Furthermore, $u$ is constant on the boundary of $A$.
\end{theorem}

\begin{proof}
We recall the construction in Section \ref{sec:sequence}
and Lemma~\ref{thm:mos}.
We have a sequence of solutions $\{u_i\}$ of the differential
equation, the support of each forcing term is the
set $A_{i-1}$, and the minimum value of $u_i$ on $A_{i-1}$
is $c_i$.  The support of the forcing terms is increasing and
the supports form a nested sequence of sets in $\Omega$.
The $c_i$, $i \geq 1$,
form an increasing sequence bounded from above by
$c_0$, and hence a convergent sequence.
Similarly, $m(A_i)$ is increasing and bounded
from above by $m(\Omega)$, and hence these form a convergent sequence.

Define $A \subseteq \Omega$ and $c \in [c_1,c_0]$ by
\[
  A  =  \cup_{i=0}^{\infty} A_i,              \quad
  c  =  \lim_{i \to \infty} c_i.
\]
Denote by $u$ the solution of the partial differential
equation $- \Delta u = \lambda \chi_A$, and $u\big|_{\partial \Omega} = 0$.
Since the projection of
the open subset $int\Sigma( c_i )$ of the graph of $u_i$ is open,
each set $nt A_i$ is open, whence
we have that each $\chi_{A_i} \in L_2(\Omega)$
and $\chi_A \in L_2(\Omega)$.  This follows from the construction,
since $int A_i = \{ x \in \Omega \mid u_i(x) > c_i \}$
and each $u_i$ is continuous by Lemma \ref{thm:slb}
since $\chi_{A_{i-1}} \in L_2(\Omega)$.

Let $v_i = u - u_i$. Since $\Delta v_i \in L_2$,
we have using elliptic regularity and the Sobolev Imbedding Theorem that
${\| v_i \|}_{H^2} \leq C_1{\| \chi_A - \chi_{A_i} \|}_{L_2}$
and ${\| v_i \|}_c \leq C_2{\| v_i \|}_{H^2}$,
where $C_1$ and $C_2$ are constants depending only on $\Omega$
and where ${\| u \|}_c$ is the maximum norm.
Since ${\| \chi_A - \chi_{A_i} \|}_{L_2}$ can
be made arbitrarily small for $i$ sufficiently large,
it follows that $u_i \to u$ uniformly on $\Omega$.
Since $u$ is continuous and $c \neq 0$, it follows that $A \neq \Omega$.
Now if $u$ does not equal $c$ everywhere on $\partial A$, then
by Lemmas~\ref{thm:maxprin} and \ref{thm:conw}
we can extend $A$ by projection onto $\Omega$, and obtain a
solution,
which is everywhere in $\Omega$ strictly larger than the solution
we have obtained for $- \Delta  u  =  \chi_A$,
contradicting the fact that $c$ is a limit point of the sequence.
Hence $u  = c$ everywhere on $\partial A$.
\end{proof}

Again, if at any $i \geq 1$ we find that $A_i = A_{i-1}$,
then we have that $u_i = c_i$ on $\partial A_{i-1}$.
Because of this, there is no possibility of having
the sets grow,
and the sequences are (trivially) $A_{i+j} = A_i$ and $c_{i+j} = c_i$.
By solution, we mean the pair $(u;\lambda)$ which satisfies
\begin{equation}
( u,v )  =  \lambda \langle H(u-c), v \rangle
   \quad \forall v \in H_0^1(\Omega), \quad\mbox{and}\quad
 u = 0 \quad  \mbox{on } \partial \Omega            ;
 \label{eq:leq}
\end{equation}
or; in other words,
the pair $(u;\lambda)$ is a weak solution of the problem.
Consequently, if $c = \lim c_i$,

\begin{theorem} \label{thm2}
The pair $(u/c; \lambda /c)$ is a solution of problem (\ref{eq:basic}).
\end{theorem}

\begin{proof}
The  result of this theorem follows immediately from Lemma~\ref{thm:awmax},
Lemma~\ref{thm:cwmax} which shows that a solution to the free boundary problem
may be rescaled, and Theorem 1 which shows that the free boundary problem
produces a solution in $\Omega$ which has
the value of $k$ on the boundary of the free boundary set $A$.
\end{proof}

Observe, that this rescaling can be done at convergence, or
we may rescale each iterate so that at each stage
in the construction,
\begin{equation}
\min_{x  \in  A_{i-1}} u_{i}(x)  =  1.
\end{equation}
This latter approach has the numerical advantage of keeping the problem
normalized.


\section{Some properties of solutions}

The reason for introducing the pair $(u;\lambda)$ is that
the free parameter in the problem in the solution algorithm
is the estimator of the free boundary set $A_0$
and not $\lambda$.
The value of  $\lambda$ is determined by the value of $c$.
This is a familiar result in the isoperimetric variational approach
where $\lambda$ is determined by the choice of the constraint set.
Even in the elementary example of the radial solution
on the unit disk, it was inconvenient to take $\lambda$ as a
freely defined parameter in the problem.

The isoperimetric variational principle is formulated by considering
the functional
\begin{equation}
J(u) = \int_{\Omega} \int_0^u H(s-1) \: ds \, d \Omega,
\end{equation}
where the domain of integration over $\Omega$ can be restricted to
the set where $H(s-1) \neq 0$, the free boundary set $A$.
Then (\ref{eq:leq}) yields a stationary value $u$ of the restriction of the
bilinear form $(u,u)$ to the set
\begin{equation}
\sigma(\mu)  =  \{u \in H_0^1(\Omega) \mid J(u) = \mu > 0 \}.
 \label{eq:sgm})
\end{equation}
In constructing the isoperimetric variational inequality (in the sense
of \cite{fab}) we constrain the value of $J(u)$ on the free
boundary set.  In choosing $c_0$ using the algorithm in
\eqref{eq:it1}--\eqref{eq:it2}, we constrain
the minimum size of the free boundary set;
that is, we find solutions for which $A_0 \subset A$, where $A_0$
is the set with which we begin the iteration.
Thus, in a sense, these approaches are related in that
in the isoperimetric approach
$J(u)$ is constrained to depend
on $m(A)$ such that $u$ satisfies (\ref{eq:sgm}) while
in the free boundary approach we constrain $m(A) \geq m(A_0)$.

In the convergence proof of Theorem \ref{thm:slw} we used the upper bound
provided by the basic solution to show that the sequence we constructed
converged.  Because the solution of Poisson's equation in any compact domain
is bounded irrespective of the choice of $\chi_{B}$, the sequence of $A_i$
will converge for any initial guess $B \subset \Omega$.
Thus, as an easy consequence of Theorem \ref{thm:slw}, we have
the following result.

\begin{corollary} \label{thm:cslw}
Let $B$ be any subset contained in $\Omega$.
If we solve $-\Delta u_0 = \chi_B$ in place of \eqref{eq:it1},
then the iterative procedure in
Section 2 converges to a solution of \eqref{eq:basic}.
\end{corollary}

Of course,
actually knowing the Green's function for $\Omega$
when $B$ is a ball,
$-\Delta(u)= \chi_B$ with $u=0$ on $\partial(\Omega)$
as discussed in  \cite[Theorem 3.3]{keadykloeden}.
In particular,
the orderly construction of the free boundary set $A$
from any initial iterate gives,


\begin{corollary} \label{thm:infsl}
There are an (uncountably) infinite number of solutions.
\end{corollary}

\begin{proof}
Let $u$ be any solution (\ref{eq:char})
with free boundary set $A$.
On $\Omega$ take as a new starting
value any set $B \supset A$.
Then the sequence of sets $\hat{A}_0, \hat{A}_1, \dots$
converges to a solution $\hat{A}$ of (\ref{eq:basic}).
Since $\hat{A}_0 \supset A$,
the monotonicity of the construction leads to a new solution
with $\hat A \supset A$.
\end{proof}

An additional advantage of the monotonic nesting of the
sets $A_i$ lies in the characterization of the topology of the solutions.
In place of dealing abstractly with a functional
in the variational approach we deal directly with iterates which
converge to the free boundary set.
A simple but important result which applies to the class of solutions
which we construct is the following.

\begin{corollary} \label{coro3}
Every connected component of the free boundary
set obtained by the iterative procedure in Section \ref{sec:sequence}
contains (at least) a local maximum of the initial
solution of (\ref{eq:de}).
\end{corollary}

This corollary is an immediate consequence of
the way we construct the free boundary set,
since $A_i \supset A_0, \forall i$ and $A_0$ has the stated property.
\smallskip

Thus the connectivity of any free boundary solution obtain
using \eqref{eq:it1}--\eqref{eq:it2} is related to the location of the
maxima of the initial iterate $u_0$ and the growth of the sets $A_i$.
Consider using the basic solution as the starting iterate.
By taking a level $c_0$ sufficiently close to the maxima of the
basic solution, the projected set $A_0$ contains components which
are not connected whenever the maxima are distinct, i.e.,
if $\Sigma(c_0)$ contains components which are not connected.
If the sequence of nested sets $A_i$ converges to the free boundary set $A$
before the components of $A_i$
merge, then the the free boundary problem can be seen to admit
solution sets which are not connected.
The connectivity of the free boundary set is seen to be a property
of the initial guess and the growth of the iterates.
Since the maxima of the the initial solution
$-\Delta u_0 = \chi_B$ are easily computable, free boundary
sets which are not connected would appear to be admissible as solutions of
(\ref{eq:basic}) only when it is possible to bound the growth of the iterates.

Finally, consider the rescaling property
of the semi-linear differential equation in Lemma~\ref{thm:cwmax}.
We consider boundary-value problems of the form
$-\Delta u  =  H(u - k/ \lambda)$.
When $\lambda > 0$ and $k \to 0$ or when $\lambda \to \infty$
and $k > 0$, we obtain the basic solution of $-\Delta u = 1$.
For solutions for which the free boundary set is small, then
$\lambda \to \infty$ as $m(A) \to 0$.
This follows since $\lambda \langle \chi_A , v \rangle \to 0$
for fixed $\lambda$ as $m(A) \to 0$,
and the solution to (\ref{eq:basic}) collapses
into the trivial solution unless $\lambda$ is made arbitrarily large.
Based on this, solutions will not be unique for large $\lambda$
since there is the possibility of the following two
solutions sharing the same value of $\lambda$:
one with the free boundary decreasing in measure, approaching a
point in $\Omega$ (necessitating that $\lambda$ become arbitrarily large)
and the other with the free boundary set approaching the domain $\Omega$.
Between these two extremes, the value of $\lambda$ has a minimum, below
which only the trivial solution is possible.  If we let $(u;\lambda)$
be the basic solution to the problem with maximum value $u_{max}$,
then a bound for the minimum
value, $\lambda_{min}$ of (\ref{eq:basic}) is
given by $\lambda_{\rm min} \geq \lambda /u_{\rm max}$.

\subsection*{Conclusion and further work}

Using a technique of successive linear approximations,
we have shown the existence of a class of solutions of a
semi-linear elliptic boundary-value problem with a discontinuity
in the forcing term \eqref{eq:basic}.
A sequence of linear elliptic boundary-value problems is constructed
whose solution converged to a desired solution of the related
free boundary-value problem.
This construction applies in any domain in $\mathbb{R}^2$ or
$\mathbb{R}^3$ in which we can obtain solutions for the
underlying linear elliptic boundary-value problem.
In addition we gain insight into the
connectivity of the free boundary sets for this class of solutions
based on the maxima of the basic
solution.


In practice, preliminary numerical results in applying the
method show that the most significant factor affecting the convergence rate
of this algorithm was the choice of the initial set $B^h$.
Taking $B^h$ based on the basic solution, invariably lead to
solutions which were obtained within a few iterations to within
machine accuracy.
The method was found to be robust
and the convergence rate of the iteration scheme was found to depend
on the shape of the domain and the size of the free boundary set
relative to the mesh being used. In particular, the use of a mesh
fitted to the equipotential lines associated with the
basic solution on the domain $\Omega_h$
is recommended to enhance the
resolution of the method for solutions with small free boundary sets.
Initial results also show that the use of adaptive gridding
of the domain is needed in order to capture the details of
small free boundary sets.

It is worth
observing there may exist solutions to which the numerical procedure
cannot converge since we have not shown that every solution
of the problem belongs to the class computable by this algorithm.
That is, it may be possible for solutions to exist which are not the limits
of monotone sequences we can construct for any choice of initial solution.

\begin{thebibliography}{00}

\bibitem{adams}
R. A. Adams, \emph{Sobolev spaces}, Academic Press, New York, 1975.

\bibitem{aaf}
R. K. Alexander and B. A. Fleishman, \emph{Perturbation and bifurcation in a
  free boundary problem}, Journal of Differential Equations \textbf{45} (1982),
  34--52.

\bibitem{berger}
M. S. Berger, \emph{Non-linearity and functional analysis}, Academic Press, New
  York, 1977.

\bibitem{baf}
M. S. Berger and L. E. Fraenkel, \emph{Non-linear desingularization in certain
  free-boundary problems}, Communications in Mathematical Physics \textbf{77}
  (1980), 149--172.

\bibitem{elcrat}
A. Elcrat, B. Fornberg, M. Horn, and K. Miller, \emph{Some steady vortex flows
  past a circular cylinder}, Journal of Fluid Mechanics \textbf{409} (2000),
  13--27.

\bibitem{fab}
L. E. Fraenkel and M. S. Berger, \emph{A global theory of steady vortex rings
  in an ideal fluid}, Acta Mathematika \textbf{132} (1974), 13--51.

\bibitem{friedman}
A. Friedman, \emph{Variational principles and free-boundary problems}, John
  Wiley \& Sons, New York, 1982.

\bibitem{gilbargandtrudinger}
D. Gilbarg and N. S. Trudinger, \emph{Elliptic partial differential equations
  of second order}, Springer-Verlag, Berlin, 1983.

\bibitem{grisvard}
P. Grisvard, \emph{Behavior of the solutions of an elliptic boundary-value
  problem in a polygonal or polyhedral domain}, Numerical Solution of PDE-III
  (B.~Hubbard, ed.), Academic Press, New York, 1976.

\bibitem{hutsonpym}
V. Hutson and J. S. Pym, \emph{Application of functional analysis and operator
  theory}, Academic Press, New York, 1980.

\bibitem{keady}
G. Keady, \emph{An elliptic boundary-value problem with a discontinuous
  non-linearity}, Proceedings of the Royal Society of Edinburgh \textbf{91A}
  (1981), 161--174.

\bibitem{keadykloeden}
G. Keady and P. E. Kloeden, \emph{An elliptic boundary-value problem with a
  discontinuous nonlinearity, part ii}, Proc.~Roy.~Soc. Edinburgh \textbf{105A}
  (1987), 23--37.

\bibitem{kan}
G. Keady and J. Norbury, \emph{A semilinear elliptic eigenvalue problem, {II}.
  {T}he plasma problem}, Proceedings of the Royal Society of Edinburgh
  \textbf{87A} (1980), 83--109.

\bibitem{littman}
W. Littman, \emph{Generalized subharmonic functions: Monotonic approximations
  and an improved maximum principle}, Ann. Scuola Norm. Sup. Pisa \textbf{17}
  (1963), no.~3, 207--222.

\bibitem{norbury}
J. Norbury, \emph{Steady planar vortex pairs in an ideal fluid}, Communications
  on Pure and Applied Mathematics \textbf{28} (1975), 679--700.

\bibitem{steiner}
G. P\'olya and G.~Szeg\H{o}, \emph{Isoperimetric inequalities in mathematical
  physics}, Princeton University Press, Princeton, 1951.

\bibitem{schecter}
M. Schecter, \emph{Modern methods in partial differential equations},
  McGraw-Hill, 1977.

\bibitem{showalter}
R. E. Showalter, \emph{Hilbert space methods for partial differential
  equations}, Pitman, London, 1977.

\bibitem{turkington}
B. Turkington, \emph{On steady vortex flow in two dimensions, {I} and {II}},
  Communications in Partial Differential Equations \textbf{8} (1983), no. 9,
  999--1071.

\end{thebibliography}



%\bibliography{kolibal}
%\bibliographystyle{amsplain}

\end{document}