\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{gathered} - \Delta \psi = \lambda H(\psi-1) \quad \mbox{in }\Omega, \\ \psi = 0 \quad\mbox{on }\partial \Omega, \end{gathered} \label{eq:basic}$$ 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 $$\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}$$ 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{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}$$ 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{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}$$ The solution to (\ref{eq:line}) is completely determined by imposing $C^1$ continuity of $\psi$ across the free boundary at $r=a$, thus $$\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}$$ 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 $$\Sigma(c) = \{(u,x) : u(x) \geq c, x \in \Omega \} ,$$ 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{gathered} -\Delta u_0 = \lambda \chi_B \quad\mbox{in }\Omega, \\ u_0 = 0 \quad\mbox{on }\partial \Omega, \end{gathered} \label{eq:de}$$ 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{gathered} -\Delta u_0 = \lambda \quad \mbox{in }\Omega,\\ u_0\big|_{\partial \Omega} = 0. \end{gathered}\label{eq:it1}$$ 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 $$A_0 = \{x \in \Omega : u_0(x) \geq c_0 \}. \label{eq:adef}$$ Continuing recursively, for $i \geq 1$ we denote by $u_i$ the solution of $$-\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}$$ and determine $$c_i = \min_{x \in A_{i-1}} u_i(x). \label{eq:cdef}$$ By projecting $\Sigma (c_i)$ onto $\Omega$, the set $$A_i = \{x \in \Omega : u_i(x) \geq c_i \} \label{eq:it2}$$ 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 $$\| u \|_{H^m} = \Big\{ \sum_{| k | \leq m} \| D^k u \|^2_{L_2} \Big\} ^{1/2},$$ 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 $$( u,v ) = \int_{\Omega} \nabla u \cdot \nabla v \, d \Omega,$$ It will also be convenient to introduce the inner product of two functions in $L_2$, $$\langle u,v \rangle = \int_{\Omega} u v \, d\Omega .$$ 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 $$(u,v) = \lambda \langle \chi_A, v \rangle, \qquad \forall v \in H_0^1(\Omega);$$ 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 $$( 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}$$ 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, $$\min_{x \in A_{i-1}} u_{i}(x) = 1.$$ 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 $$J(u) = \int_{\Omega} \int_0^u H(s-1) \: ds \, d \Omega,$$ 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 $$\sigma(\mu) = \{u \in H_0^1(\Omega) \mid J(u) = \mu > 0 \}. \label{eq:sgm})$$ 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}