\magnification = \magstephalf \hsize=14truecm \hoffset=1truecm \parskip=5pt \nopagenumbers \input amssym.def \font\eightrm=cmr8 \font\eighti=cmti8 \font\eightbf=cmbx8 \headline={\ifnum\pageno=1 \hfill\else% {\tenrm\ifodd\pageno\rightheadline \else \leftheadline\fi}\fi} \def\rightheadline{EJDE--1995/11\hfil A Numerical Scheme\hfil\folio} \def\leftheadline{\folio\hfil P.W. Bates, X. Chen, and X. Deng \hfil EJDE--1995/11} \voffset=2\baselineskip \vbox {\eightrm\noindent\baselineskip 9pt % Electronic Journal of Differential Equations Vol. {\eightbf 1995}(1995) No. 11, pp. 1--28.\hfill\break ISSN 1072-6691: URL: http://ejde.math.swt.edu (147.26.103.110)\hfil\break telnet (login: ejde), ftp, and gopher access: ejde.math.swt.edu or ejde.math.unt.edu} \footnote{}{\vbox{\hsize=10cm\eightrm\noindent\baselineskip 9pt % 1991 {\eighti Subject Classification:} 35R35, 65C20, 65M06, 65R20, 82C26. \hfil\break {\eighti Key words and phrases:} Boundary Integral, Coarsening, Free Boundary Problem, Motion by Curvature, Ostwald Ripening. \hfil\break \copyright 1995 Southwest Texas State University and University of North Texas.\hfil\break Submitted June 14, 1995. Published August 18, 1995.\hfil\break Supported by: NSF grant DMS 9305044 (PWB), NSF grant DMS 9404773 and Alfred P. Sloan Research Fellowship (XC).} } \bigskip\bigskip \centerline{A NUMERICAL SCHEME FOR THE TWO} \centerline{PHASE MULLINS-SEKERKA PROBLEM} \smallskip \centerline{Peter W. Bates, Xinfu Chen, and Xinyu Deng} \bigskip\bigskip {\eightrm\baselineskip=10pt \narrower \centerline{\eightbf Abstract} An algorithm is presented to numerically treat a free boundary problem arising in the theory of phase transition. The problem is one in which a collection of simple closed curves (particles) evolves in such a way that the enclosed area remains constant while the total arclength decreases. Material is transported between particles and within particles by diffusion, driven by curvature which expresses the effect of surface tension. The algorithm is based on a reformulation of the problem, using boundary integrals, which is then discretized and cast as a semi-implicit scheme. This scheme is implemented with a variety of configurations of initial curves showing that convexity or even topological type may not be preserved. \bigskip} \input epsf \def\gammat{\Gamma_{\!t}} \def\pupn{\partial u \over \partial n} \def\pwpn{\partial W_g \over \partial n} \def\intgamma{\int_\Gamma} \def\intgammao{\int_{\Gamma_0}} \def\intgammat{\int_{\gammat}} \def\xtoinf{|x|\rightarrow\infty} \def\ov2pi{1\over2\pi} \def\xyeps{x_0-y\pm\epsilon n(x_0)} \def\xyepsf{\xyeps\over{|\xyeps|^2}} \def\xzeroyf{{x_0-y}\over{|x_0-y|^2}} \def\xzeroyepsf{{x_0-y}\over{|\xyeps|^2}} \def\epsxyepsf{\epsilon n(x_0)\over{|\xyeps|^2}} \def\epsxyepsfnon{\epsilon\over{|\xyeps|^2}} \def\xyf{{x-y}\over{|x-y|^2}} \def\tfpmeps{t^2+{\left(f(t)\pm\epsilon\right)}^2} \def\sqrfprim{\sqrt{1+f^{'2}(t)}} \def\xf{x\over|x|^2} \def\oxsqr{O\left(1/|x|^2\right)} \def\ox{O\left(1/|x|\right)} \def\zrzh{z_R(t+h)-z(t+h)} \centerline{\bf 1. Introduction} \medskip Many formulations of time-dependent, multibody free-boundary problems involve the solution of Laplace's equation in an irregular domain at each instant of time as the free boundary evolves. For such problems boundary integral techniques have been used and are a natural choice since they are adaptive in that information only at the free surface is used to predict its motion. Thus the dimension of the problem is reduced by one. Moreover, the discretisation of the free surface can be done to resolve areas of high curvature, a much more difficult task for other methods. In addition, boundary integral techniques automatically account for far field asymptotic behaviour. Boundary integral techniques have been used successfully to study a variety of free surface flows. Examples include the propagation of waves on a fluid-fluid interface [4], the Rayleigh-Taylor instability problem [3, 18], Hele-Shaw flow [9, 15], and crystal growth [14]. McFadden, Voorhees, Boisvert and Meiron examined a multibody free boundary problem as it applies to Ostwald ripening [17, 21, 22]. The model is a quasistatic version of one proposed by Mullins and Sekerka to study the morphological development of a particle [19]. In this case, the free boundary consists of the interface separating two distinct thermodynamic phases, and the dynamics of the process is due to the tendency of the system to minimize the total interfacial surface area by diffusing material from smaller particles to larger ones through the surrounding region of second-phase material. The computational results in [17, 21], based on a boundary integral formulation, give a detailed description of the local behavior of single particles and small collections of particles. However, diffusion within particles is neglected and the numerical analysis employs an explicit time-stepping procedure. Because of the explicit nature of the scheme, to maintain stability the algorithm requires the length of the time-step to be scaled as the cube of the distance between mesh points. Therefore, this procedure is unable, in practice, to compute close to the disappearance of a particle. Implementation was done on a supercomputer with vectorized code. Even then for a four particle system over five and a half hours of CPU time was required to compute until the smallest particle reached a mean radius of $.3$, its initial radius being $.9$. Our purpose is to give a numerical treatment of the two phase Mullins-Sekerka problem in two dimensions employing an algorithm which does not suffer from severe stability constraints on the time step. We present a semi-implicit scheme based on a different boundary integral formulation of the problem. Implementation is done on a workstation, performing experiments with a variety particle configurations, including a four particle system examined in [21]. While the single phase model seems more appropriate when considering solid particles of solute in a saturated solvent, the two phase model may be more realistic when considering diffusion of atomic species in binary solids. The particular model we treat can be described as follows: Given a simple closed curve (or finite collection of nonintersecting simple closed curves), $\Gamma$, in the plane, find a harmonic function in ${\Bbb R}^2\backslash\Gamma$ which on $\Gamma$ is equal to the curvature of $\Gamma$. This function is continuous but not smooth across $\Gamma$ (unless $\Gamma$ is a collection of circles of equal radii). Now evolve $\Gamma$ with normal velocity equal to the jump in the normal derivative across $\Gamma$ of the harmonic function. The problem can also be posed in a bounded domain, $\Omega$, by imposing homogeneous Neumann conditions on $\partial\Omega$. When $\Omega={\Bbb R}^2$ we require a certain decay on the gradient as $\vert x\vert\to\infty$. Pego [20] derived this model, using formal asymptotic expansions, to describe the evolution of level sets of layered solutions to the Cahn-Hilliard equation. In [1] this connection was rigorously proved assuming the existence of smooth solutions to the Mullins-Sekerka flow. At the same time uniqueness of solutions was established. In [5], the existence of weak solutions was proved and later in [6] the existence of smooth solutions was established (see also Constantin and Pugh [8] for results concerning a related problem). In [5] it is shown that the evolution is curve shortening and area preserving so one is naturally drawn to the questions of whether the flow preserves convexity of a single particle, whether or not a single particle can split in two, and whether or not two particles can coalesce. In [16], Mayer proved that for the one (exterior) phase problem, convexity could be lost. In our work, numerical evidence is presented which suggests that this is also the case for the two phase problem. Experiments are also presented which suggest that pinching-off and coalescence are not excluded by this model. In section 2 we present the model in its original form and establish an equivalent formulation using boundary integrals. This is then used to devise a semi-implicit algorithm for advancing the curve. This involves linearizing the operator which maps normal velocity to curvature of the curve at the next time step under the corresponding explicit scheme. In section 3 we discretize the algorithm, computing geometric quantities such as tangent vectors and curvatures using natural differencing ideas. In section 4 we present the results of experiments performed by using the procedure developed, starting with an accuracy check with a configuration for which we know the solution, namely, the evolution of concentric circles. Here we use a variety of time step sizes and numbers of mesh points to aid in our selection of these parameters for subsequent experiments and also to estimate the convergence rate as mesh size and time step tend to zero. However, no attempt is made here to rigorously establish convergence of our general scheme. The next experiment shows the evolution of a nonconvex ``rose'' curve as it quickly becomes convex and tends to a circle, keeping the enclosed area constant. The third experiment starts with a convex particle, a tube with circular endcaps, which at first loses its convexity and then regains it, eventually becoming circular. We then conduct experiments to show that a single particle can ``pinch off'' to become two and two particles can coalesce to become one. Our final experiment uses a configuration of four particles arranged as in an experiment reported in [21] for the one phase problem. Our code, written in $C$, is contained in an appendix to this paper. After this work was completed we were made aware of several interesting results for related problems. In [12] the authors present a way to handle integrals with singular parts arising in fluid interface problems such as Hele-Shaw. A detailed analysis of singularity formation is given in [2] where numerical and analytical methods provide convincing evidence for the formation of singularities in Hele-Shaw flow. Also, the second author, with Hou and Zhu, have incorporated some ideas from [12] with other time saving techniques to substantially accelerate an algorithm similar to that presented here. \noindent{\bf Acknowledgement:} The first author would like to thank Giorgio Fusco for helpful suggestions and discussions. We would also like to thank the referee for a careful reading of the manuscript. Much of this work was reported in the Master's thesis [10] of the third author. \bigskip \centerline{\bf 2. Equivalent formulations and algorithm }\medskip Let $\Omega$ be a bounded and simply connected domain in ${\Bbb R}^2$ and $\Gamma_0$ be a smooth simple closed curve (or finite collection of such) in $\Omega$. Consider the free boundary problem of finding a function $u(x,t), x \in \Omega, t \geq 0$, and a free boundary $\Gamma_{\!0,T} = \cup_{0\leq t\leq T} \left( \gammat \times \{t\}\right) $ for some $T>0$ satisfying $$\cases{\ a)\quad \Delta u(\cdot,t) = 0&in $\Omega \backslash \gammat$, $0 \leq t \leq T$, \cr \ \cr \ b)\quad {\pupn} = 0&on $\partial\Omega \times [0, T]$, \cr \ \cr \ c)\quad u = K&on $\gammat$, $0 \leq t \leq T$, \cr \ \cr \ d)\quad -[{\pupn}]_{\gammat} = V&on $\gammat$, $0 \leq t \leq T$, \cr \ \cr \ e)\quad \Gamma_{\!0,T} \cap \{t = 0\} = \Gamma_0, \cr} \eqno(2.1)' $$ where $n$ is the unit outward normal to $\partial\Omega$ or to $\gammat$, $[{\pupn}]_{\Gamma_{\!\!t}}$ is the sum of the {\it outward} normal derivatives of $u$ from each side of $\gammat$ enclosing a bounded region (which is also equal to the jump of the normal derivative of $u$ across $\gammat$), $K$ is the curvature, taking the sign convention that convex curves have positive curvature, while $V$ is the normal velocity of $\gammat$ and the normal velocity of an expanding curve enclosing the bounded region is positive. For simplicity, we only consider the case in which $\Omega = {\Bbb R}^2$. In this case, we replace the boundary condition (2.1b)$'$ by $$ \nabla u(\cdot, t) = \oxsqr \qquad\qquad \hbox{as }\xtoinf$$ and call the resulting system (2.1). We first derive an integral representation for the solution to Eqs. (2.1). \medskip \noindent {\bf Lemma 2.1.} {\sl Let $\Gamma = \cup_{i=1}^M\Gamma^i$ $(M\geq 1)$, where $\Gamma^i$ are disjoint simple closed curves such that $\Gamma$ separates ${\Bbb R}^2$ into one unbounded and $m$ bounded regions. Let $n$ be the outward unit normal to $\Gamma$. For each $g \in L^2(\Gamma)$, define $$ W_g(x) = {\ov2pi} \int_\Gamma \ln |x-y|\>g(y)\>ds_y,\qquad x \in {\Bbb R}^2\backslash \Gamma.$$ Then the following holds: (1) $\Delta W_g = 0 \qquad$ in ${\Bbb R}^2 \backslash \Gamma$; (2) $\displaystyle{-\left[{\pwpn}\right]}_\Gamma \equiv {\pwpn}^{out} - {\pwpn}^{in} = g$ \qquad on $\Gamma$; (3) if in addition we assume that$\>\displaystyle{\intgamma g(y)\>ds_y\>\>=\>\>0}$, \quad then \qquad$\displaystyle{|{\nabla W_g}| = \oxsqr}$ \quad and \quad $\displaystyle{W_g = \ox}$ \qquad as $\xtoinf$; (4) The mapping $g\in\{ g \in L^2(\Gamma)|\intgamma g = 0\}\longmapsto W_g$ is negative definite, i.e. \qquad $\displaystyle{(g,W_g) \equiv \intgamma g\>W_g}\> < 0$\qquad if $g\not\equiv0$.} \noindent {\bf Proof:} The first assertion of the lemma follows from the fact that $$ \Delta_x \ln |x-y|=0 \qquad \hbox{if } x\not=y.$$ The second is a standard calculation in potential theory (see for instance [11] section 3D or [13]). To see (3), note that when $|x|$ is large, $$\eqalign{\nabla W_g &= {\ov2pi}\intgamma {\xyf}\>g(y)\>ds_y\cr &= {\ov2pi}\intgamma {\left({\xyf}-{\xf}\right)}\>g(y)\>ds_y,\cr}$$ since $\intgamma g(y)\>ds_y=0$. Also, $$\eqalign{{\xyf}-{\xf} &={{|x|^2-|x-y|^2}\over{|x-y|^2\cdot|x|^2}}\>x-{y\over{|x-y|^2}}\cr &= \oxsqr.\cr}$$ Hence, $$\nabla W_g(x)=\oxsqr.$$ Similarly, $$\eqalign{W_g &= {\ov2pi}\>\intgamma \left[{\ln|x-y|\>g(y)-\ln|x|\>g(y)}\right]\>ds_y\cr &= {\ov2pi}\>\intgamma \ln\left({|x-y|\over|x|}\right)\>g(y)\>ds_y\cr &= {\ov2pi}\>\intgamma \ln\left({1+\ox}\right)\>g(y)\>ds_y\cr &= \ox.\cr}$$ The third assertion of the lemma then follows. Finally, set $u=W_g$, then $$\eqalign{-\intgamma W_g(y)g(y)\>ds_y &= \intgamma u \left[{\pupn}\right]_{\Gamma} \>ds_y\cr &= \int\!\!\!\int_{{\Bbb R}^2} {u\Delta u} +\int\!\!\!\int_{{\Bbb R}^2} {\nabla u\cdot \nabla u}\cr &= \int\!\!\!\int_{{\Bbb R}^2} {|\nabla u|^2}\>dxdy > 0,\cr}$$ where in the second equation we have used the assertion of (3). This completes the proof of the lemma.\quad\vrule height4pt width3pt The previous lemma specifies a zero mean value jump across $\Gamma$ of the normal derivative of a proposed function which is harmonic on each side of $\Gamma$ and produces such a harmonic function. Clearly, any constant can be added to the harmonic function produced by the lemma giving another suitable function. The next lemma shows that any harmonic function on ${\Bbb R}^2\backslash\Gamma$ having a sufficiently regular trace on $\Gamma$ has this trace realized by the harmonic function constructed above from the jump in the normal derivative, up to an additive constant. \medskip \noindent {\bf Lemma 2.2.} {\sl Suppose that $\Gamma$ is as in Lemma 2.1. Suppose that $u$ defined on ${\Bbb R}^2$, $f\in H^1(\Gamma)$ and $g\in L^2(\Gamma)$ are related by $$ \left.{\matrix{\Delta u &=0\hfill \hbox{ \quad \ \ in } {\Bbb R}^2\backslash\Gamma\hfill\cr \vert\nabla u\vert &=O\left({{1\over \vert x\vert^2}}\right)\hfill \hbox{ as }\vert x\vert\to\infty\hfill\cr u &=f\hfill \hbox{ \ on }\Gamma\hfill\cr -\left[{{\partial u\over\partial n}}\right]_{\Gamma}&= g\hfill \hbox{ \ on }\Gamma.\hfill}}\right\}\eqno(2.2)$$ Then there exists a constant $c$ such that $$f(x)={1\over 2\pi}\int_{\Gamma} \ln\vert x-y\vert g(y)ds_y+c\hbox{ for } x\in\Gamma.\eqno(2.3)$$ Furthermore,} $$\int_{\Gamma} g(y)ds_y=0.\eqno(2.4)$$ \noindent {\bf Proof:} Define $W_g$ as in Lemma 2.1 and $$c(x)=u-W_g.$$ Then $\Delta c=0$ in ${\Bbb R}^2\backslash\Gamma$ and $c$ is continuous across $\Gamma$ since both $u$ and $W_g$ are. Furthermore, $\left[{{\partial c\over \partial n}}\right]_{\Gamma}=0$ and hence $\Delta c=0$ in ${\Bbb R}^2$. Finally, $\nabla c= \oxsqr$ as $\vert x\vert\to\infty$ and so $c$ is a constant. This yields (2.3). To establish (2.4), let $B_R$ be a disk of radius $R$, then, as $R\to\infty$, $$\int_{\Gamma} g=-\int_{\Gamma} \left[{{\partial u\over\partial n}}\right] =\int_{\partial B_R} {\partial u\over\partial n}-\int\!\!\!\!\int_{B_R\backslash\Gamma} \Delta u=2\pi R\cdot O(1/R^2)\to 0.\eqno \quad\vrule height4pt width3pt$$ The above lemmas can now be combined to show the equivalence of (2.1) and the system (2.3)--(2.4) with $f=K$ and $g=V$, the curvature and normal velocity, respectively, of $\Gamma$. \medskip \noindent {\bf Theorem 2.3} {\sl If $\Gamma_{0,T}\equiv\cup_{0\leq t\leq T} \left({\Gamma_t\times\{t\}}\right)$ is a continuous family of $C^3$ curves which satisfy (2.3) and (2.4) with $\Gamma=\Gamma_t$ for some $g=G(x,t)$ and $c=c(t)$, where $f=K(x,t)$ is the curvature of $\Gamma_t$ at $x$, then $\Gamma_{0,T}$ is the interface associated with the solution to (2.1) and $g$ is the normal velocity, $V$, of $\Gamma_t$. Conversely, if $(u,\Gamma_{0,T})$ is a solution to (2.1) then (2.3)--(2.4) hold for each $t\in [0,T]$ with $\Gamma=\Gamma_t, g=V$ and $f=K$.} \noindent {\bf Proof:} Let $(\Gamma_t, g, c)$ be a solution to (2.3)--(2.4) at each time $t\in [0,T]$ with $f=K(\cdot, t)$. Then Lemma 2.1 shows that defining $u(\cdot, t)$ on ${\Bbb R}^2\backslash\Gamma_t$ by $$u(x,t)={1\over 2\pi} \int_{\Gamma_t} \ln\vert x-y\vert g(y,t)ds_y+c(t)$$ gives a solution to (2.1) with the normal velocity of $\Gamma_t$ being given by $g$. Solutions to (2.1) with initial curve $\Gamma_0$ are unique by the results in $[1]$ and so $\Gamma_{0,T}$ is that solution. Conversely, if $(\Gamma_{0,T}, u)$ is a solution to (2.1) then the solution is smooth by the results in [6]. By Lemma 2.2, (2.3)--(2.4) hold for some $c$ with $\Gamma=\Gamma_t, f=K(\cdot, t)$, and $g=V(\cdot, t)$, for each $t\in [0,T]$.\quad\vrule height4pt width3pt Before we present a numerical scheme for solving (2.3)--(2.4), we first solve an inverse problem: \medskip \noindent {\bf Lemma 2.4} {\sl Given $f\in H^1(\Gamma)$, with $\Gamma$ as in Lemma 2.1, there exists a unique $g\in L^2(\Gamma)$ and constant $c$ such that (2.3) and (2.4) hold.} \noindent {\bf Proof:} We first prove uniqueness. Since (2.3), (2.4) is a linear system it suffices to show that $f\equiv 0$ implies $g\equiv 0$ and $c=0$. Let $w=W_g+c$, then $$\Delta w=u\, \hbox{ in }\, {\Bbb R}^2\backslash\Gamma$$ and on $\Gamma, w=f=0$. Furthermore, by Lemma 2.1, $\vert\nabla w\vert=\oxsqr$ as $\vert x\vert\to\infty$. It follows that $w\equiv 0$ on ${\Bbb R}^2$ and by Lemma 2.1 (2), $g\equiv 0$ on $\Gamma$. Consequently, $c=0$ also. To prove existence of a solution let $u$ satisfy $$\eqalign{\Delta u&=0 \hbox{ on } {\Bbb R}^2\backslash\Gamma,\cr u&=f \hbox{ on } \Gamma,\cr \vert\nabla u\vert&= \oxsqr \hbox{ as } \vert x\vert\to\infty.}$$ Since $f\in H^1, u$ exists, is unique, and $u\in H^2$. Set $$g= -\left[{{\partial u\over\partial n}}\right]_{\Gamma}.$$ Then, $g\in L^2(\Gamma)$ and by Lemma 2.2, $\int_{\Gamma}g=0$ and $c\equiv u-W_g$ is constant.\quad\vrule height4pt width3pt From the previous discussion, (2.1) can now be solved step by step as follows: Assume that $\Gamma_t$ is known. Calculate $f=K(\cdot,t)$ on $\Gamma_t$, then find $g=V(\cdot, t)$ and $c(t)$ by solving $$\eqalign{&K(x,t) = {\ov2pi}\>\int_{\Gamma_t} \ln|x-y|V(y,t)\>ds_y + c(t),\cr &\int_{\Gamma_t} V(y,t) \>ds_y = 0.\cr}\eqno(2.5)$$ Once we know $V$, we can advance the curve by $$x(t+dt)=x(t)+Vndt.\eqno(2.6)$$ This is an explicit scheme. The weakness of the scheme is that the time step has to be extremely small due to instablity. A small time step, however, introduces more machine error. We avoid these problems through the following semi-implicit scheme. It linearizes the mapping which uses the explicit method to get the curvature at one time step from the normal velocity at the previous step: Write $$\eqalign{&\Gamma_t =\left\{ x=x(s,t)\>\>\vert \>\>s\in S^1=R^1 \>\> (\hbox{mod} \>\> 2\pi)\right\},\cr &k(s,t) = K(x(s,t), t),\cr &V(s,t) = {\partial x\over\partial t}\cdot n.\cr}$$ Assume $x(s,t)$ is known for some t, we find $V$ by solving $$\left\{\eqalign{&k(s,t)+BV\Delta t = {\ov2pi}\>\int_{\Gamma_t} \ln|x(s,t)-y|V\>ds_y+c,\cr &\int_{\Gamma_t} V(y,t) \>ds_y = 0,\cr}\right.\eqno(2.7)$$ where $$BV={\partial K(x(s,t)+hnV)\over\partial h}\Big\vert_{h=0}$$ defines a linear operator. We then advance $\Gamma$ according to (2.7). Note that $k(s,t)+BV\Delta t$ is an approximation of the curvature at $t+\Delta t$. \bigskip \centerline{\bf 3. Discretization }\medskip Assume that $\Gamma_t=\cup_{m=1}^M{\Gamma^m_t}$ consists of $M$ disjoint simple closed curves and take $N$ points from each $\Gamma^m_t$, labelling them by $z_{(m-1)N+1}$, $z_{(m-1)N+2}$, $\dots$, $z_{mN}$ listed counterclockwise. Then, to deal with the periodicity, we define the following permutation functions: $$\eqalign{L[i]&=\cases{\ i-1&if ${{i-1}\over N} \not= $ integer, \cr \ i-1+N&if ${{i-1}\over N} = $ integer, \cr}\cr R[i]&=\cases{\ i+1&if ${i\over N} \not= $ integer, \cr \ i+1-N&if ${i\over N} = $ integer, \cr}\cr} \eqno(3.1)$$ to denote the indices of the points clockwise and counterclockwise from $z_i$, respectively. Assume that $z_L$, $z$, $z_R$, are three counterclockwise consecutive points on $\Gamma$, we interpolate $\Gamma$ near $z$ as a segment of the circle passing through $z_L$, $z$, $z_R$. We use the unit tangent, $\tau$, the outward unit normal, $n$, and curvature $K$ of the circle passing through $z_L$, $z$, $z_R$ as the corresponding quantities of $\Gamma$(see Fig. 1). If $z_L$, $z$, $z_R$ are colinear, then the curvature is zero and the tangent is the obvious one. We denote $$\eqalign{T_L&={z-z_L\over|z-z_L|}, \qquad T_R={z_R-z\over |z_R-z|},\cr d_L&=|z-z_L|, \qquad d_R=|z_R-z|,\cr \hbox{ and }\qquad\qquad d_{RL}&= |z_R-z_L|.}\eqno(3.2)$$ Note that if a tangent vector has coordinates $\tau=(\tau^x,\tau^y)$ then $n=(\tau^y,-\tau^x)$ is a unit normal. This motivates the definitions $$N_L=(T^y_L,-T^x_L)\hbox{ and } N_R=(T^y_R,-T^x_R)$$ as approximate normals, which are outward by our ordering of the $z_j$'s. The following is easily demonstrated and we omit the proof: \centerline{\epsfbox{Fig1.ps}} \centerline{Fig. 1 Interpolation} \medskip \noindent{\bf Lemma 3.1.} {\sl Let $C$ be the circle passing through noncolinear $z_L,z,z_R$ and let $\tau,n$ and $K$ be the counterclockwise unit tangent, outward unit normal, and curvature of $C$ at $z$, respectively. Then,} $$\eqalign{\tau&={d_RT_L+d_LT_R\over d_{RL}},\cr n&={d_RN_L+d_LN_R\over d_{RL}},\cr K&={2T_L\cdot N_R\over d_{RL}}=-{2T_R\cdot N_L\over d_{RL}}.}\eqno(3.3)$$ Let $\Gamma=\cup_{j=1}^{MN}{\Gamma_j}$, where $\Gamma_j$ is a segment of $\Gamma $ containing $z_j$ in its interior. We take these segments small so that a given function $g$ on $\Gamma$ is almost constant on $\Gamma_j$. Define $$W_g(z') = {\ov2pi}\int_{\Gamma} \ln|z'-z|g(z)\>ds_z\eqno(3.4)$$ and set $d_i = \int_{\Gamma_i} ds_i$, $W^i={1\over d_i}\int_{\Gamma_i} W_g(z')\>ds_z$, and $g_j=g(z_j)$. Then $$W^i \simeq \sum_{j=1}^{MN} a_{ij}g_j \quad \hbox{for}\quad i= 1,2,\dots,MN, \eqno(3.5)$$ where $$a_{ij} = {1\over{2\pi d_i}} \int_{\Gamma_{\!i}}\!\int_{\Gamma_{\!j}} \ln|z'-z| \>ds_{z'}ds_z.\eqno(3.6)$$ The remaining work is for the approximation of $a_{ij}$. \vbox{ \epsfbox{Fig2.ps} \centerline{Fig. 2. Approximation of $a_{ij}$}} We are seeking second order approximation, so we can replace $\Gamma_i$ and $ \Gamma_j$ in (3.6) by line segments $\hat \Gamma_i$ and $\hat \Gamma_j$ (Fig. 2). Define $\tau_i$ as the unit tangent at point $z_i$, as in (3.3), and take $$\eqalign{d &= d_{ij} = |z_i-z_j|,\cr d_i &= {1\over 2} \left(d_{R[i]}+d_{L[i]}\right),\cr \tau_{ji} &= \cases{\ \tau_i&if $i=j,$ \cr \ {{z_j-z_i}\over d}&if $i\not= j,$ \cr}\cr A &= \cos\alpha = -\tau_i\cdot \tau_{ji},\cr B &= \cos\beta = \tau_j\cdot \tau_{ji}.\cr}$$ We have $$ \eqalign{ a_{ij} &= {1\over {2\pi} d_i} \int_{\Gamma_{\!i}}\!\int_{\Gamma_{\! j}} \ln|z'-z|\>ds_{z'}ds_z\cr &\cong {1\over {2\pi} d_i}\int_{-{d_{L[i]}\over 2}}^{d_{R[i]}\over 2}\!\int_{-{d_{L[j]}\over 2}}^{d_{R[j]}\over 2} \ln|d+Ax+By|\>dy dx \cr &= \Biggl[{1\over{4{\pi}d_iAB}}[d+Ax+By]^2\left(\ln|d+Ax+By|-{3\over 2}\right) \Biggr]\Biggl|_{y=-{d_{L[j]}\over 2}}^{y={d_{R[j]}\over 2}} \Biggr|_{x=-{d_{L[i]}\over 2 }}^{x={d_{R[i]}\over 2}} \cr &={1\over {2{\pi}AB\left(d_{R[i]}+d_{L[i]}\right)}} \Biggl\{a_1^2\left(\ln|a_1|-{3\over 2}\right)+a_2^2\left(\ln|a_2|-{3\over 2}\right)}$$ $$-a_3^2\left(\ln|a_3|-{3\over 2}\right)-a_4^2\left(\ln|a_4|-{3\over 2}\right) \Biggr\} \equiv \hat a_{ij},\eqno(3.7)$$ where $$ \eqalign{a_1 &= d+A{d_{R[i]}\over 2}+B{d_{R[j]}\over 2}, \cr a_2 &=d-A{d_{L[i]}\over 2}-B{d_{L[j]}\over 2}, \cr a_3 &= d+A{d_{R[i]}\over2}-B{d_{L[j]}\over 2}, \cr a_4 &= d-A{d_{L[i]}\over2}+B{d_{R[j]}\over 2}. \cr}$$ A more detailed analysis shows that $|a_{ij}-\hat a_{ij}| = O(r^2)$, where $r=\max |z_i-z_{R[i]}|$. Note that $W^i\simeq W_g(z_i)$ and so with $V_j=V(z_j)$ and $K_j=K(z_j)$ for fixed $t$, system (2.5) is discretized to give the following system, where we denote $\hat a_{ij}$ by $a_{ij}$. $$ \eqalign{&\sum_{j=1}^{MN}a_{ij}V_j + C = K_i\quad \hbox{for} \quad i=1,\dots, MN,\cr &\sum_{j=1}^{MN}d_jV_j = 0,\cr}\eqno(3.8)$$ where $$ \eqalign{U &\equiv (V_1, \dots, V_{MN}, C)^T \hbox{ are unknowns,}\cr K_i &= {-2T_{R[i]}\cdot N_{L[i]}\over d_{RL[i]}},\cr d_j &= {1\over 2} \left(d_{R[j]}+d_{L[j]}\right),\cr a_{ij} &= {1\over {2{\pi}AB\left(d_{R[i]}+d_{L[i]}\right)}} \Biggl\{a_1^2\left( \ln|a_1|-{3\over 2}\right)+a_2^2\left(\ln|a_2|-{3\over 2}\right) \cr &\qquad-a_3^2\left(\ln|a_3|-{3\over 2}\right)-a_4^2\left(\ln|a_4|-{3\over 2}\right)\Biggr\}, \cr}$$ where $a_1$, ..., $a_4$ are as above with $d=d_{ij}=|z_i-z_j|$. Hence, we find the solutions $U$, by solving (3.8) and the explicit scheme updates the location of the interface through the formula $$ z(t+h) = z(t)+hn(z(t))V,\eqno(3.9)$$ where $h$ is the time step and $n$ is the outward unit normal. In (3.8), if we let $K_i$ be the curvature evaluated at $z_i(t)$ then the scheme is explicit and unstable unless $h$ is sufficiently small. If we evaluate $K_i$ at $z_i(t+h)$ given in (3.9), then the scheme is implicit and stable for larger $h$. In this case, $K_i$ depends on $V$ in a non-linear and non-local manner. For ease of computation, we take a semi-implicit scheme by extracting the linear part of this dependence and ignoring the remainder. From experiments we performed, we see that this is enough for the stability of the scheme, as well as the accuracy (see section 4). To extract the linear part of the dependence of curvature on velocity, we compute the first derivative of $K$ with respect to $h$. Since $$K=-{{2T_R\cdot N_L}\over d_{RL}}$$ then $$ {\partial K \over \partial h} = -{{2{{\partial T_R}\over{\partial h}}\cdot N_ L}\over d_{RL}}-{{2T_R\cdot{{\partial N_L}\over{\partial h}}}\over d_{RL}}+{{2T_ R\cdot N_L {\partial d_{RL}\over\partial h}}\over d_{RL}^2}.$$ Since $$\eqalign{T_R\cdot {\partial N_L\over\partial h} &=(T_R^x, T_R^y)\left ( {\partial N_L^x\over\partial h}, {\partial N_L^y\over\partial h}\right )\cr &=(T_R^x, T_R^y) \left ( {\partial T_L^y\over\partial h}, -{\partial T_L^x\over\partial h}\right )\cr &=T_R^x{\partial T_L^y\over\partial h} - T_R^y{\partial T_L^x\over\partial h}\cr &=-\left(N_R^y{\partial T_L^y\over\partial h} + N_R^x{\partial T_L^x\over\partial h}\right)\cr &=-N_R\cdot {\partial T_L\over\partial h}\cr}$$ we get $$ {\partial K \over \partial h} = -{{2{{\partial T_R}\over{\partial h}}\cdot N_ L}\over d_{RL}}+{{2N_R\cdot{{\partial T_L}\over{\partial h}}}\over d_{RL}}-{{K{ \partial d_{RL}\over\partial h}}\over d_{RL}}.\eqno(3.10)$$ From (3.2) and (3.9), We have $$ \eqalign{{\partial T_R\over\partial h} &= {\partial \over\partial h}\left(\zrzh\over {|\zrzh|}\right)\cr &= {{{\partial \over \partial h}\left(\zrzh\right)}\over {|\zrzh|}}\cr &-{\zrzh\over{|\zrzh|}^3} \left[(\zrzh) \cdot {\partial \over \partial h} \left(\zrzh\right)\right]\cr &= {{n_RV_R-nV}\over{|\zrzh|}}\cr &-{\zrzh\over{{|\zrzh|}^3}} \left[(\zrzh)\cdot (n_RV_R-nV)\right],\cr}$$ where we use $n$, $n_L$ and $n_R$ to denote unit normals at $z$, $z_L$ and $z_R$, respectively (see Fig. 1). Sending h $\rightarrow$ 0, since $T_R = {z_R-z\over |z_R-z|}$ and $N_R(N_R\cdot x)+T_R(T_R\cdot x) = x,$ we obtain $$ \eqalign{{\partial T_R\over\partial h}\bigg|_{h=0} &= {1\over {|z_R-z|}} \left \{(n_RV_R-nV)-T_R(T_R\cdot (n_RV_R-nV))\right\}\cr &= {N_R\over d_R} \left\{N_R\cdot (n_RV_R-nV)\right\}.\cr}\eqno(3.11)$$ Similarly, $${\partial T_L\over\partial h}\bigg|_{h=0} = {N_L\over d_L} \left\{N_L\cdot (nV- n_LV_L)\right\}.$$ We also compute $$\eqalign{&{\partial d_{RL}\over\partial h} = {1\over d_{RL}} \left\{(z_R(t+h)- z_L(t+h))\cdot (n_RV_R-n_LV_L)\right\}\cr &{\partial d_{RL}\over\partial h}\bigg|_{h=0} = {1\over d_{RL}} \left\{(z_R-z_L) \cdot (n_RV_R-n_LV_L)\right\}.\cr}$$ Then (3.10) can be written as $$ \eqalign{{\partial K \over \partial h}\bigg|_{h=0} &= -{2\over {d_R d_{RL}}}(( n_RV_R-nV) \cdot N_R) (N_R \cdot N_L)\cr &\quad +{2\over {d_L d_{RL}}}((nV-n_LV_L) \cdot N_L) (N_L \cdot N_R)\cr &\quad-{K\over d_{RL}^2}(z_R-z_L)\cdot (n_RV_R-n_LV_L)\cr &=V_L\left(-{2(N_L \cdot n_L) (N_R \cdot N_L) \over{d_L d_{RL}}}+K{(z_R-z_L) \cdot n_L\over d_{RL}^2}\right)\cr &\quad+V\left({2(N_R \cdot n) (N_R \cdot N_L) \over{d_R d_{RL}}}+{2(N_L \cdot n) (N_R \cdot N_L) \over{d_L d_{RL}}}\right)\cr &\quad+V_R\left(-{2(N_R \cdot n_R) (N_R \cdot N_L) \over{d_R d_{RL}}}+K{(z_L-z_R) \cdot n_R\over d_{RL}^2}\right).\cr}\eqno(3.12)$$ Hence, we have $$K(t+h)\simeq K(t)+hBV,\eqno(3.13)$$ where $$h B = \{b_{ij}\}\quad \hbox{for} \quad i,j=1,\dots, MN $$ and $$ \eqalign{&b_{iL[i]} = h\left(-{2(N_L \cdot n_L) (N_R \cdot N_L) \over{d_L d_{RL}}}+K{(d_RT_R+d_LT_L) \cdot n_L\over d_{RL}^2}\right)\cr &b_{ii} = h\left({2(N_R \cdot n) (N_R \cdot N_L) \over{d_R d_{RL}}}+{2(N_L \cdot n) (N_R \cdot N_L) \over{d_R d_{RL}}}\right)\cr &b_{iR[i]} = -h\left({2(N_R \cdot n_R) (N_R \cdot N_L) \over{d_R d_{RL}}}+K{( d_RT_R+d_LT_L) \cdot n_R\over d_{RL}^2}\right)\cr &\qquad\qquad \hbox{for } i=1,\dots, MN,\cr &b_{ij}= 0 \qquad\hbox{for } j\not= i, \> R[i], \hbox{or } L[i]; \>\> i,j=1, \dots, MN.\cr}\eqno(3.14)$$ Now (3.8) can be modified as $$ \left\{\eqalign{&\sum_{j=1}^{MN}(a_{ij}-\beta b_{ij})V_j + C = K_i\quad \hbox{for}\quad i=1,\dots, MN,\cr &\sum_{j=1}^{MN}d_jV_j = 0\cr}\right.\eqno(3.15)$$ where $\beta\in [0,1]$ is a factor indicating the extent of the implicit nature of the scheme. Several experiments using different values for $\beta$ (not reported here) convinced us that taking $\beta=1$ provides a stable scheme while maintaining accuracy, therefore all our subsequent experiments take $\beta=1$. We may write (3.15) in matrix form $$GU = P,\eqno(3.16)$$ where $$G=\pmatrix{a_{ij}-b_{ij}&1\cr d&0\cr}, \quad P=\pmatrix{K\cr 0\cr}, and \quad U=\pmatrix{V\cr C\cr}.$$ After solving (3.16) for $U$ we update $z$ using (3.9). Our algorithm is based on the premise that each boundary component $\Gamma_t^m$ is described by $N$ mesh points. During initialization, the $N$ points are chosen to be equispaced in arc length of the boundary. After each iteration, new mesh points will be generated, representing the evolution of the boundary. These new points may not be equispaced and, unless we are careful, may concentrate at certain locations, leading to computational errors. To avoid this problem, we shall redistribute these newly generated points so that they are almost equispaced. \centerline{\epsfbox{Fig3.ps}} \centerline{Fig. 3 Reparameterization} Let $z_0$ be the center of the circle passing through three newly created points: $z$, $z_L$, and $z_R$, and let $z^*$ be the mid-point of the arc $\widehat {z_Lz_R}$. We shall take $z^*$ as the new location for the point $z$. The outward unit normals $n$, $N_L$, and $N_R$ are given by the definition following (3.2), while $n_{LR}$ represents the outward unit normal at $z^*$ (refer to Fig. 3). To solve for $z^*$, we write $$ \eqalignno{z_0 &= z-{1\over K}n &(3.17)\cr z^* &= z_0 +{1\over K}n_{LR} = z+{1\over K}\left(n_{LR}-n\right). &(3.18)\cr}$$ Replacing $n$ and $K$ in (3.18) using (3.3), we have $$\eqalign{z^*&=z+{d_{RL}\over 2T_L\cdot N_R} \left(n_{LR}-{d_RN_L+d_LN_R\over d_{RL}}\right)\cr &=z+{d_{RL}n_{LR}-d_RN_L-d_LN_R\over 2T_L\cdot N_R}.}\eqno(3.19)$$ Since $$d_{RL}n_{LR} = {\left(d_{RL}T_{RL}\right)}^\perp={\left(d_RT_R+d_LT_L\right)}^ \perp= d_RN_R+d_LN_L,$$ (3.19) can be written as $$\eqalign{z^* &= z+{{d_RN_R+d_LN_L-d_RN_L-d_LN_R}\over{2T_L\cdot N_R}}\cr &= z+{(d_R-d_L)(N_R-N_L)\over{2T_L\cdot N_R}}.\cr}\eqno(3.20)$$ If we write $$N_R-N_L = aT_L +bT_R,$$ then we find $$ \eqalign{a=b &={{1-N_L\cdot N_R}\over {T_L\cdot N_R}}={{1-T_L\cdot T_R}\over {T_L\cdot N_R}}={{1-{(T_L\cdot T_R)}^2}\over {\left(1+(T_L\cdot T_R)\right)(T_L \cdot N_R)}}\cr &= {(T_L\cdot N_R)^2\over {\left(1+(T_L\cdot T_R)\right)(T_L\cdot N_R)}}= {{T_L \cdot N_R}\over {1+(T_L\cdot T_R)}}\cr}\eqno(3.21)$$ and so $$N_R-N_L = {{(T_L\cdot N_R)(T_L+T_R)}\over{1+(T_L\cdot T_R)}}.\eqno(3.22)$$ Substituting into (3.20) gives $$ z^* = z+{{(d_R-d_L)(T_L+T_R)}\over{2\left(1+T_L\cdot T_R\right)}}.\eqno(3.23 )$$ Therefore, after obtaining $z_i=z(t+h)$ , we rearrange $z_i$ as follows before calculating the positions predicted by the next time step. Calculate $$\eqalign{d_{R[i]}&=\vert z_{R[i]}-z_i\vert, \quad d_{L[i]}=\vert z_i-z_{L[i]}\vert,\cr T_{R[i]}&={z_{R[i]}-z_i\over d_{R[i]}}, \quad T_{L[i]}={z_i-z_{L[i]}\over d_{L[i]}},\cr \hbox{then set } \hfil z_i: &= z_i+{{(d_{R[i]}-d_{L[i]})(T_{L[i]}+T_{R[i]})}\over{2\left(1+T_{L[i]}\cdot T_{R[i]}\right)}}.\cr}\eqno(3.24)$$ So our numerical scheme has two steps: Given a curve discretization $\Gamma_t$, compute a discretization of $\Gamma_{t+h}$ by using the semi-implicit algorithm to find $V$, then redistribute the mesh points on $\Gamma_{t+h}$ and repeat the whole process. The results of implementing this in several examples are given below. \bigskip \centerline{\bf 4. Experiments }\medskip First, to test the accuracy of our numerical method, we choose a case which we can solve analytically, namely, the case for two concentric circles. Take concentric circles of radii $R_1dR_1$$ and so by using (4.5) we have $$ t = \int_{R_1}^{R_1^0} {{R_1\ln{R_2\over R1}}\over {{1\over R_1}+{1\over R_2} }}\>dR_1 = {1\over 2}\int_{R_1}^{R_1^0} {{ R_1^2 \ln{{A^2+R_1^2}\over R_1^2}\sqrt{R_1^2+A^2}}\over {R_1+\sqrt{R_1^2+A^2}}}\>dR_1.\eqno(4.6)$$ Let $R_1 = Ar$ then (4.6) can be written as $$t={A^3\over 2}\int_{R_1\over A}^{R_1^0\over A} {{r^2\sqrt{1+r^2} \ln\left(1+{1 \over r^2}\right)}\over{r+\sqrt{1+r^2}}}\> dr.\eqno(4.7)$$ Or we may write $r={1\over k}$ in (4.7), to obtain the curvature-time relationship $$t= {A^3\over 2}\int_{Ak^0}^{Ak} \ln\left(1+k^2\right){\sqrt{1+k^2}\over{\left( 1+\sqrt{1+k^2}\right)k^4}}\> dk,\eqno(4.8)$$ where $$ k^0 = {1\over R_1^0}.$$ Given curvature values $k_1$ and $k_2$ obtained from the numerical simulation of (2.1) over time step $\Delta t$, setting $a=Ak_1$ and $b=Ak_2$, we will compare the results with those obtained by integrating (4.8). To calculate the integral in (4.8), we use Simpson's rule. Denote $$ f(k) = {\sqrt{1+k^2}\over {k^4\left(1+\sqrt{1+k^2}\right)}} \ln\left(1+k^2\right).\eqno(4.9)$$ The discretization of the integral will use step length $h < {\Delta t\over 10}$ by setting $n=[10(b-a)/\Delta t]+1$ and $h=(b-a)/n$. Then the integral over the interval [$a$, $b$] is $$\int_a^b f(k)dk \simeq \sum_{i=0}^{n-1} {h\over 6}\left(f(a+ih)+4f\left(a+ih+{h\over 2}\right)+f(a+ih+h)\right).\eqno(4.10)$$ The above formula has accuracy given by $$\eqalign{\int_\alpha^{\alpha+h} f(k)dk = &{h\over 6}\left(f(\alpha)+4f\left(\alpha+{h\over 2}\right)+f(\alpha+h)\right)\cr &+{1\over 2880} f^{(4)}(\tilde\theta)h^5 \hbox{ for some }\ \ \tilde\theta\in(\alpha,\alpha+h).\cr}$$ Clearly, $f^{(4)}(k)$ is bounded for $k\geq 1$. Hence, $t$ can be computed to an accuracy of order $O(h^4)$ by $$t = {A^3\over 2}\Biggl[\sum_{i=0}^{n-1} {h\over 6}\left(f(a+ih)+4f\left(a+ih+{h\over 2}\right)+f(a+ih+h)\right)\Biggl].\eqno(4.11)$$. We apply our algorithm, (3.16), taking concentric circles of radii $1$ and $3$ and terminate when the radius of the small circle is approximately $0.1$. We display the curvature-time relationship of the ``exact'' solution (obtained using (4.11)) against that produced using (4.16) for $N=8$, $16$, $32$, $64$, $128$ as $\Delta t=0.01$, $0.002$, and $0.0004$ respectively in Fig. 4. \smallskip \centerline{\epsfbox{outcmp8b1.ps} \epsfbox{outcmp16b1.ps}}\smallskip \centerline{\epsfbox{outcmp32b1.ps} \epsfbox{outcmp64b1.ps}}\smallskip \centerline{\epsfbox{outcmp128b1.ps}}\smallskip \centerline{Fig. 4 Accuracy Check for for Concentric Circles} Fig.~4 reveals that our simulation becomes more accurate as the time step decreases and as more mesh points are used, which is to be expected. To make a detailed comparison, first we fixed a time (0.3 in the figures below) and for each value of $\Delta t$ estimated the ``spatial error'', $se$, as the difference between the curvature of the small circle when 128 mesh points were used and that when $N(=8, 16, 32, 64)$ points were used. Then assuming the spatial error satisfies the relationship $$ se = c{\left({1\over N}\right)}^\alpha,\eqno(4.12)$$ for two choices, $N_1$ and $N_2$, we have $$ \alpha={{\ln\left({se_1\over se_2}\right)}\over {\ln\left({N_2\over N_1}\right)}}.\eqno(4.13)$$ This allows us to estimate the spatial convergence rate $\alpha$ according to (4.13): \settabs\+\indent&(n1/n2)/$\Delta t$\qquad&2.275846725508301\quad&2.2758 46725508301\quad&2.275846725508301\quad&\cr \vbox{\hrule width4.8in height1pt} \+&(N1/N2)/$\Delta t$&0.01&0.002&0.0004\cr \vbox{\hrule width4.8in height1pt} \+&8/16&2.173 &2.167 &2.165\cr \+&16/32&2.121 &2.119 &2.119\cr \+&32/64&2.341 &2.340 &2.340\cr \vbox{\hrule width4.8in height1pt} \centerline{Table 1 Spatial Convergence Rate Estimation} Similarly, compared with the data in the case of $\Delta t=0.0004$, we calculated the relative errors of the other cases and estimate the temporal convergence rate: $$\vbox{ \settabs\+\indent&($n/(\Delta t_1/\Delta t_2)$\qquad&\cr \vbox{\hrule width2.5in height1pt} \+&$N/(\Delta t_1/\Delta t_2)$&0.01/0.002\cr \vbox{\hrule width2.5in height1pt} \+&8&1.091\cr \+&16&1.091\cr \+&32&1.091\cr \+&64&1.091\cr \+&128&1.092\cr \vbox{\hrule width2.5in height1pt} }$$ \centerline{Table 2 Temporal Convergence Rate Estimation} Tables 1 and 2 suggest that for the algorithm given by (4.16), the spatial and temperal convergence rates are of orders around 2 and 1, respectively. Now we report the results of implementing the algorithm with a variety of initial curves $\Gamma_0$. Since a component of the interface accelerates as it shrinks to a point, we make a dynamical choice of time step $\Delta t$ according to the maximum interfacial velocity $$\Delta t = \Delta t_0 *\min\left\{{1\over V_{max}}, 1\right\},$$ where $\Delta t_0$ is the initial setting of the time step and $V_{max}= \max(|V_i|)$, for $i=1$, $\dots$, $MN$ at each time $t$. We first consider the relaxation of a single body as it becomes circular. The initial shape is given parametrically by $x=(2+0.5\sin3\theta)\cos\theta$, $y=(2+0.5\sin3\theta)\sin\theta$. $N=132$ mesh points are distributed around the boundary. The interface shape and curvature as a function of arc length at four subsequent times are displayed in Fig.~5, where the arc length is measured from the point where the interface intersects the positive $x$ axis. \settabs 2\columns \+\epsfbox{outinrose132.ps} & \epsfbox{outkrose132c0.ps}\cr \medskip \+\epsfbox{outfrose132c1.ps} & \epsfbox{outkrose132c1.ps}\cr \medskip \+\epsfbox{outfrose132c2.ps} & \epsfbox{outkrose132c2.ps}\cr \medskip \+\epsfbox{outfrose132c3.ps} & \epsfbox{outkrose132c3.ps}\cr\smallskip \centerline{Fig. 5 Interface Evolution for Rose Curve} Notice that the initially nonconvex body becomes convex fairly rapidly and then remains convex with the interfacial velocity decreasing as the boundary becomes closer to being circular. It was proved in [5] that a curve sufficiently close to being circular converges to a circle. It would be natural to conjecture that convexity is preserved by this evolution, however, the next example suggests this is not the case. In [16] Mayer proved that convexity could be lost when the evolution was governed by the one phase flow. The initial curve in [16] is basically that of our next experiment. We consider a shape given by a long thin tube with semicircular end caps. The straight part has length of 16, while the radius of the circular part is 0.125. We take $N=120$ mesh points around the boundary. Shown in Fig. 6 are the interface shape and curvature as a function of arc length at five subsequent times, where the arc length is measured from the point (8. 0163, -0.1239) (the body center is located at (0, 0)). \+\epsfbox{outgintube120.ps}& \epsfbox{outktube120c0.ps}\cr \smallskip \+\epsfbox{outftube120c1.ps}& \epsfbox{outktube120c1.ps}\cr \smallskip \+\epsfbox{outftube120c2.ps} & \epsfbox{outktube120c2.ps}\cr \smallskip \+\epsfbox{outftube120c3.ps} & \epsfbox{outktube120c3.ps}\cr \smallskip \+\epsfbox{outftube120c4.ps}& \epsfbox{outktube120c4.ps}\cr \smallskip \centerline{Fig. 6 Interface Evolution for a Thin Tube}\medskip Convexity is lost immediately as the ends become bulbous but eventually convexity is regained and the curve tends to a circle. A careful analysis shows that even though convexity is lost, the initial motion widens the tube in the center but it widens more quickly near the endcaps. This raises the question of whether or not a curve can ever pinch off, increasing the number of components. Our next simulation considers an initial curve which is a dumbbell formed by taking two ellipses joined by a narrow tube. If the ellipses are small then the curve does not form self-intersections. However, for large ellipses, such as shown in Figure 7, self intersections are formed. Of course, the model ceases to be valid as soon as singularities form but our simulation suggests that the model allows pinching-off. Even though geometric singularities occur, it can be shown that no singularities form in the integral when one part of the evolving curve becomes tangent to another part. However, for computational ease $10^{-8}$ is added to the terms $|a_1| - |a_4|$ in (3.7), thereby avoiding possible problems arising due to discretization. Figure 7 superimposes the initial curve and its state at two later times (smooth curves). It is the nonlocal nature of the problem that allows the curvature to drive a curve to self-intersect. It is well known that the explicit motion by curvature in the plane does not allow self intersections and also preserves convexity. For the case at hand it is interesting to observe the undulations in the connecting tube prior to the occurance of self intersection. One is led to speculate whether simultaneous multiple self intersections could occur in the connecting tube. \epsfbox{d10d2.ps} \centerline{Fig. 7 Pinch-Off} \epsfbox{Fig82.ps} \centerline{Fig. 8 Coalescence} The reverse change in topology, coalescence of two particles, also seems possible as shown by the simulation starting with two equal ellipses having vertices close together with colinear minor axes. The symmetry ensures that neither particle gains area at the expense of the other, so one expects the particles to evolve to become equal circles, which is a stationary configuration. However, the particles are so close that these circles could not exist without overlapping unless their centers are further apart than the centers of the original ellipses. In fact the centers do move apart but not far enough (see Figure 8). Figures 9 and 10 show the morphological evolution of four particles in an initial configuration considered in [21], but here we use our algorithm for the two phase flows, while in [21] the evolution was according to single phase diffusion. The initial radii, from left to right, are 1, 0.9, 0.8 and 0.8. The centers of the particles are located at (-2.5, 0), (0, 0), (2, 1.2) and (2, -1.2), respectively. $N=64$ mesh points are taken equally spaced around the boundary of each particle. Originally, the center of the overall mass of the particles is at (0.019, 0). In Fig. 9, we use a small diamond to indicate the center of mass as the particles evolve. \+\epsfbox{outgin4c64_1.ps}& \epsfbox{outf4c64_1c1.ps}\cr \+\epsfbox{outf4c64_1c2.ps}& \epsfbox{outf4c64_1c3.ps}\cr \+\epsfbox{outf4c64_1c4.ps}& \epsfbox{outf4c64_1c5.ps}\cr \centerline{Fig. 9 Evolution of Four Particles}\smallskip Notice that the leftmost (largest) particle grows slowly with time, maintaining approximately circular symmetry, while the middle particle expands more quickly accompanied by the shrinking of the smaller particles. This is because there is a strong local diffusional interaction between the middle particle and nearby small particles. With a dynamic change in the time step, our algorithm is able to follow to the disappearance of the smallest particles. The removal of vanishing particles is done by deleting particles whose average curvature is greater than 300. It is interesting to note the motion of the center of mass during the evolution. For the single phase problem, one can show that the centroid remains fixed. Fig.~10 gives a clearer description of shape distortions before the smaller particles disappear. It is clear that the middle particle recedes slightly on its left side, material being transported to the leftmost particle, which has a corresponding advance on its right side. By the time the smallest particles disappear, the middle particle has become largest and so this is the only survivor (see Figure 9). \bigskip \centerline{\bf Appendix }\medskip This contains two programs. The first, called {\tt calcinit.c} constructs initial curves and after compiling, execution is done by typing {\tt calcinit file1 file2}, where the user gives names to {\tt file1} and {\tt file2}. The second, called {\tt hs.c} executes our algorithm. It is executed with the command {\tt hs file1 file3 file4 file5 file6}. In constructing the curve the user is asked to choose the number of curves, their type from a menu (e.g. ellipse, tube, etc.), and then dimensions (e.g. semiminor and semimajor axes) and center locations are required. When executing {\tt hs}, the desired values of $N, \Delta t$ and the number of time steps are requested. The output is placed in {\tt file4} and can be viewed using gnuplot or some other data plotting program. \epsfbox{outf4c64_1cal.ps} \centerline{Fig. 10 Evolution of Four Particles} \bigskip \centerline{\bf References }\medskip \item{[1]} Alikakos, N.D., Bates, P.W. and Chen, X., ``Convergence of the Cahn-Hilliard equation to the Hele-Shaw Model,'' {\it Arch. Rat. Mech. Anal.} {\bf 128} (1994), 165--205. \item{[2]} Almgren, R., Singularity formation in Hele-Shaw bubbles, Preprint May 1995. \item{[3]} Baker, G.R., Meiron, D.I., and Orszag, S.A., ``Boundary integral methods for axisymmetric and three-dimensional Rayleigh-Taylor instability problems,'' {\it Physica D} {\bf 12} (1984), 19--31. \item{[4]} Baker, G.R., ``Generalized vortex methods for free-surface flows,'' in Meyer, R.E. (ed.), Waves on fluid interfaces, Academic Press, New York, (1983), 53--81. \item{[5]} Chen, X., ``The Hele-Shaw problem and area-preserving, curve shortening motions,'' {\it Arch. Rat. Mech. Anal.} {\bf 123} (1993), 117--151. \item{[6]} Chen, X., Hong, J., Yi, F., ``Existence, uniqueness, and regularity of classical solutions of the Mullins-Sekerka problem,'' preprint. \item{[7]} Chen, X., Hou, T.Y., and Zhu, J., Preprint, 1995. \item{[8]} Constantin, P., and Pugh, M., ``Global solutions for small data to the Hele-Shaw problem,'' {\it Nonlinearity} {\bf 6} (1993), 393--415. \item{[9]} DeGregoria, A.J., and Schwatz, L.W., ``A boundary-integral method for two-phase displacement in Hele-Shaw cells,'' {\it J. Fluid Mech.} {\bf 164} (1986), 383--400. \item{[10]} Deng, X., ``A numerical analysis approach for the Hele-Shaw problem,'' M.S. Thesis, Brigham Young University, April 1994. \item{[11]} Folland, G. B., Introduction to partial differential equations, Princeton University Press, Princeton, NJ, 1976. \item{[12]} Hou, T.Y., Lowengrub, J.S., and Shelley, M.J., Removing the stiffness from interfacial flows with surface tension, {\it J. Comp. Phys.} {\bf 114} (1994), 312--338. \item{[13]} Kellogg, O.D., Foundations of potential theory, Springer-Verlag, New York, 1929. \item{[14]} Kessler, D.A., Kopik, J., and Levine, H., ``Numerical simulation of two- dimensional snowflake growth,'' {\it Phys. Rev. A} {\bf 30} (1984), 2820--2823. \item{[15]} Kessler, D.A., and Levine, H., ``Theory of the Saffman-Taylor `finger' pattern. I \& II,'' {\it Phys. Rev. A} {\bf 33} (1986), 2621--2639. \item{[16]} Mayer, U.F., ``One-sided Mullins-Sekerka flow does not preserve convexity,'' {\it Electr. J. Diff. Eqns.} {\bf 1993} (1993), No. 8, 1--7. \item{[17]} McFadden, G.B., Voorhees, P.W., Boisvert, R.F., and Meiron, D.I., ``A boundary integral method for the simulation of two-dimensional particle coarsening,'' {\it J. Sci. Comp.} {\bf 1} (1986), 117--143. \item{[18]} Menikoff, R., and Zemach, C., ``Methods for numerical conformal mapping,'' {\it J. Comp. Phys.} {\bf 36} (1980), 366--410. \item{[19]} Mullins, W.W. and Sekerka, R.F., ``Morphological stability of a particle growing by diffusion and heat flow,'' {\it J. Appl. Phys.} {\bf 34} (1963), 323--329. \item{[20]} Pego, R., ``Front migration in the nonlinear Cahn-Hilliard equation,'' {\it Proc. Roy. Soc. London A} {\bf 422} (1989), 261--278. \item{[21]} Voorhees, P.W., McFadden, G.B., and Boisvert, R.F., ``Numerical simulation of morphological development during Ostwald ripening,'' {\it Acta Meta.} {\bf 36} (1988), 207--222. \item{[22]} Voorhees, P.W., ``Ostwald ripening of two-phase mixtures,'' {\it Annu. Rev. Mater. Sci.} {\bf 22} (1992), 197--215. \bigskip \vbox{Peter W. Bates \hfil\break Mathematics, Brigham Young University \hfil\break Provo, UT 84602 \hfil\break E-mail address: peter@math.byu.edu} \bigskip \vbox{Xinfu Chen \hfil\break Mathematics, University of Pittsburgh \hfil\break Pittsburgh, PA 15260 \hfil\break E-mail address: xinfu$+$@pitt.edu} \bigskip \vbox{Xinyu Deng \hfil\break Mathematics, Brigham Young University \hfil\break Provo, UT 84602 \hfil\break E-mail address: cindy@math.byu.edu} \vfil\eject \centerline{\bf Erratum} \medskip {\bf September 26, 1995}. I did not write Theorem 2.3 correctly so I enclose the following new version with the necessary changes in the first few lines of the proof.\hfil\break Best, Peter Bates. \bigskip \noindent {\bf Theorem 2.3} {\sl If $\Gamma_{0,T}\equiv\cup_{0\leq t\leq T}\left({\Gamma_t\times\{t\}}\right)$ is a continuous family of $C^3$ curves which satisfy (2.3) and (2.4) with $\Gamma=\Gamma_t$ for $g=V(x,t)$, the normal velocity of $\Gamma_t$, and some $c=c(t)$, where $f=K(x,t)$ is the curvature of $\Gamma_t$ at $x$, then $\Gamma_{0,T}$ is the interface associated with the solution to (2.1). Conversely, if $(u,\Gamma_{0,T})$ is a solution to (2.1) then (2.3)--(2.4) hold for each $t\in [0,T]$ with $\Gamma=\Gamma_t, g=V$ and $f=K$.} \noindent {\bf Proof:} Let $(\Gamma_t, g, c)$ be a solution to (2.3)--(2.4) at each time $t\in [0,T]$ with $f=K(\cdot, t)$ and suppose that $g = V$, the normal velocity of $\Gamma_t$. Then Lemma 2.1 shows that defining $u(\cdot, t)$ on ${\Bbb R}^2\backslash\Gamma_t$ by $$u(x,t)={1\over 2\pi} \int_{\Gamma_t} \ln\vert x-y\vert g(y,t)ds_y+c(t)$$ gives a solution to (2.1) with $-\left[{{\partial c\over \partial n}}\right]_{\Gamma_t}$ being given by $V$. \bye