\documentclass[reqno]{amsart}
\usepackage{graphicx}
\AtBeginDocument{{\noindent\small
{\em Electronic Journal of Differential Equations},
Vol. 2002(2002), No. 95, pp. 1--29.\newline
ISSN: 1072-6691. URL: http://ejde.math.swt.edu or http://ejde.math.unt.edu
\newline ftp ejde.math.swt.edu (login: ftp)}
\thanks{\copyright 2002 Southwest Texas State University.}
\vspace{1cm}}
\begin{document}
\title[\hfilneg EJDE--2001/95\hfil Heteroclinic orbits]
{Heteroclinic orbits, mobility parameters and stability for thin film
type equations}
\author[Richard. S. Laugesen \& Mary C. Pugh\hfil EJDE--2001/95\hfilneg]
{Richard. S. Laugesen \& Mary C. Pugh}
\address{Richard. S. Laugesen \hfill\break
Department of Mathematics, University of Illinois, Urbana, IL 61801, USA}
\email{laugesen@math.uiuc.edu}
\address{Mary C. Pugh\hfill\break
Department of Mathematics, University of Toronto, Toronto, ON M5S 3G3, Canada}
\email{mpugh@math.toronto.edu}
\date{}
\thanks{Submitted February, 28, 2002. Published November 5, 2002.}
\subjclass[2000]{35K55, 37C29 ,37L15, 76D08}
\keywords{Nonlinear PDE of parabolic type, heteroclinic orbits, \hfil\break\indent
stability problems, lubrication theory}
\begin{abstract}
We study the phase space of the evolution equation
\[
h_t = -(h^n h_{xxx})_x - \mathcal{B} (h^m h_x)_x ,
\]
where $h(x,t) \geq 0$. The parameters $n>0$, $m \in \mathbb{R}$,
and the Bond number $\mathcal{B}>0$ are given.
We find numerically, for some ranges of $n$ and $m$, that
perturbing the positive periodic steady state in a certain direction
yields a solution that relaxes to the constant steady state. Meanwhile
perturbing in the opposite direction yields a solution that appears to
touch down or `rupture' in finite time, apparently approaching a
compactly supported `droplet' steady state.
We then investigate the structural stability of the evolution by
changing the mobility coefficients, $h^n$ and $h^m$. We find evidence
that the above heteroclinic orbits between steady states are perturbed
but not broken, when the mobilities are suitably changed.
We also investigate touch--down singularities, in which the solution
changes from being everywhere positive to being zero at isolated
points in space. We find that changes in the mobility exponent $n$
can affect the {\em number} of touch--down points per period,
and affect whether these singularities occur in finite or infinite time.
\end{abstract}
\maketitle
\numberwithin{equation}{section}
\newcommand{\full}[2]{\frac{d #1}{d #2}}
\section{Introduction}
\label{introduction}
We study the evolution equation
\begin{equation} \label{evolve}
h_t = - (h^n h_{xxx})_x - \mathcal{B} (h^m h_x)_x ,
\end{equation}
%
where $n>0, m \in \mathbb{R}$, and the Bond number $\mathcal{B}>0$. This is the
one-dimensional version of
\begin{equation} \label{evolve2}
h_t = - \nabla \cdot( f(h) \nabla \Delta h) - \nabla \cdot (g(h) \nabla h)
\end{equation}
%
with $f(h) = h^n, g(h) = -\mathcal{B} h^m$. Such equations have been used to
model the dynamics of a thin film of viscous liquid. The air/liquid
interface is at height $z = h(x,y,t) \geq 0$ and the liquid/solid
interface is at $z=0$. The one dimensional equation (\ref{evolve})
applies if the liquid film is uniform in the $y$ direction.
The fourth order term in equation (\ref{evolve2}) reflects surface
tension effects \cite{myers,oron}. Typically $f(h) \sim h^n$ as $h \rightarrow 0$, for some $1 \leq n \leq 3$, and this motivates our choice of
$f(h)=h^n$ in (\ref{evolve}). Notice also that the fourth order term $-(h^n
h_{xxx})_x$ in (\ref{evolve}) is linearly stabilizing around the
constant steady state.
The second order term in (\ref{evolve2}) can reflect gravity, van der
Waals interactions, thermocapillary effects or the geometry of the
solid substrate. Typically $g(h) \sim \mathcal{B} h^m$ as $h \rightarrow 0$,
for some $m, \mathcal{B} \in \mathbb{R}$. Choosing $\mathcal{B} > 0$, the second order term $-\mathcal{B}
(h^m h_x)_x$ in (\ref{evolve}) is linearly destabilizing around the
constant steady state (it is like a backwards heat equation term). The
competition between this destabilizing term and the stabilizing fourth
order term generates interesting dynamics. The dynamics are less
interesting when $\mathcal{B} \leq 0$, which we do not consider in this paper.
\begin{figure}[ht]
\begin{center}
\includegraphics[width=.4\textwidth]{steady_states.ps}
\end{center}
% \renewcommand{\figurename}{Fig.$\!$$\!$} \setcaptionwidth{5in}
\caption{\label{steady} Four types of steady state. The $x$-axis
is space. The two steady states on the left extend smoothly to be
$P$-periodic. The two on the right have less regularity and
are called `droplet' steady states. The third profile has zero
contact angles and the fourth has nonzero contact angles. Both
have length less than $P$, and so are possible long--time
limits of a $P$-periodic solution of the evolution equation.}
\label{four-steady}
\end{figure}
\subsection*{Definitions}
%
In \cite{LP1} we studied the four types of steady state shown in
Figure~\ref{four-steady}. There are two smooth types of steady state:
constant and (nonconstant) positive periodic. Two other types of steady state with lower regularity are the `droplet' steady states. They either have zero contact angles or nonzero contact angles.
Positive periodic steady states are classical solutions of the steady
state equation (that is, (\ref{evolve}) with $h_t \equiv 0$). When we
refer to the period, $P$, of a nonconstant positive periodic steady
state, we mean the shortest period. The area, $A$, is then defined
to be $\int_0^{P} h_{\rm ss} \; dx$. Constant steady states have no
shortest period, however we will always discuss them in the context of
an initial value problem, for which the period and area are
unambiguous. Because constant steady states are periodic steady
states of the most trivial kind, in the following when we refer to
positive periodic steady states we implicitly mean {\it nonconstant}
ones.
A droplet steady state $h_{\rm ss}$ is by definition positive on some interval $(a,b)$ and zero elsewhere, with $h_{\rm ss} \in C^1[a,b]$; $h_{\rm ss}$ satisfies the steady state equation classically on $(a,b)$, and has equal acute contact
angles: $0 \leq h_{\rm ss}'(a) = - h_{\rm ss}'(b) < \infty$. (Throughout the
paper, if a function has only one independent variable then we use
${}'$ to denote differentiation with respect to that variable:
$h_{\rm ss}' = (h_{\rm ss})_x$.) The area, $A$, of a droplet is $A = \int_a^b
h_{\rm ss} \; dx$ and the length is $b-a$.
As suggested by Figure~\ref{four-steady}, we are interested in droplet
steady states whose length is equal to or shorter than $P$. This
is because we want to study the initial value problem with
positive initial data of period $P$ and area $A$. A steady
droplet with area $A$ whose length is less than or equal to $P$
must be considered as a possible long--time limit of the evolution. A
`configuration' of droplet steady states is defined to be a collection
of steady droplets whose supports are disjoint. Any such
configuration whose total area is $A$ and whose total length is
less than $P$ must also be considered as a possible long--time
limit of the initial value problem.
\subsection*{Background, and Overview of Paper}
%
In \cite{LP2} we proved linear stability results for steady states.
Section~\ref{bifurcation} summarizes these results by means of a
family of bifurcation diagrams, and also presents a weakly nonlinear stability analysis.
Our investigations in the rest of the paper are guided by the study of a dissipated energy for the evolution equation. This energy is defined for $P$-periodic functions
$\ell$ on $\mathbb{R}$ to be
%
\begin{equation} \label{energy_is}
\mathcal{E}(\ell) = \int_0^{P} \left[ \frac{1}{2} (\ell'(x))^2
- \frac{\mathcal{B}}{(m-n+1)(m-n+2)} \ell^{m-n+2}(x)
\right] dx.
\end{equation}
%
This energy is strictly dissipated: if $h(x,t)$ is a smooth, spatially
$P$-periodic solution of (\ref{evolve}) then $\frac{d\ }{dt}
\mathcal{E}(h(\cdot,t)) \leq 0$, with equality if and only if $h$ is a steady
state of the evolution (cf.\ \cite[\S2.1]{LP3}). Like the evolution
equation, the energy $\mathcal{E}(\ell)$ is insensitive to translation in $x$.
The evolution (\ref{evolve}) describes {\em gradient flow} for the energy $\mathcal{E}$, with
respect to the following weighted $H^{-1}$ inner product. Let $h(x,t)$
be a positive smooth function that is $P$-periodic in $x$. For each $t$,
define an inner product on the space of continuous $P$-periodic
functions with mean value zero, by
%
\begin{equation} \label{inner_prod}
\langle u,v \rangle := \int_0^{P} U'(x) V'(x) h(x,t)^n \, dx ,
\end{equation}
%
where $U$ and $V$ are $P$-periodic, have zero mean, and are
determined from $u$, $v$, and $h(\cdot,t)$ via: $u=(h^n U_x)_x$ and
$v=(h^n V_x)_x$. Then the evolution equation (\ref{evolve}) is
equivalent to:
\[
\frac{\delta \mathcal{E}}{\delta \phi} =
\lim_{s \to 0} \frac{\mathcal{E}(h + s \phi) - \mathcal{E}(h)}{s}
= - \langle h_t , \phi \rangle \quad
\]
for all continuous $P$-periodic $\phi$ having mean value zero.
The variation of the energy is most negative in the direction $\phi =
h_t$, so that the evolution equation for $h(x,t)$ describes flow by
steepest descent on the energy surface of $\mathcal{E}$, with respect to the
time- and solution-dependent inner product (\ref{inner_prod}). The
steady states are critical points of this energy surface, with
unstable steady states being saddle points and asymptotically stable
steady states being minima. The above gradient flow formulation was
observed by Fife \cite{FifeGradOrig} and by Taylor and Cahn \cite{TC}.
Gradient flow ideas for related equations have been used in
\cite{BBW,WB00,Otto}. We refer the readers to Fife's survey article
on pattern formation in gradient systems \cite{FifeReview}.
% Also, while the energy $\mathcal{E}$ depends on the ratio of the coefficients,
% $g(h)/f(h) = h^{m-n}$, only $f(h) = h^n$ is relevant for the
% %evolution on the energy landscape.
Note that the perturbing function $\phi$ is required to have mean value zero
because the equation (\ref{evolve}) preserves area, under
periodic (or Neumann) boundary conditions. That is, given initial data
with area $A$, only those points on the energy landscape
representing functions with the {\em same} area will be accessible
to the solution.
In \cite{LP3} we studied the energy landscape of (\ref{energy_is}) under the evolution (\ref{evolve}), and
conjectured the existence of heteroclinic orbits from certain
high-energy steady states to certain low-energy steady states. In
this paper we find such orbits numerically, in Section~\ref{pert_hets}, by computing solutions of the initial value problem for (\ref{evolve}).
We further ask how robust these orbits are under changes in the
mobility coefficients $h^n$ and $h^m$. If we vary $m$ and $n$ while
keeping $m-n$ fixed, the steady states and their energy stability are
unchanged (since $\mathcal{E}$ depends only on $m-n$). We find that the
heteroclinic orbits are perturbed but do not break, when $m$ and $n$
are varied in this way.
In Sections \ref{time_pinch} and \ref{one_two_sec} we study orbits
from a positive periodic steady state to a droplet steady state. Here
the long--time limit is not strictly positive, raising the question of
whether or not the solution will be positive (and hence classical) for all
time, or whether there will be `touch--down singularities', places and times at which the solution is zero (and hence weak). We find that
the value of $m-n$ can affect whether touch--down singularities are
present (\S\ref{time_pinch}), as well as the number of such
singularities that arise per period (\S\ref{one_two_sec}).
Section~\ref{simulations} presents a detailed numerical study of the
evolution equation (\ref{evolve}), for initial data close to a steady
state and for a number of different exponents $n$ and $m$. Our
stability and energy level results from \cite{LP2,LP3} lead to many
predictions for the behavior of the solutions, both short and long
time, and these predictions are borne out by our simulations. However there are also situations for which we can make no prediction, and where the numerically observed behavior is rather intriguing.
In Section~\ref{conclusions} we sum up the
paper. Appendix~\ref{numerics} discusses our numerical methods and
gives the specific parameters used in the simulations.
% \section{Review of steady states}
% \label{power_law_review}
%
% We quickly review some basic facts about positive periodic steady
% states. For more on steady states and their properties, and for
% justifications of the following remarks, one can consult
% \cite[\S2.3]{LP3} and \cite{LP1}.
%
% \vskip 6pt
%
% We start with a non-constant, positive, $P$-periodic steady state
% $h_{\rm ss} \in C^4(\mathbb{R})$ of the evolution equation (\ref{evolve}). One can
% integrate up the steady state equation to obtain
% \begin{equation} \label{ss1}
% h_{\rm ss}''' + \mathcal{B} h_{\rm ss}^{q-1} h_{\rm ss}' = 0 ,
% \end{equation}
% where
% \[
% \fbox{$q := m-n+1$.}
% \]
% A second integration shows the steady states satisfy a nonlinear
% oscillator equation, when $q\neq 0$:
% \begin{equation} \label{ss3}
% h_{\rm ss}'' + \frac{\mathcal{B} h_{\rm ss}^q - D}{q} = 0
% \end{equation}
% for some constant $D$. (For $q=0$ the analogous equation is $h_{\rm ss}'' +
% \mathcal{B} \log (h_{\rm ss}/D) = 0$.) Rescaling $h_{\rm ss}$ and $x$ leads to
% \begin{eqnarray} \label{ss4}
% k'' + \frac{k^q - 1}{q} & = & 0, \quad q \neq 0, \\
% \label{ss5}
% k'' + \log k & = & 0, \quad q = 0 .
% \end{eqnarray}
%
% Every positive periodic steady state $h_{\rm ss}$ can, after translation, be
% rescaled to such a function $k=k_\alpha$, where $\alpha=k(0) \in
% (0,1)$ and $k_\alpha'(0)=0$. Conversely, for each $q \in \mathbb{R}$ and
% $\alpha \in (0,1)$ there exists a unique smooth positive periodic
% $k_\alpha$ satisfying equations (\ref{ss4}--\ref{ss5}) and with
% $k_\alpha(0)=\alpha, k_\alpha'(0)=0$. The same holds for $\alpha=0$,
% when $q>-1$, although $k_0$ may be only $C^1$-smooth at the $x$
% intercepts (see \cite[Theorem~3.2]{LP1}).
%
% We write
% \[
% P_\alpha=P(\alpha) \quad \text{and} \quad
% A_\alpha=A(\alpha)=\int_0^{P(\alpha)} k_\alpha(x) \, dx
% \]
% respectively for the least period of $k_\alpha$ and for the area under
% its graph. As seen in \cite{LP2,LP3}, the function
% \[
% E(\alpha) := P(\alpha)^{3-q} A(\alpha)^{q-1}
% = P(\alpha)^2 [A(\alpha)/P(\alpha)]^{q-1}
% \]
% plays a large role in the stability theory of the steady states, in
% part due to its rescaling invariance:
% \begin{equation} \label{invariant}
% E := \mathcal{B} P^{3-q} A^{q-1} = P(\alpha)^{3-q} A(\alpha)^{q-1}
% = E(\alpha),
% \end{equation}
% where $P$ and $A$ are the period and area of the original steady
% state $h_{\rm ss}$.
\section{Bifurcation diagrams and weakly nonlinear analysis}
\label{bifurcation}
First, we rescale the equation and present the weakly nonlinear
analysis near the constant steady state. Then, we present bifurcation
diagrams that summarize the linear stability of constant steady states
and positive periodic steady states.
\subsection{Non-dimensionalizing the equation}
\label{non-dimensional}
A solution of the evolution (\ref{evolve}) reflects five free
parameters: $m$, $n$, $\mathcal{B}$, the period $P$, and the area $A$.
First, we rescale space, time, and the solution itself:
\[
\zeta = \frac{x}{P},
\quad
t' = P^{n+4} A^{-n} t,
\quad \mbox{and} \quad
\frac{A}{P} \; \eta (\zeta, t' )
= h(x,t).
\]
The rescaled solution, $\eta$, has period $1$ and area $1$ and
satisfies the
rescaled evolution equation
\begin{equation} \label{evolve_rescaled}
\eta_{t'} = - ( \eta^n
\eta_{\zeta \zeta \zeta})_\zeta - E (\eta^m \eta_\zeta)_\zeta
\end{equation}
where
\begin{equation} \label{bifurc_param}
E = \mathcal{B} A^{q-1} P^{3-q}
\quad \mbox{and} \quad
q := m-n+1.
\end{equation}
%
%The energy $\mathcal{E}$ was defined for $P$-periodic functions in
%(\ref{energy_is}). If $\ell$ is a $P$-periodic function with area
%$A$, then rescaling it to $\rho$, a $1$-periodic function with area
%$1$, rescales $\mathcal{E}$ as follows:
%
%\begin{equation} \label{energy_rescaled}
%A^2 P^{-3} \; \mathcal{E}_{E}(\rho) :=
%A^2 P^{-3} \int_0^1 \left[ \frac{1}{2} (\rho'(x))^2
%- \frac{E}{q(q+1)} \rho^{\; q+1}(x) \right] \; dx
%\end{equation}
%The energy $\mathcal{E}$ for the unrescaled equation contains $\mathcal{B}$ in its
%integrand and the unrescaled evolution equation (\ref{evolve}) can be
%viewed as a gradient descent flow with respect to it. The rescaling
%produces a constant multiple $A^2 P^{-3}$ and the integrand of
%(\ref{energy_rescaled}) now contains the bifurcation parameter $E$.
%For this reason, in the study of the unrescaled equation, it is
%natural that it is the particular combination of $\mathcal{B}$, $P$, and
%$A$ that makes $E$ that is crucial in the study of the steady
%states \cite{LP1}, their linear stability \cite{LP2}, and their energy
%$\mathcal{E}$ \cite{LP3}. Given initial data with period $P$ and area
%$A$, we use $P$, $A$, and $\mathcal{B}$ to find the value of $E$.
%In turn, from $E$ we have information about the energy landscape,
%allowing us to make predictions about the behavior of the solution of
%the initial value problem (\ref{evolve}) when the initial data is near
%a steady state. In Section \ref{simulations} we discuss these
%predictions and present numerical simulations confirming the
%predictions.
%
In the simulations that follow, in Sections~\ref{mobility} and \ref{simulations}, we always take $P = A = 1$. Then the
original evolution equation (\ref{evolve}) and the rescaled evolution
(\ref{evolve_rescaled}) are identical, and $\mathcal{B} = E$ is the
bifurcation parameter. So in the following we will
refer to $h$, $\mathcal{B}$ and equation (\ref{evolve}), rather than to
$\eta$, $E$ and equation (\ref{evolve_rescaled}).
\subsection{Weakly nonlinear analysis}
For the weakly nonlinear analysis, we consider values of $\mathcal{B}$ such
that the constant steady state has one mode which is barely linearly
stable or is barely linearly unstable. Since all other modes are
strongly damped, this provides a separation of timescales, allowing
one to find a reduced representation of the PDE in terms of an ODE
governing the amplitude of the unstable mode (cf.\
\cite[\S5.1]{Manneville}).
Linearizing equation (\ref{evolve}) about $\overline{h} = 1$ and considering
perturbations $\epsilon \cos(k 2 \pi x + \phi)$, one finds a critical
value $k_c = \sqrt{\mathcal{B}}/(2 \pi)$. If $k_c \leq 1$ then there are no
unstable modes. If $k_c > 1$ then there is a finite collection of
linearly unstable modes. We assume $\mathcal{B} = 4 \pi^2 + Q \delta^2$ where
$Q = \pm 1$. Here, $\delta$ is a small parameter and varying $\delta$
results in $k_c$ moving through the wave number $1$. We then
introduce a slow time-scale $\tau = \delta^2 t$ and expand the
solution in orders of $\delta$: $ h(x,\tau) = 1 + \delta \; h_1(x,\tau) +
\delta^2 h_2(x,\tau) + O(\delta^3).$ For simplicity, we assume the
solution is even. By the usual arguments, $h_1(x,\tau) = A(\tau)
\cos(2 \pi x)$ and $h_2(x,\tau) = B(\tau) \cos(4 \pi x)$. Putting
this ansatz into the evolution equation (\ref{evolve}) and expanding
in orders of $\delta$, we find there are no $O(1)$ or $O(\delta)$
terms. At $O(\delta^2)$ one determines the amplitude $B(\tau)$ in
terms of $A(\tau)$. At $O(\delta^3)$ one finds that $A(t)$ satisfies
\[
\full{A}{\tau} = 4 \pi^2 Q A(\tau) - \kappa A(\tau)^3 ,
\quad \mbox{where} \quad
\kappa = \frac{8}{3} \pi^4 (q-1) (7/4 - q )
\]
and we recall $q=m-n+1$. The dynamics of the amplitude $A(\tau)$ depend on the signs of the
Landau constant $\kappa$ and of the linear term. If $\kappa > 0$ then
for $Q = 1$ the steady state $A_0(\tau) \equiv 0$ is linearly unstable
and $A(\tau)$ saturates to the linearly stable steady state $A_c(\tau)
\equiv \sqrt{\sigma/\kappa}$. This corresponds to a supercritical
bifurcation. If $\kappa < 0$ then for $Q = -1$ the steady state
$A_0(\tau) \equiv 0$ is linearly stable and the steady state
$A_c(\tau) \equiv \sqrt{-\sigma/\kappa}$ is linearly unstable. This
corresponds to a subcritical bifurcation. And so,
\[
1 < q < 7/4 \Longrightarrow \mbox{{\it supercritical
bifurcation}},
\quad
q < 1 \; \mbox{or} \; 7/4 < q \Longrightarrow \mbox{{\it subcritical
bifurcation}}.
\]
Subcritical bifurcations are often seen in systems that can have
finite-time pinching (rupture) singularities {\it e.g.}
\cite[\S3.2]{BBW}, \cite[\S{IV}]{OB}. The above weakly nonlinear
analysis was done for $q=-3$ in \cite[\S2.3]{WB00}.
\subsection{Bifurcation diagrams}
Figure~\ref{bifurcs} gives bifurcation diagrams for representative
$q$-values. To construct a bifurcation diagram, we fix a value for
$q=m-n+1$ and then compute a collection of positive periodic steady states
$h_{\rm ss}$ all with period $1$ and area $1$. Each steady state satisfies
the steady version of (\ref{evolve}) for some value of $\mathcal{B}$. We then
plot the amplitude of $h_{\rm ss}$ versus the bifurcation parameter $\mathcal{B}$. We
use our linear stability results from \cite{LP2} to determine whether
to plot with solid lines (linearly stable) or with dashed lines
(linearly unstable). The horizontal axes of these diagrams show the
linear stability of the (zero-amplitude) {\em constant} steady
state. These stability results are all with respect to zero-mean
perturbations having the same period as the steady state.
Perturbations with longer period always lead to linear instability
\cite{LP2}.
\begin{figure}[ht]
\begin{center}
\includegraphics[width=.2\textwidth]{qm3_bifurc.ps}
\hspace{.5cm}
\includegraphics[width=.2\textwidth]{q5_bifurc.ps}
\hspace{.5cm}
\includegraphics[width=.2\textwidth]{q1_bifurc.ps}
\\
%% /data/users/mpugh/Rick/Programs/Bifurcations
\includegraphics[width=.2\textwidth]{q15_bifurc.ps}
\hspace{.5cm}
\includegraphics[width=.2\textwidth]{q175_bifurc.ps}
\hspace{.5cm}
\includegraphics[width=.2\textwidth]{q176_bifurc.ps}
\\
\includegraphics[width=.2\textwidth]{q1768_bifurc.ps}
\hspace{.5cm}
\includegraphics[width=.2\textwidth]{q1775_bifurc.ps}
\hspace{.5cm}
\includegraphics[width=.2\textwidth]{q25_bifurc.ps}
\end{center}
\caption{\label{bifurcs} The horizontal axis is the bifurcation
parameter $\mathcal{B}$ and the vertical axis is the amplitude of the steady
state, $(h_{max}-h_{min})/2$. Dashed: linearly unstable; dotted:
linearly neutrally stable; solid: linearly stable. There is an
interval of $\mathcal{B}$ values for which (nonconstant) positive periodic
steady states exist. They are linearly unstable for $q \in
(-\infty,1)$, are linearly stable for $q \in (1,7/4)$, and are
linearly unstable for $q \in (q^*,\infty)$ where $q^* \approx 1.794$.
For $q \in (7/4,q^*)$, there can be two such steady states, one
linearly unstable and one linearly stable. One can prove that for $q
\leq -1$, as $B \to 0$ the amplitude of the solution converges to $3/4$.
The solid and dashed lines on the horizontal
axis represent the linear stability of the constant (zero amplitude)
steady states. }
\end{figure}
Figure \ref{bifurcs} provides nine graphs to help the reader visualize the
dependence of the diagram on $q$. (In what follows, we do not discuss
simulations for $q=1.76$ or $1.775$.) By examining the bifurcation
diagrams near the point $(4 \pi^2,0)$, one observes the subcritical
and supercritical bifurcations predicted by the weakly nonlinear
stability analysis.
The bifurcation diagrams also encode existence information for steady
states. Consider Figure~\ref{bifurcs}b for $q=1/2$. It starts at $\mathcal{B}
= 33.07$ and ends at $\mathcal{B} = 4 \pi^2$ with amplitude zero. If $\mathcal{B}$
equals $31$ or $42$, for example, then there is no (nonconstant)
positive $1$-periodic steady state with area $1$. If $\mathcal{B} = 34$ then
there is a unique nonconstant positive $1$-periodic steady state with
area $1$. The monotonicity of the bifurcation diagram corresponds to
uniqueness of the positive periodic steady state with specified period
and area, if it exists. If the diagram is non-monotonic, there may be
zero, one, or two (nonconstant) positive $1$-periodic steady states
with area $1$, depending on the value of $\mathcal{B}$. Non-monotonicity holds
for $q \in (7/4,1.794)$ where $1.794$ approximates a criticial
exponent $q^*$ (see \cite[\S5.1]{LP1}).
% Note that qualitatively, the diagrams change continuously as $q$
% increases. Also, one can show that choosing a period other than $2
% \pi$ or a Bond number other than $\mathcal{B} = 1$ simply dilates the $y$-axis
% of the diagrams.
%
% the x-axis is P_k^{3-q} A_k^{q-1} and is unaffected by the choice of
% B and P_h.
%
% the y-axis is 1/2 (k_max-k_min) (P_k/P_h)^(2/(q-1)) B^(1/(1-q)) and
% changing B and P_h corresponds to a dilation of the y-axis.
%
\section{Dynamics: the effect of changing the mobility exponents, $n$ and $m$}
\label{mobility}
Here we vary $n$ and $m$ in $h_t = - (h^n h_{xxx})_x - \mathcal{B} (h^m h_x)_x$
while keeping $m-n$ fixed. We call the exponents $n$ and $m$ mobility
parameters, since they determine the diffusion coefficients of the
fourth and second order terms in the equation. Since $q=m-n+1$ is
being kept fixed, the energy landscape and steady states of the
evolution are unchanged. The linear stability properties of the
constant and positive periodic steady states are also unchanged
\cite{LP2}. We ask three questions about the effects of changing $n$
and $m$ in this way:
%
\begin{enumerate}
\item Can a heteroclinic orbit between steady states be {\em broken},
or is it merely perturbed?
\item Can the {\em type} of a singularity be altered ({\it e.g.} from
finite-time to infinite-time)?
\item Can the {\em number} of singularities be altered ({\it e.g.}
from one to two per period)?
\end{enumerate}
%
% While one often thinks of the initial and terminal points of a
% heteroclinic connection as being {\em isolated} equilibrium points,
% our equilibria are not isolated. First, translates of a steady state
% are themselves steady states. (Of course this translational freedom
% disappears when considering even perturbations of even initial data.)
% Second, when we consider configurations of droplet (compactly
% supported) steady states, the droplets can be translated, shrunk or
% expanded as long as the total area (area) is unchanged and the
% droplets remain disjoint.
In the following, we are interested in heteroclinic orbits between
smooth steady states and between smooth steady states and droplet
steady states. For this reason, we take $q < 1$ or $q > 7/4$ since
for these $q$-values one can find linearly unstable positive periodic
steady states (see Figure~\ref{bifurcs}). For specificity, we present
results for $q = 5/2$ simulations. We observed similar phenomena for
other values of $q$. We did not have to work hard to capture the $t
\to \infty$ limit and so we expect that our infinite--time limits are not
saddle points. The reader who is curious about infinite--time limits
that are saddle points should refer to
\cite{Champneys99,Beyn02,Doedel97}.
\subsection{Perturbing heteroclinic orbits between smooth steady states}
\label{pert_hets}
First we fix $q=5/2$ and find a positive periodic steady state $h_{\rm ss}$
having period $1$ and area $1$. It is a steady state of
(\ref{evolve}) with $\mathcal{B} = 35.32$ and is linearly unstable, by
bifurcation diagram 2i. We perturb $h_{\rm ss}$, computing the solution of
$h_t = - (h h_{xxx})_x - \mathcal{B} (h^{5/2} h_x)_x$ with initial data $h_{\rm ss} +
.001 h_{\rm ss}^{\prime \prime}$. We find that the local minimum of the
solution remains fixed in space and, after a short transient, it
increases to $1$. The maximum behaves similarly, decreasing to
$1$. That is, the solution relaxes to the constant steady state as $t
\to \infty$. This suggests there is a heteroclinic orbit connecting
the positive periodic steady state to the constant steady state.
We then vary the mobility coefficients, taking $n = 0$, $1$, $2$, and
$3$, choosing $m$ in turn so that $q = m-n+1 = 5/2$. That is, we
compute solutions for the four evolution equations, all with
the
same initial data. We find that all four solutions relax to the
constant steady state: the apparent heteroclinic orbit is not
broken by this change in mobility.
But there are differences attributable to the change in mobilities: in
Figure~\ref{time_scales} we plot $h_{min}(t)$ and $h_{max}(t)$ versus
time for the four solutions. The larger the exponent $n$, the longer
it takes for the solution to relax to the constant steady state. We
explain this as follows. Because the mean value of $h(\cdot,t)$ is
$1$, $h_{min}(t)<1$ for all $t$ and $h_{max}(t)>1$ for all $t$. Thus
at each time, for $x$ near the minimum point one has
$1=h(x,t)^0>h(x,t)^1>h(x,t)^2>h(x,t)^3$, suggesting that the larger
$n=0,1,2,3$ is, the slower the fourth-order diffusion will be (near the minimum). Similarly, since the maximum is larger than $1$ we have
$1=h(x,t)^0 n_0$ then solutions are
positive for all time while solutions can touch down in finite time if
$n < n_0$. Here, we seek the analogous critical exponent $n_0(q)$.
The mobility coefficients in $h_t = -(h^n h_{xxx})_x - \mathcal{B} (h^m h_x)_x$
can affect whether a positive solution can become zero somewhere in
finite time. For example, if $7/2 < n \leq m < n+2$ then it has been
proved that the solution stays positive for all $t > 0$, by
\cite[\S4.2]{BPLW}. Note that $m3$,
Bertozzi and Pugh \cite{BPLW} conjecture that solutions could blow up,
with $\|h(\cdot,t)\|_{H^1} \to \infty$ in finite time. However, the
methods of \cite{BPFS} do apply to show that that if $n>7/2$ then the
solution remains positive for as long as it exists. Thus the critical
exponent for touch--down singularities satisfies $n_0(q) \leq 7/2$.
For $q=1$ and $n=1$, Goldstein et al.\ \cite[\S4]{GPS97} presented
simulations
%with $\mathcal{B}>1$
suggesting a finite-time singularity is possible. Bertozzi and Pugh
\cite{BPLW} presented numerical simulations for $q=1$ and $n=3$ in
which the solutions remain positive for all time, though they appear
to converge to droplets as $t \to \infty$. So it seems that $1 <
n_0(1) < 3$.
For $q=5/2$, we find that when $n=1$ the solution appears to touch
down in finite time, implying $1 \leq n_0(5/2) \leq 7/2$. To further
approximate $n_0(5/2)$, we performed simulations with a variety of
$n$-values, all with the initial data $h_{\rm ss} - .001 h_{\rm ss}''$. Our
findings suggest
\[
1.65 < n_0(5/2) < 1.6625.
\]
\begin{figure}[ht]
\begin{center}
\includegraphics[width=3.25in]{q25_find_nc.ps}
%/data/users/mpugh/Rick/Programs/Evolution/NewCode/Q2.5/HxxPerts/Mobilities
%
% used file `make_hmin_fig.m' to create figure
% n=1: started with 256 meshpoints, stopped with 8192, t= 43.60597999999820
% n=1.25: started with 256, stopped with 8192, t= 53.94800000000270
% n=1.5: started with 256, stopped with 8192, t= 71.62500000009160
% n=1.55: started with 256, stopped with 8192, t= 76.64000000020459
% n=1.6: started with 256, stopped with 8192, t= 82.49000000026810
% n=1.65: started with 256, stopped with 8192, t= 89.84000000040000
% n=1.6625: started with 256, stopped with 16384, t= 180.1800000000000
% refined to 32,768 and 65,536
% n=1.675: started with 256, stopped with 8192, t= 179.8500000000000
% n=1.7: started with 256, stopped with 8192, t= 180.5800000000000
% n=1.75: started with 256, stopped with 8192, t= 5182.700000069730
% n=2: started with 256, stopped with 8192, t= 55941.09999944990
%
% \hspace{1.5cm}
% \includegraphics[width=.35\textwidth]{q5_find_nc.ps}
%/data/users/mpugh/Rick/Programs/Evolution/NewCode/Q.5/HxxPerts/Mobilities
%
% used file `make_hmin_fig.m' to create figure
% n = 1, started with 256, stopped with 8192, t = 32.68134999999560
% n = 1.5, started with 256, stopped with 8192, t = 28.00800000000240
% n = 1.8, started with 256, stopped with 8192, t = 37.24000000012180
% n = 1.85, started with 256, stopped with 8192, t = 154.5800000000630
% n = 1.9, started with 256, stopped with 8192, t = 176.5400000001720
% n = 1.95, started with 256, stopped with 8192, t = 166.9900000000000
% n = 2, started with 256, stopped with 8192, t = 275.3600000000000
\end{center}
\caption{\label{find_nc} We fix the initial data and compute the
solution of (\ref{evolve}) where $q=5/2$ and $n=1$, $1.25$, $1.5$,
$1.55$, $1.6$, $1.65$, $1.6625$, $1.675$, $1.7$, $1.75$,and $2$. We
plot $\log_{10}(h_{min}(t))$ versus $t$. (In the plot, the graphs
move rightwards as $n$ increases.) For the first five values of $n$
the graph appears to go to $-\infty$ in finite time.}
\end{figure}
The evidence is presented in Figure~\ref{find_nc}, where we plot
$\log_{10} h_{min}(t)$ versus $t$. If $h_{min}(t)$
is decreasing at an exponential rate then the graph will be linear.
If $h_{min}(t)$ is decreasing to zero in finite time with an algebraic
rate then the graph will drop to $-\infty$ at that time with a
vertical slope. From Figure~\ref{find_nc}, if $n=2$ then $h_{min}$
decreases monotonically in time, eventually decreasing with an
exponential rate. For $n=1.75$, $1.7$, $1.675$, and $1.6625$,
$h_{min}$ decreases and then increases, ultimately decreasing with an
exponential rate. (We ran the simulations many decades beyond those
shown, to verify the exponential rate of decrease.) The solutions
with $n < 1.6625$ appear to be touching down in finite time. However
the $n=1.6625$ simulation gives a note of caution; it is possible that
the simulations with $n < 1.6625$ would run until $h_{min}$ became
quite small but would then increase somewhat and ultimately decrease
exponentially. We stopped our simulations when the solutions were no
longer numerically resolved (see \S\ref{stop}--\ref{point_double}).
To check that the bound on $n_0(5/2)$ is not dependent on our
choice of initial data, we chose ten large random perturbations with
$||v||_\infty = .178 \sim {h_{\rm ss}}_{min}$. Seven of the resulting solutions
relaxed to the constant steady state and three appeared to touch down
in finite time. We then fixed those initial data and varied the value
of $n$. For all three initial data, we found the
above upper and lower bounds
on $n_0(5/2)$.
The critical exponent $n_0(q)$ certainly depends on the value of $q$.
For example, for $q=1/2$, we observe similar phenomena to the $q=5/2$
case presented in Figure \ref{find_nc}, and conclude
\[
1.8 \leq n_0(1/2) < 1.85.
\]
To sum up, we find that changing the mobility coefficients ($m$, $n$)
while keeping the energy landscape fixed can change the nature of
trajectories across the energy landscape: if $n > n_0(q)$ then the
solution is smooth and classical at all times, while if $n < n_0(q)$
then there can be times when the solution is weak.
\subsection{Splitting singularities}
\label{one_two_sec}
% Our work in \S\ref{pert_hets} suggests that heteroclinic orbits can be
% preserved under some changes of the mobility. On the other hand,
% qualitative features of a solution can change significantly when the
% mobility is changed. For example, in \S\ref{time_pinch} we
% demonstrated that the mobility can affect the regularity of the
% solutions; for sufficiently large $n$, solutions are classical for all
% time while for smaller $n$, solutions can become nonnegative weak
% solutions in finite time. In this section,
%
Here we demonstrate another effect of changing the mobility: a
solution that evolves towards touch--down at just one local minimum
(per period) can change into a solution with two local minima (per
period).
\subsubsection*{$\mathbf{q=5/2}$}
\label{one_two_2.5}
If $n > 0$ and $q = 5/2$ then positive periodic smooth solutions of
(\ref{evolve}) remain bounded in $H^1$ for as long as they exist,
since $\|h(\cdot,t)\|_{H^1} \leq M < \infty$ by \cite{BPLW}. Because
the energy $\mathcal{E}$ decreases in time, we expect that these solutions will
converge to a steady state, as $t \to \infty$. The positive periodic
steady state is linearly unstable, and so we expect solutions to
converge either to the constant steady state or to a configuration of
steady droplets.
We take the same initial data as in \S\ref{time_pinch}. There, we
found that if $n < 1.6625 \approx n_0(5/2)$ then solutions appear to
touch down in finite time, and if $n > n_0(5/2)$ then touch--down is
in infinite time. Here we investigate the nature of the touch--down
more closely.
We find for $n = 1$ that the solution touches down in finite time at
one point per period, consistent with a long--time limit of one
droplet per period (left side of Figure~\ref{mobility_q2.5_n1}). For
$n=2$ the solution appears to be positive at all times and to touch
down at two points per period in the long--time limit (right side of
Figure~\ref{mobility_q2.5_n1}). This is consistent with a long--time
limit of two steady droplets per period. But it is impossible to
contain two zero contact angle steady droplets in an interval of
length $1$ (see \S\ref{mobility_q2.5_params}). In fact, we find that
the small `proto-droplet' flanked by the local minima is actually
draining, with its maximum decreasing to zero like $t^{-.4024}$.
\begin{figure}[ht]
\begin{center}
\includegraphics[width=.35\textwidth]{q25_n1.ps}
% /data/users/mpugh/Rick/Programs/Evolution/NewCode/Q2.5/HxxPerts/Mobilities/Expon1
\hspace{1.5cm}
\includegraphics[width=.35\textwidth]{q25_n2.ps}
% /data/users/mpugh/Rick/Programs/Evolution/NewCode/Q2.5/HxxPerts/Mobilities/Expon2
\end{center}
\caption{\label{mobility_q2.5_n1}$q=5/2$. Dashed line: initial data.
The same initial data is used for both simulations. Solid lines: the
solution at various times. Left: $n=1$; $h_{min}(t)$ occurs at $x=0$.
As $t$ increases, $h_{min}(t)$ decreases and $h_{max}(t)$ increases.
At all times there is one minimum per period. Right:
$n=2$. Again, as $t$ increases $h_{min}(t)$ and $h_{max}(t)$ decrease
and increase respectively. Initially there is only one local minimum
but after some time it splits into two. The minima flank a small
`proto--droplet' and suggest a possible long--time limit of two
droplets per period. But, in fact, the proto--droplet drains away as
$h_{min} \to 0$.}
%%/hans/data/users/mpugh/Rick/Programs/Evolution/NewCode/Q2.5/HxxPerts/Mobilities/Expon1
%/hans/data/users/mpugh/Rick/Programs/Evolution/NewCode/Q2.5/HxxPerts/Mobilities/Expon2
\end{figure}
Our simulations suggest that a second critical exponent, which we call
$n_1(q)$, governs the number of local minima per period. If
$n < n_1(q)$ then there is one minimum per period. If $n >
n_1(q)$ then there are two local minima per period, with their
positions moving in time. That is, the single local minimum splits
into two as $n$ increases through $n_1(q)$. Goldstein et al.\
\cite[\S4C]{GPS97} observed something in a similar spirit for $h_t = -
(h h_{xxx})_x - \mathcal{B} (h h_x)_x$, namely a single symmetric singularity
that splits into a pair of asymmetric singularities as $\mathcal{B}$ increases
past $\mathcal{B} \approx 1.35$.
\begin{figure}[ht]
\begin{center}
\includegraphics[width=3.25in]{q25_one_two_pinch.ps}
%/data/users/mpugh/Rick/Programs/Evolution/NewCode/Q2.5/HxxPerts/Mobilities
% n = 1: 256 meshpoints, stopped with 8192, t= 43.60597999999820
% n = 1.25: 256, stopped with 8192, t= 53.94800000000270
% n = 1.5: 256, stopped with 8192, t= 71.62500000009160
% n = 1.6: 256, stopped with 8192, t= 82.49000000026810
% n = 1.7: 256, stopped with 8192, t= 180.5800000000000
% \hspace{1.5cm}
% \includegraphics[width=.35\textwidth]{q5_one_two_pinch.ps}
%/data/users/mpugh/Rick/Programs/Evolution/NewCode/Q.5/HxxPerts/Mobilities
% n = 1: 256 to 8192, stops at t=32.68134999999560
% n = 1.25: 256 to 8192, stops at t=24.60470000000170
% n = 1.4: 256 to 8192, stops at t=26.42700000000160
% n = 1.6: 256 to 8192, stops at t=30.08750000000080
% n = 1.8: 256 to 8192, t = 37.24000000012180
% n = 2: 256 to 8192, t = 275.3600000000000
\end{center}
\caption{\label{one_two} $q=5/2$. Fixed initial data, late--time
profiles for the solution of (\ref{evolve}) computed with different
values of $n$. Top: $n=1$ (single minimum) and $1.25$ (two minima). Bottom: $n =
1.5$, $1.6$, $1.7$, where the distance between the two minima increases with
$n$.}
\end{figure}
We find that
\[
1 < n_1(5/2) < 1.25 .
\]
In Figure~\ref{one_two} we plot late--time profiles for a range of
$n$. In the top plot, we plot the solutions near $x=0$
for $n=1$ and $1.25$. The $n=1$ profile has only one local minimum,
while the $n=1.25$ profile has two. In the bottom plot, we plot the
solutions near $x=0$ for $n=1.5$, $1.6$ and $n=1.7$. Each profile has
two local minima, with the distance between the minima increasing with
$n$. For each $n \geq 1.25$, we find that the proto-droplet is
draining with its maximum decreasing to zero with a rate that depends
on $n$.
% The $n=3$ and $4$ evolutions appear to be very similar to the $n=2$
% evolution. We did not do any simulations beyond $n=4$ since the
% larger the value of $n$, the longer the simulation had to run before
% we could observe anything tangible --- the mobility coefficient $h^n$
% is very small near the local minimum, where $h$ is small. We could
% not compensate by taking large time--steps, because $h^n$ could be
% quite large near the local maximum, since $h_{max}(t) > 1$. (Large
% differences in time--scales are difficult to handle numerically.)
The phrase `splitting singularities' is perhaps a bit of a misnomer.
The two local minima described above do not appear to touch down in a
way that yields isolated singularities. Rather, the solution between
them appears to be touching down --- they are the endpoints of a
developing dry interval. A similar phenomenon was observed by
Constantin et al.\ \cite[\S{III},IV]{First} with $n=1$ and $\mathcal{B}=0$ (and
with different boundary conditions): their proto-droplet decayed like
$1/t$.
\subsubsection*{$\mathbf{q=1/2}$}
\label{one_two_.5}
%$\mathcal{B} = 36.29$
For this value of $q$, we took initial data $h_{\rm ss} - .001 h_{\rm ss}''$ and observed
phenomena similar to those in the $q=5/2$ case, finding
$ 1 < n_1(1/2)<1.25$.
% see the plots to the right of Figure~\ref{one_two}. In the right top
% plot, we present the final resolved solutions for $n=1$ and $1.25$
% with $32,\!768$ meshpoints. The local minima for $n=1.25$ are much
% closer to each other than in the $q=5/2$ plot. In fact, at $8192$
% meshpoints it appeared that there would be only one local minimum; as
% more time passed it split into two. In the right bottom plot of
% Figure~\ref{one_two} we plot solutions with $n=1.4$, $1.6$, $1.8$, and
% $2$. The $n=1.4$, $1.6$, and $1.8$ solutions are with $8192$
% meshpoints, and the $n=2$ solution is with $16,\!384$ meshpoints. The
% profiles are all resolved and were chosen to have comparable
% $h_{min}$. The distance between the local minima in the plot increase
% monotonically with $n$.
%
We observe behavior similar to that shown in
Figure~\ref{mobility_q2.5_n1} for $q = 5/2$. Specifically, the
long--time limit appears to be one droplet. This is interesting
since, unlike for $q = 5/2$, if $q=1/2$ then it is possible to have
two disjoint steady droplets in an interval of length $1$ (see
\S\ref{mobility_q.5_params}). And so the proto-droplet could, in
theory, converge to a steady droplet. Nonetheless, like for $q=5/2$,
the proto-droplet appears to drain away.
We have not been able to find examples of $q$, $n$, and initial data
that yield a solution whose long--time limit is a configuration of more than one steady droplet per period.
\section{Dynamics: $q$ and its effect on heteroclinic orbits}
\label{simulations}
%
%
We now consider the evolution equation with a variety of $q$-values,
taking a wide range of initial data near steady states. The resulting
solutions display a diversity of behaviors. Our
stability theorems \cite{LP2,LP3} often allow us to predict the
observed short--time behavior, and our theorems on the energy levels
of steady states \cite{LP3} often allow us to guess the long--time
limit of the solution. As part of this work we predict (and find
strong numerical evidence for) heteroclinic connections between
certain steady states. Interestingly, there are several cases where we \emph{cannot} predict the long-time behavior of the solution; in particular see the cases $q=1$ and $3/2$ below.
% Recall the steady states depend on the parameter
% \[
% q = m-n+1.
% \]
We present results for seven values of $q$:
\[
q = -3, \; 1/2, \; 1, \; 3/2, \; 1.768, \; 5/2, \; 4,
\]
chosen from the intervals
$\{(-\infty,-1],(-1,1),(1,7/4],(7/4,1.794),(1.795,3),[3,\infty) \}$ in
which our theorems from \cite{LP1,LP2,LP3} suggest the solutions will
display distinct behaviors. (How the above intervals were chosen will
become clear in what follows. Also, we did study other values of $q$ and
found that the phenomena reported here are robust.) For $q=-3$ we take $n=3$ and $m=-1$. For the other six $q$-values we take $n=1$ and $m=q$. Our
numerical simulations are not greatly affected if we change $n$ and
$m$ in a manner that keeps $q$ fixed, except for the features reported in Section~\ref{mobility}.
%For the $q=1$ case, numerical simulations were earlier presented with
%$m=n=1$ in \cite{GPS97}, with $m=n=2$ in \cite[\S8]{Gruen00}, and with
%$m=n=3$ in \cite{BPLW}.
\subsection{$\mathbf{q=-3}$: the van der Waals case}
\label{qm3}
{\it Characteristic features for $q \in (-\infty,-1]$: positive
periodic steady states are linearly unstable and there are no droplet
steady states with acute contact angles.} (See bifurcation
diagram~\ref{bifurcs}a and \cite[\S2.2]{LP1}.) We study the equation
%
\begin{equation} \label{vanderwaal}
h_t = - (h^3 h_{xxx})_x - \mathcal{B} (h^{-1}h_x)_x ,
\end{equation}
%
for which $q=-3$. The equation was proposed by Williams and Davis
\cite{WD82} to model a thin liquid film with net repulsive van der
Waals interactions, and more recently it has been studied by Zhang and
Lister \cite{ZL99} and by Witelski and Bernoff \cite{WB99,WB00}.
\subsubsection{$q=-3$. Perturbing the positive periodic
steady state}
\label{qm3pp}
First, a positive periodic steady state for $\mathcal{B} = .08930$ is
constructed. It is linearly unstable (see Figure \ref{bifurcs}a).
There is at least one linearly unstable eigenfunction for this steady
state \cite{LP2} and hence the unstable manifold is at least
one-dimensional. The weakly nonlinear stability analysis suggests
that, at least for nearly-constant positive periodic steady states,
the unstable manifold is exactly one dimensional.
\noindent $\bullet$ {\it Even perturbations.} The steady state is
even and the evolution equation preserves this, after an even perturbation. First we perturb $h_{\rm ss}$
with the even, zero-mean perturbation $\pm \epsilon h_{\rm ss}''$. Since the
perturbation $+ \epsilon h_{\rm ss}''$ lowers the maximum and raises the minimum,
one might hope the resulting solution would converge to the constant
steady state $h \equiv 1$. If this happens for all small $\epsilon$, then
this would be strong evidence for existence of a heteroclinic orbit
connecting $h_{\rm ss}$ to the constant steady state.
\begin{figure}[ht]
\begin{center}
\includegraphics[width=3.25in]{qm3_relaxing.ps}
%/data/users/mpugh/Rick/Programs/Evolution/NewCode/Qm3/HxxPerts/Pos
% done with adaptive timestepping in
%/data/users/mpugh/Rick/Programs/Evolution/NewCode/Qm3/HxxPerts/Pos/Adaptive
% verified that h_min decreases a little before increasing.
\end{center}
\caption{\label{qm3_relax} $q=-3,n=3$. Dashed: initial data $h_{\rm ss}
+.0001 h_{\rm ss}''$. Solid: the solution at a number of times. The local
extrema are fixed in space and, after a short transient, $h_{min}$
and $h_{max}$ converge monotonically to $1$ and the solution $h$
relaxes to the constant steady state.}
\end{figure}
%
Perturbing in the direction $+ .0001 h_{\rm ss}''$, we find that after a
short transient, the local minimum increases (and maximum decreases)
to $1$. The extremal points remain fixed in space, while the solution
relaxes to the constant steady state as $t \to \infty$; see
Figure~\ref{qm3_relax}. (This was
observed previously in \cite[Figure 4b]{WB00}.) We repeated the
simulation for smaller values of $\epsilon$ and found that all the
perturbations yielded solutions that relaxed to the constant steady
state as $t \to \infty$. This is convincing evidence that
heteroclinic orbits connecting $h_{\rm ss}$ to the constant steady state
exist. There are also theoretical reasons to suspect these
heteroclinic orbits exist: (i) $h_{\rm ss}$ is energy unstable in the
directions $\pm h_{\rm ss}''$ by \cite[Theorem~2]{LP3}, (ii) the energy of
$h_{\rm ss}$ is higher than that of the constant steady state $h \equiv 1$ by
\cite[Theorem~6]{LP3} (also observed numerically by Witelski and
Bernoff \cite[\S3]{WB00}), and (iii) the constant steady state is a
local minimum of the energy $\mathcal{E}$ by \cite[Theorem~10]{LP3}.
Next we perturb in the opposite direction with $- .0001 h_{\rm ss}''$, so
that the maximum is raised and the minimum lowered. Since the
perturbation decreases the energy $\mathcal{E}$, we might expect the solution to
subsequently converge to a droplet steady state or to a configuration
of droplet steady states. If such a droplet exists it must have
$90^\circ$ contact angles, by \cite[\S2.2]{LP1}. We find that after a
short transient, the minimum height of the solution decreases in time,
appearing to decrease to zero in finite time. The top left plot in
Figure~\ref{n32768_evolve} presents the evolution of the solution near
$x=0$. The local extrema are fixed in space, with the solution
appearing to touch down at one point per period. (This was shown
previously in \cite[Figure 4c]{WB00}.) Computing the derivative $h_x$
of the solution, we find that its maximum and minimum values grow as
time passes, as in the bottom left plot of Figure~\ref{n32768_evolve}.
These extremum points of $h_x$ move in time, moving toward $x=0$ as
the singular time approaches. This is consistent with a solution that
touches down with $90^\circ$ contact angles in finite time. The right
plots in Figure~\ref{n32768_evolve} show a late--time profile of $h$.
%
% But note that we did not design the code to study the formation of
% finite--time singularities; the solution has decreased by only a
% factor of $6.04$ due to limitations of our uniform--mesh code (see
% \S\ref{stop}).
\begin{figure}[ht]
\begin{center}
% /data/users/mpugh/Rick/Programs/Evolution/NewCode/Qm3/HxxPerts/Neg
% \includegraphics[width=.35\textwidth]{alt_n131072_evolve.ps}
% /data/users/mpugh/Rick/Programs/Evolution/NewCode/Qm3/HxxPerts/Neg/Adaptive
\includegraphics[width=.35\textwidth]{new_n32768_evolve.ps}
% done with adaptive timestepping in
% /data/users/mpugh/Rick/Programs/Evolution/NewCode/Qm3/HxxPerts/Neg/Adaptive
% verified that h_min increases a little before decreasing
\hspace{1cm}
% /data/users/mpugh/Rick/Programs/Evolution/NewCode/Qm3/HxxPerts/Neg
% \includegraphics[width=.35\textwidth]{final_h.ps}
% /data/users/mpugh/Rick/Programs/Evolution/NewCode/Qm3/HxxPerts/Neg/Adaptive
\includegraphics[width=.35\textwidth]{new_final_h.ps}
\end{center}
\caption{\label{n32768_evolve}Left top: dashed line, initial data $h_{\rm ss}
- .0001 h_{\rm ss}''$; solid lines, solution $h$ at later times. The local
minimum is fixed in space and, after a short transient, decrease
monotonically to zero. Left bottom: dashed line, initial slope; solid
lines, $h_x$ at same times as above. As $t$ passes, $||h_x||_\infty$
increases and the positions of the extrema move in toward $x=0$. This
suggests $90^\circ$ contact angles are forming as the singular time
approaches. Right top: solution at a late time. Right bottom: a
close--up near $x=0$ at the same time.}
\end{figure}
The work of Zhang and Lister \cite[\S5]{ZL99} on similarity solutions
suggests that $h(x,t) \sim (t_c - t)^{1/5} H(x/(t_c - t)^{2/5})$ as
touchdown approaches; here $t_c$ is the time of touchdown and $H$ is a
particular positive function with $H(\eta) \sim (0.807) \; \mathcal{B}^{1/4}
|\eta|^{1/2}$ for large $\eta$.
%
%The rescaled function $\mathcal{B}^{-1/4} h(x,\mathcal{B}^{-3/4} t)$
%satisfies (\ref{vanderwaal}) with $\mathcal{B}=1$, as considered
%by Zhang and Lister.
%
Our computations are consistent with this ansatz. See
\cite{WB99,WB00} for more on the similarity solutions of
(\ref{vanderwaal}).
%
%We found the outer time--scale of $h_{min} \sim (t_c-t)^{1/5}$
%(really $.2000$) and we found the inner time--scale (if we
%assume the outer is exact) of $\eta \sim x/(t_c-t)^{2/5}$
%(really $-0.5902$) by fitting ${h_{xx}}_{max}$.
%But the best we could do with the spatial structure
%was $H(\eta) \sim \eta^{.5207}$\dots}
After the singular time, one possible behavior of the solution shown
in Figure~\ref{n32768_evolve} is that the solution becomes a
nonnegative weak solution. Then it might relax, as a weak solution, to
a droplet steady state with $90^\circ$ contact angles. Alternatively,
the solution might, at some later time, become positive and classical
again, ultimately relaxing to the constant steady state. This is
certainly possible since each profile shown in
Figure~\ref{n32768_evolve} has higher energy $\mathcal{E}$ than the constant
steady state. We note that pursuit of this question will require
further analysis, because the current weak existence theory for the
evolution equation requires $ n \leq m$, whereas here we have $m < n$.
Further, the current weak solution theory does not admit $90^\circ$
contact angles unless they occur for a set of times of zero
measure. It may be that completely new techniques will have to be
developed, to study this equation.
\noindent $\bullet$ {\it Other perturbations.}
To check the degree to which the behaviors described above depend on
the choice of perturbation, we performed a number of runs with random
perturbations. (See \S\ref{rand_perts}.) We found that all the
solutions either relaxed to the constant steady state or else appeared
to touch down in finite time. (The gross dynamics are as in Figures
\ref{qm3_relax} and \ref{n32768_evolve}; the finer dynamics concern
the positions of the local extrema as a function of time.)
\subsubsection{$q=-3$. Perturbing the constant steady state}
%We continue to take $q=-3, n=3, m=-1$.
Suppose $\mathcal{B} < 4 \pi^2$, so that the constant steady state $h \equiv 1$
is linearly stable with respect to zero--mean perturbations of period
$1$. By \cite[Theorem~10]{LP3}, the constant steady state is a strict
local minimum of the energy, and is dynamically stable.
This is the uninteresting case.
\begin{figure}[ht]
\begin{center}
\includegraphics[width=3.25in]{qm3_pert_from_mean.ps}
% /data/users/mpugh/Rick/Programs/Evolution/NewCode/Qm3/PerturbingTheConstant
% adaptive timestepping in
% /data/users/mpugh/Rick/Programs/Evolution/NewCode/Qm3/PerturbingTheConstant/Adaptive
\end{center}
\caption{\label{qm3_meanpert}$q=-3,n=3$. Dashed: initial data $1 -
.0002 \cos(2 \pi x)$. Solid: solution at various times. Top: the
short--time dynamics. The local minimum decreases, while the local
maximum increases for a while. The top then `flattens' and two local
maxima form one to each side of the flat region. Bottom: later-time
dynamics. The solution appears to touch down at one point per period
and continues to have two local maxima per period. The local minimum
is at $x=0$ and, after a short transient, decreases monotonically to
zero.}
\end{figure}
Now suppose $\mathcal{B} > 4 \pi^2$. Then the constant steady state $h \equiv 1$
is a saddle point for the
diagram \ref{bifurcs}a), we suspect that a perturbation of the
constant steady state will yield a solution that touches down in
finite or infinite time.
To investigate, we first take $\mathcal{B} = 2.467 < 4 \pi^2 $. For
all initial data we considered, we found that the resulting solutions
appear to relax to the constant steady state. This is as predicted.
We then took $\mathcal{B} = 631.7 > 4 \pi^2$. Here, the constant steady state
has $39$ linearly unstable eigenmodes and is a saddle point of the
energy. For initial data $1 - .0002 \cos (2 \pi x)$, the solution
appears to touch down in finite time, with one touchdown per period.
Figure~\ref{qm3_meanpert} shows this evolution over two periods.
\subsection*{Remarks}
In our studies of positive periodic steady states and constant steady
states, we checked that the observed phenomena persist with smaller
perturbations: the behaviors are robust. For this reason, in the
remainder of the article we will not discuss smaller perturbations.
Further, we will not discuss random perturbations or odd
perturbations, since we found that the observed phenomena were like
those observed for even perturbations. (Except that for even
perturbations the local extrema are fixed in space while for other
perturbations they move slightly as the solution evolves.)
We also will not present any further discussions of perturbations of
the constant steady state. We found that the behaviors were always
those predicted by the bifurcation diagram --- if the constant steady
state is linearly stable then small perturbations converge back to it,
while if the constant steady state is linearly unstable then
perturbing it yields a solution that converges to a stable positive periodic steady state, if one exists. If none exists, then we found that the solution either touches down in finite or infinite time ($q<3$) or
else it blows up in finite time ($q \geq 3$).
% The final profile presented does not look like any known steady state,
% and so we expect the solution to continue evolving after touching
% down.
%
%is clearly not near a constant steady steady state because the
%constant steady state satisfies a nonlinear
%oscillator formulation, which precludes a positive local minimum
%unless it is part of a positive periodic solution.
% \noindent $\bullet$ {\it Random perturbations.}
% Using random perturbations, we verified that the above results are
% robust: the steady state $\overline{h} = 2$ is asymptotically stable and
% $\overline{h} = 1/2$ is unstable.
%
% \vskip 6pt
% Incidentally, we verified that the evolution is exponential in time
% near the periodic and constant steady states. This is consistent with
% the nonlinear behavior being dominated by the linear theory when the
% solution is sufficiently near a steady state. We found there was a
% short transient before the exponential behavior began, suggesting that
% the direction $h_{\rm ss}''$ is near but not equal to the first
% eigendirection. During the transient time, the solution is locating
% this eigendirection.
\subsection{$\mathbf{q=1/2}$}\
\label{q.5}
{\it Characteristic features of $q \in (-1,1)$: positive periodic
steady states are linearly unstable. A `Mountain pass' scenario can
occur --- the energy of the non-constant positive periodic steady
state is higher than the energies of the constant steady state and the
zero-contact angle droplet steady state.} (See bifurcation
diagram~\ref{bifurcs}b and remarks after \cite[Theorem~11]{LP3}.)
\vskip 3pt
We compute solutions of
\[
h_t = - (h^1 h_{xxx})_x - \mathcal{B} (h^{1/2} h_x)_x.
\]
%
%
First, a positive periodic steady state for $\mathcal{B} = 35.59$ is
constructed. It is linearly unstable (see Figure \ref{bifurcs}b). As
in \S\ref{qm3pp}, we perturb $h_{\rm ss}$ with $\pm .0001 h_{\rm ss}''$.
For the initial data $h_{\rm ss} + .0001 h_{\rm ss}''$, we see a heteroclinic connection
to the constant steady state, very much like in the $q=-3$ case shown in
Figure~\ref{qm3_relax}.
For the initial data $h_{\rm ss} - .0001 h_{\rm ss}''$ we find the solution appears
to touch down in finite time. Like the $q=-3$ simulation in Figure~\ref{n32768_evolve},
$h_{min}(t)$ is located at $x=0$, and after a short transient the
minimum decreases monotonically in time. But unlike the $q=-3$
simulation, the profile seems to touch down with {\em zero} contact
angles. There does exist a zero-angle droplet steady state $\hat{h}_{\rm ss}$
that has the same area as $h_{\rm ss}$, has length less than $1$, and has
lower energy than $h_{\rm ss}$, by \cite[Theorem~7]{LP3}. Presumably this droplet steady state is the intended long-time limit of the solution, up to translation. But this cannot currently be proved, because the known zero-angle weak existence theory requires $0 -1$.
\subsection{$\mathbf{q=1}$}\
\label{q1}
{\it Characteristic features for $q =1$: all positive periodic steady
states are linearly neutrally stable.} (See bifurcation
diagram~\ref{bifurcs}c and \cite[Lemma~4]{LP3}.)
\vskip 3pt
The non-constant positive periodic steady states are neutrally
stable, when $q=1$. We take $n=m=1$ and compute solutions of
\[
h_t = -(h^1 h_{xxx})_x - \mathcal{B} (h^1 h_x)_x.
\]
Numerical simulations for $q=1$ have been presented before by other
authors: with $m=n=1$ in \cite{GPS97}, with $m=n=2$ in
\cite[\S8]{Gruen00}, and with $m=n=3$ in \cite{BPLW}. But the latter
two articles do not consider Bond numbers for which periodic steady
states might be observed. In the first article, Goldstein et al.\
\cite[Fig. 3a]{GPS97} found that fairly large multi-modal
perturbations of positive periodic steady states yield relaxation to
(generally different) steady states.
First, we constructed a positive periodic steady state for $\mathcal{B} =
39.48$. In the left plot of Figure~\ref{q1_fig}, we present two
simulations confirming that this steady state is dynamically stable,
with a small perturbation yielding convergence to a nearby positive
periodic steady state. In both cases, the solution relaxes to a
positive periodic steady state with a local minimum close to $x=0$ and
an amplitude close to $.8$. We have no rule for predicting the
amplitude of the long--time limit and, unless the perturbation is
even, we have no way of predicting the position of the local minimum.
The long--time dynamics will be especially difficult to predict when
$q=1$ since there are infinitely many $1$-periodic steady states all
having area $1$. (In the $q \neq 1$ case there are at most two such
steady states.)
\begin{figure}[ht]
\begin{center}
\includegraphics[width=.35\textwidth]{q1_rand_perts.ps}
%/data/users/mpugh/Rick/Programs/Evolution/NewCode/Q1/EvenPerts/
\hspace{1.5cm}
\includegraphics[width=.35\textwidth]{q1_flat_bottom.ps}
%/data/users/mpugh/Rick/Programs/Evolution/NewCode/Q1/EvenPerts/ShorterPerts
\end{center}
\caption{\label{q1_fig}$q = 1$, $n=1$; dashed line: the initial data.
Left: the solid line is a late--time profile of the solution. We
found that the solution relaxes to a steady state close to the
original steady state $1 - .8 \cos(2 \pi x)$. Top left: $h_0(x) = 1 -
.8 \cos(2 \pi x) + .3 v(x)$ where $v$ is a random zero--mean
perturbation. Bottom left: $h_0(x) = 1 - .8 \cos(2 \pi x) - .19
\exp(-100 \sin^2(\pi x)) + .19 \exp(-100 \sin^2(\pi (x-1/2)))$.
Right: $h_0(x) = 1 - 1.143 \cos(2 \pi x) + .2714 \cos(4 \pi x)$.
Solid lines are $h$ at various times. The extrema are fixed in space
and, after a short transient, decrease/increase monotonically with
$h_{min}(t)$ touching down in finite time.}
\end{figure}
Since these simulations suggest that the positive periodic steady
states are dynamically stable, one might guess that solutions cannot touch
down in finite time. (This is what we observe later for $q=3/2$ and
$q=1.768$.) And as the bottom left plot
of Figure~\ref{q1_fig} suggests, initial data that has a {\em sharp}
local minimum will likely not evolve towards touch--down; the local
minimum retracts in time, as expected for a solution of a
surface--tension driven flow. On the other hand, initial data that is
very flat near its local minimum does appear to lead to touch--down in
finite time, as shown in the top right plot of Figure~\ref{q1_fig}.
The bottom right plot shows this evolution near the touch--down point.
%
\subsection{$\mathbf{q=3/2}$}\
\label{q1.5}
{\it Characteristic features for $q \in (1,7/4]$: positive periodic
steady states are linearly stable.} (See bifurcation diagrams
\ref{bifurcs}d--e.)
\vskip 3pt
We compute solutions of
$$
h_t = - (h^1 h_{xxx})_x - \mathcal{B} (h^{3/2} h_x)_x.
$$
We constructed a positive periodic steady state $h_{\rm ss}$ for $\mathcal{B} = 40.07$.
It is linearly stable (see Figure \ref{bifurcs}d).
% data is in:
%
%/hans/data/users/mpugh/Rick/Programs/Evolution/NewCode/Q1.5/HxxPerts/Expon1/Pos
%
%/hans/data/users/mpugh/Rick/Programs/Evolution/NewCode/Q1.5/HxxPerts/Expon1/Neg
%
\begin{figure}[ht]
\includegraphics[width=3.25in]{q15_long_pert.ps}
% /data/users/mpugh/Rick/Programs/Evolution/NewCode/Q1.5/HxxPerts/Expon1/LongerPert
\caption{\label{q1.5_oddpert} $q=3/2$, $n=1$. Top plot: dashed line is
initial data $h_{\rm ss} - .0001 \cos(\pi x)$; heavy solid line is final
resolved solution; light solid lines are $h$ at various times. After
a short transient, $h(1,t)$ increases monotonically. At late time, a
pair of local minima form to either side of $x=1$ and touch down in
finite time. Bottom plot: close--up of the final resolved solution.}
\end{figure}
Also, every perturbation of the same and shorter period increases the
energy, by \cite[Theorem~5]{LP3}, and so we expect to observe
relaxation back to a translate of $h_{\rm ss}$. This is precisely what our
simulations show. We have not been able to predict the amount of
translation that occurs, but there is some hope of progress here,
since impressive results on a similar translation problem have been
obtained in \cite{bricmont,carlen} for the Cahn--Hilliard equation
$h_t = - h_{xxxx} - ((1 - 3 h^2) h_x)_x$ on the whole real line.
Next consider zero-mean perturbations of \emph{longer} period, to
which $h_{\rm ss}$ is linearly unstable by \cite[Theorem~1]{LP2}. We ask, to
what long-time behavior does this instability give rise? The
perturbation $-.0001 \cos(\pi x)$, for example, raises the local
minimum of $h_{\rm ss}$ at $x=1$ and lowers it at $x=0$. The top plot of
Figure~\ref{q1.5_oddpert} presents the resulting evolution of the
$2$-periodic solution. The solution appears to touch down in finite
time, though interestingly, it does not do so at $x=0$. Instead the
touchdown is driven by a dramatic increase in the solution near
$x=1$. The bottom plot of Figure~\ref{q1.5_oddpert} shows a close--up
of the final resolved solution. The smaller droplet, centered on
$x=0$, is not close to a steady droplet since it contains a local
minimum within itself --- an impossibility for a steady droplet.
Therefore we expect the solution would continue to evolve as a
nonnegative weak solution, relaxing either to a single steady droplet
or to some (unknown) configuration of steady droplets.
%
\subsection{$\mathbf{q=1.768}$}\
\label{q1.768}
{\it Characteristic features for $q \in (7/4,1.794)$: some positive
periodic steady states are linearly stable, while others
are linearly unstable; and there can be more than one positive
periodic steady state with the same period and area.} (See bifurcation
diagrams \ref{bifurcs}f--h, and \cite[\S5.1]{LP1}.)
\vskip 3pt
\begin{figure}[ht]
\includegraphics[width=3.25in]{q1768_heteroclinics.ps}
%/hans/data/users/mpugh/Rick/Programs/Evolution/NewCode/Q1.768/Case2/Expon1
\caption{\label{q1.768_het}$q = 1.768$ and $n=1$. The light solid
lines are $h$ at a sequence of times. The extrema are fixed in space.
(a) Dashed: initial data $h_{{\rm ss}2}+.0001h_{{\rm ss}2}''$; heavy
solid: constant steady state. Solution relaxes to constant. After a
short transient, $h_{max}$ and $h_{min}$ converge monotonically to $1$
and $h$ relaxes to the constant steady state. (b) Dashed: initial
data $h_{{\rm ss}2}-.0001h_{{\rm ss}2}''$; heavy solid: $h_{{\rm
ss}1}$. After a short transient, $h_{max}$ increases and $h_{min}$
decreases monotonically, with $h$ relaxing to $h_{{\rm ss}1}$.}
\end{figure}
We compute solutions of
\[
h_t = -(h^1 h_{xxx})_x - \mathcal{B} (h^{1.768} h_x)_x.
\]
Now the possibility arises of a heteroclinic connection between two
fundamentally different positive periodic steady states. For Bond
number $\mathcal{B} = 39.46$ we found two distinct
positive $1$-periodic steady states, $h_{{\rm ss}1}$ and $h_{{\rm
ss}2}$, that have area $1$. We denote the steady state that has
larger amplitude by $h_{{\rm ss}1}$, and the other by $h_{{\rm ss}2}$.
We expect $h_{{\rm ss}1}$ to be linearly stable and $h_{{\rm ss}2}$ to
be unstable, by \cite[Theorem~9]{LP3}, with $h_{{\rm ss}1}$ having
lower energy. That is, $h_{{\rm ss}1}$ lies on the stable branch of
the bifurcation diagram~\ref{bifurcs}g and $h_{{\rm ss}2}$ lies on the
unstable branch. The constant steady state is linearly stable because
$\mathcal{B} < 4\pi^2$.
%A_unst = 4.66055351849197, mean value = 0.74175013001232
%A_stab = 4.66055351848994, mean value = 0.74175013001200
Our simulations confirmed these predictions. First, all small
perturbations of $h_{{\rm ss}1}$ resulted in solutions that relaxed
back to $h_{{\rm ss}1}$. All perturbations of $h_{{\rm ss}2}$ yielded
solutions that either connect to the constant steady state or to
$h_{{\rm ss}1}$. Figure \ref{q1.768_het} shows a typical pair
of solutions.
\subsection{$\mathbf{q=5/2}$}
\label{q2.5}
{\it Characteristic features of $q \in (1.795,3)$: positive periodic
steady states are linearly unstable. `Mountain pass' scenario can
occur.} (See bifurcation diagram~\ref{bifurcs}i and remarks after
\cite[Theorem~11]{LP3}.)
\vskip 3pt
We compute solutions of
\[
h_t = - (h^1 h_{xxx})_x - \mathcal{B} (h^{5/2} h_x)_x .
\]
% We constructed a positive periodic steady state for $\mathcal{B} = 35.32$. It
% is linearly stable (see Figure~\ref{bifurcs}i). Our simulations
% confirmed our predictions: all perturbations lead to solutions that
% either relax to the constant steady state or that appear to touch down
% in finite time.
We ran one simulation, for $\mathcal{B} = 35.32$, and found behaviors that were
qualitatively the same as for $q= 1/2$, in \S\ref{q.5}.
\subsection{$\mathbf{q=4}$}\
\label{q4}
{\it Characteristic features of $q \in [3,\infty)$: positive periodic
steady states are linearly unstable, and if a positive periodic steady
state and a zero-angle droplet steady state have the same area, then
the period of the former is less than the length of the latter.
(See \cite[Theorem~7]{LP2} and the proof of \cite[Theorem~7]{LP3}.)}
\vskip 3pt
\begin{figure}[ht]
\begin{center}
\includegraphics[width=3.25in]{q4_evolve.ps}
% used adaptive time-stepping
%/data/users/mpugh/Rick/Programs/Evolution/NewCode/Q4/HxxPerts/Run0
\caption{\label{q4_blow}$q=4,n=1$. Dashed: initial data $h_{\rm ss} -.001
h_{\rm ss}''$. Solid: $h$ at a sequence of times. The extrema are fixed in
space and, after a short transient, $h_{max}$ increases monotonically
towards infinity. After a short transient, $h_{min}$ decreases
monotonically to a positive value.}
\end{center}
\end{figure}
We compute solutions of
%
\begin{equation} \label{blow}
h_t = - (h^1 h_{xxx})_x - \mathcal{B} (h^4 h_x)_x .
\end{equation}
%
This equation is `super-critical' in the sense of Bertozzi and Pugh
\cite{BPLW}, since $m>n+2$
({\it i.e.} $q>3$). According to their conjecture in \cite{BPLW},
positive periodic solutions can blow up in finite time, with
$\|h(\cdot,t)\|_\infty \to \infty$. Bertozzi and Pugh made the same
conjecture for compactly supported weak solutions on the line, and
proved in \cite{BPFS} that blow-up can occur in finite time when $n=1$
and $m \geq n+2 = 3$. Specifically, they proved for such cases that
if the compactly supported initial data $h_0$ has negative energy
$\mathcal{E}(h_0)<0$,
%\[
% \mathcal{E}(h_0) = \int \left[ \frac{1}{2} h_0'(x)^2
% - \frac{\mathcal{B}}{q(q+1)} h_0(x)^{q+1} \right] dx < 0,
%\]
then the compactly supported weak solution blows up in finite time
with its $L^\infty$ and $H^1$ norms going to infinity.
Here we present computational evidence that positive periodic
solutions of (\ref{blow}) can also blow up in finite time. Further,
we find initial data that has positive energy yet still appears to
yield finite--time blow-up, suggesting that $\mathcal{E}(h_0) < 0$ is not
necessary for blow-up, in the periodic case.
We took $\mathcal{B} = 22.60$, and found a linearly unstable positive periodic
steady state. We considered initial data $h_{\rm ss} \pm .001 h_{\rm ss}''$. The
initial data $h_{\rm ss} + .001 h_{\rm ss}''$ yielded a solution that relaxed to the
constant steady state.
The initial data $h_{\rm ss} - .001 h_{\rm ss}''$ yielded a solution that appears to
blow up in finite time (see Figure~\ref{q4_blow}). A self-similarity
ansatz suggests that $h(x,t) \sim (t_c - t)^{-1/7} H((x-1/2)/(t_c -
t)^{3/14})$ as blowup approaches. Here $t_c$ is the time of blowup and
$H$ is a positive function with $H(\eta) \sim C |\eta|^{-2/3}$ for
large $\eta$. Our simulations are consistent with this ansatz.
%
% Again, if we make the ansatz then we can estimate the blowup time
% $t_c$, since taking the ratio of the computed values of $h(1/2,t)$ at
% two late times $t_1$ and $t_2$ gives the value of $(t_c -
% t_1)^{-1/7}/(t_c - t_2)^{-1/7}$, from which $t_c$ can be determined.
% We find that $t_c$ is slightly larger than the final resolved time.
%
Self--similar blow-up for super-critical exponents has also been found
for $h_t = -(h^3 h_{xxx})_x -\mathcal{B} (h^6 h_x)_x$ (presented by Bertozzi
and Pugh at the APS Division of Fluid Dynamics meeting, November
1997).
% \noindent $\bullet$ {\it Odd perturbations.}
%
% The odd perturbation $10^{-3} h_{\rm ss}'$ yielded an evolution qualitatively
% like that shown in Figure~\ref{q4_blow}. The solution appears to blow up in
% finite time, with the location of the maximum moving in time.
%
% \noindent $\bullet$ {\it Random perturbations.}
%
% Random perturbations led to the same type of dynamics as seen with
% $\pm h_{\rm ss}''$ perturbations.
%
% Checked with six random perturbations of magnitude .001
% \subsubsection{$q=4$, perturbing the constant steady state}
%
% Our numerics robustly confirm our predictions: for $\mathcal{B}
% \overline{h}^{q-1}<1$, small zero-mean perturbations of $\overline{h}$
% yield solutions that relax to the constant steady state, while for $\mathcal{B}
% \overline{h}^{q-1}>1$ such perturbations yield solutions that appear
% to blow up in finite time.
% Checked with perturbations of magnitude -.001, cos(x) and four random
% perturbations apiece.
%
% $h=2$ blows up
% $h=.5$ relaxes
\section{Conclusions and Future Directions}
\label{conclusions}
We have numerically studied the evolution equation $h_t = -(h^n
h_{xxx})_x - \mathcal{B} (h^m h_x)_x$. Our work suggests that the energy landscape
through which solutions travel is fairly simple and that understanding
the relative energy levels of the steady states gives considerable
insight into the dynamics of the solutions. In particular we have
found strong evidence for the existence of heteroclinic connections
between steady states, connections we conjecture to exist based on our
theorems \cite{LP3} on the relative energy levels of steady states. It
is an open problem to prove analytically that these connections exist.
In \S\ref{mobility} we presented numerical results on the persistence
of heteroclinic connections under changes in the mobility parameters
$n$ and $m$. We changed $n$ and $m$ in a way that preserved
$q=m-n+1$, preserving the energy $\mathcal{E}$ and hence the steady states and
energy landscape. But the timescale of the dynamic solution did change
noticeably in response to changes in the mobilities, even though the
shape of the solution changed little.
We would expect these structural stability observations to continue
to hold if the mobilities were changed in a way that, while not fixing
$q$, perturbed it only a little (so that the energy landscape is also
perturbed only a little).
In \S\ref{mobility} we further investigated critical mobility
exponents, such as the critical $n$ above which solutions remain
positive for all time (in other words, the critical exponent for film
rupture or pinch-off). An interesting question for the future is to
find formulas for the critical mobility exponents. These critical
exponents determine important qualitative features of the evolution,
and determining them would shed considerable light not only on the
equation studied here but also on related equations that arise from
physical models.
Lastly, in \S\ref{simulations} we demonstrated that the physical quantities
$P$, $A$, $\mathcal{B}$, and $m-n$ appear to fully determine the
large-scale features of the evolution, since they determine the steady
states and the energy landscape via the bifurcation parameter $E =
\mathcal{B} A^{q-1} P^{3-q}$. Finer details, such as the motion of the
extrema in time, depend on the initial data. Also, we have not been able to predict the amount of translation of a solution that might occur in the long-time limit.
Another open problem is to determine precisely the long-time limit of an
evolution that approaches a steady droplet configuration. In part the difficulty is that translates of steady droplets are also steady droplets, so that even if one knows the length and area of each droplet in the
configuration, their locations relative to each other must still be
determined. Thus droplet attractors might have rather high dimension.
In conclusion, we hope our numerical investigations of the power law
evolution (\ref{evolve}) will provide resources, ideas and motivation
for researchers studying $h_t = - (f(h)h_{xxx})_x - (g(h) h_x)_x$ with
non-power law coefficient functions $f$ and $g$. Some such
numerical studies exist already. For example, the papers \cite{Gruen00,OB}
consider an $f$ that is degenerate ($f(0)=0$) and $g$'s that are not
power laws, and there is of course a large literature on the
(non-degenerate) Cahn--Hilliard equation.
\subsection*{{\bf Acknowledgments}}
R. S. Laugesen was partially supported by NSF grant number DMS-9970228,
by a grant from the University of Illinois Research Board, and by a
fellowship from the University of Illinois Center for Advanced Study.
He is grateful for the hospitality of the Department of Mathematics at
Washington University in St. Louis.
M. Pugh was partially supported by NSF grant number DMS-9971392, by the
MRSEC Program of the NSF under Award Number DMR-9808595, by the ASCI
Flash Center at the University of Chicago under DOE contract B341495,
and by an Alfred P. Sloan fellowship. The computations were done using
a network of workstations paid for by an NSF SCREMS grant,
DMS-9872029. Part of the research was conducted while enjoying the
hospitality of the Mathematics Department and the James Franck
Institute of the University of Chicago.
M. Pugh thanks Todd Dupont and Bastiaan Braams for illuminating
conversations regarding numerical issues.
\appendix
%
\section{Numerical methods and parameters used in simulations}
\label{numerics}
%
The numerical simulations are done using a finite-difference evolution
code. Throughout this section the exponent $\ell$ denotes
the $\ell^{\rm th}$ time-step and $h^\ell$ is the numerical
approximation of the solution at time $t_\ell = \ell \Delta t$. The diffusion coefficients are represented as functions $f$ and $g$.
\subsection{The evolution code}
\label{the_code}
%
We use an adaptive time--stepping scheme (discussed in \S\ref{acc}) based on a Crank-Nicolson scheme:
%
\begin{eqnarray} \nonumber %\label{crank_nic}
\frac{h^{\ell+1}-h^\ell}{\Delta t}
&=& - \frac12 \left(f(h^{\ell+1/2})h^{\ell+1}_{xxx}\right)_x
- \frac12 \left(f(h^{\ell+1/2})h^{\ell}_{xxx}\right)_x \\
&& - \frac12 \left(g(h^{\ell+1/2})h^{\ell+1}_{x}\right)_x
\nonumber - \frac12 \left(g(h^{\ell+1/2})h^{\ell}_{x}\right)_x.
\end{eqnarray}
%
The diffusion coefficients are evaluated at $h^{\ell+1/2}$, which we
find by linearly extrapolating the solutions $h^{\ell-1}$ and $h^\ell$
to time $t_\ell + \Delta t/2$. The local truncation error is $O(\Delta t^3)$ in
time intervals where $\Delta t$ is fixed and is $O(\Delta t^2)$ at times when
the timestep is changed.
Finding $h^{\ell+1}$ reduces to solving a linear problem, which we
write in a residual formulation:
\begin{equation} \label{resid_crank_nic}
\mathcal{L} z = - \Delta t \left[ f(h^{\ell+1/2}) (h^{\ell}_{xxx} +
r(h^{\ell+1/2}) h^{\ell}_x) \right]_x
\end{equation}
where $z = h^{\ell+1}-h^{\ell}$, $r(y)=g(y)/f(y)$, and the linear
operator $\mathcal{L}$ is
\[
\mathcal{L} z := z
+ \Delta t \; \frac12 \left[ f(h^{\ell+1/2})(z_{xxx} + r(h^{\ell+1/2}) z_x ) \right]_x.
\]
The scheme uses $h^{\ell-1}$ and $h^\ell$ to compute $h^{\ell+1}$; for
the first step we take $h^{-1} = h^0$.
Performing a linear stability analysis about a constant steady state,
we find the scheme is stable. Specifically, perturbations whose wave
numbers are outside the unstable band of the linearized PDE do not
grow.
% We now perform a linear stability analysis of the scheme about the
% constant steady state $h \equiv \overline{h}$, for power law coefficients
% $f(y) = y^n$ and $g(y) = \mathcal{B} y^m$. The linearized equation is $h_t = -
% \overline{h}^{\: n} h_{xxxx} - \mathcal{B} \overline{h}^{\: m} h_{xx}$. Initial data $h^0(x) = \epsilon
% \cos(k x)$ yields $h^{\ell +1} = \sigma(k)^{\ell+1} \; \epsilon \cos(k x)$
% where
% \[
% \sigma(k) =
% \frac{1 - \frac{\Delta t}{2} \overline{h}^n k^2 (k^2 - \mathcal{B} \overline{h}^{m-n})}
% {1 + \frac{\Delta t}{2} \overline{h}^n k^2 (k^2 - \mathcal{B} \overline{h}^{m-n})} =:
% \frac{1-\mu(k)}{1+\mu(k)}.
% \]
% Hence there is a band of unstable modes: if $0 < k^2 < \mathcal{B} \overline{h}^{m-n}$
% ({\it i.e.} $\mu(k) < 0$) then $|\sigma(k)| > 1$ and the initial data
% $h_0(x) = \epsilon \cos(k x)$ yields a solution that grows exponentially in
% time. We consider a numerical scheme linearly stable if perturbations
% {\it outside} this band are not amplified:
% \[
% k^2 \geq k_c^2 := \mathcal{B} \overline{h}^{m-n} \Longrightarrow
% | \sigma(k) | \leq 1.
% \]
% It follows immediately from the form of the growth factor that the
% Crank-Nicolson scheme is linearly stable, as expected. Such a linear
% stability analysis provides a useful guide, though it is directly
% relevant only for small perturbations of flat steady states. \\
For the spatial discretization, the key issue is to implement the
scheme in a way that preserves steady states. A steady state $h_{\rm ss}$
satisfies $(h_{\rm ss})_{xxx} + r(h_{\rm ss}) (h_{\rm ss})_x = 0$ \cite{oron92,LP1}. Such an
`analytic steady state' will be $O(\Delta x^2)$-close to the
`finite-difference' steady state $\tilde{h}$ satisfying the following
discretization:
%
\begin{equation} \label{fd_ss}
\tilde{h}_{i+2}
- 3 \tilde{h}_{i+1}
+ 3 \tilde{h}_{i}
- \tilde{h}_{i-1}
+ \Delta x^{\: 2} \frac{ r(\tilde{h}_{i+1})+ r(\tilde{h}_{i})}{2}
\left( \tilde{h}_{i+1} - \tilde{h}_{i} \right) = 0 ,
\quad i=1\dots N.
\end{equation}
%
The meshpoints are $x_1 = \Delta x, x_2 = 2 \Delta x ,\cdots ,x_{N} = P$, where
$P$ is the length of the interval, and we denote the function values
at the meshpoints with subscripts: $\tilde{h}_1, \tilde{h}_2, \dots
\tilde{h}_{N}$. The function is periodic: $\tilde{h}_0 =
\tilde{h}_{N}$.
To implement the residual formulation (\ref{resid_crank_nic}), we
apply the $O(\Delta x^2)$ approximation
%
\begin{eqnarray}
\left[ \left(f(h) (z_{xxx} + r(h) z_x) \right)_x \right]_i
& \simeq &
\frac{f_{i+}}{\Delta x}\left(\frac{z_{i+2}-3z_{i+1}+3z_i-z_{i-1}}{\Delta x^3}
+ r_{i+}\frac{z_{i+1}-z_i}{\Delta x} \right) \nonumber \\
&&
- \frac{f_{i-}}{\Delta x}\left(\frac{z_{i+1}-3z_{i}+3z_{i-1}-z_{i-2}}{\Delta x^3}
+ r_{i-}\frac{z_{i}-z_{i-1}}{\Delta x} \right) , \label{resid}
\end{eqnarray}
%
where the subscripts $i+$ and $i-$ denote the right-average and
left-average:
\[
r_{i+} := \frac{r(h_{i+1})
+ r(h_i)}{2} , \quad \quad
r_{i-} := \frac{r(h_{i})
+ r(h_{i-1})}{2} .
\]
This approximation yields a $O(\Delta x^2)$-accurate $N \times N$ matrix
approximation $\tilde{\mathcal{L}}$ of the operator $\mathcal{L}$:
\[
\tilde{\mathcal{L}} \vec{z} = \Vec{RHS}(h^{\ell-1},h^{\ell}) .
\]
$\tilde{\mathcal{L}}$ is pentadiagonal periodic. The righthand side of
(\ref{resid_crank_nic}) is discretized analogously.
It follows immediately from (\ref{fd_ss}) that the time-stepping
scheme preserves finite-difference steady states. Also, by factoring
$f$ out we have also isolated the pressure gradient term in the
equation, $(h_{xx} + \mathcal{B}/q \; h^q)_x$.
%
\subsection{Computing the finite-difference steady state}
\label{compute_id}
%
Here, we only discuss the case of power law coefficients $f(h) = h^n$
and $g(h) = \mathcal{B} h^m$. Given $N$ uniformly distributed meshpoints
between $0$ and $P= 1$, we seek a finite-difference steady state
$\tilde{h}$ that solves the $N$ equations (\ref{fd_ss}) to the level
of round-off error.
% allowing us to add our own perturbations to it above the level of
% round-off.
We do this with Newton--Raphson iteration. To do this, we need a good
first guess for $\tilde{h}$. In the following, we describe how we
find a first guess and then how we execute the iteration.
Given the exponent $q \neq 0$ we first compute the rescaled steady
state $k=k_\alpha$ at the $N$ points $x=P/N, 2P/N, \cdots, P$. As
described in \cite[\S6.1]{LP1}, we do this by viewing the steady state
equation
%
\begin{equation} \label{k_ss_eqn}
k_{xx} + \frac{k^q - 1}{q} = 0, \quad \quad k(0) = \alpha, \quad
k_x(0) = 0 ,
\end{equation}
%
as an initial value problem in $x$. We verify that $k$ is spectrally
accurate by using a discrete fast Fourier transform to check that the
power-spectrum is fully resolved (cf.\ \S\ref{point_double}). For
numerical purposes, $k$ is an exact solution of equation
(\ref{k_ss_eqn}), and so in the following we refer to it as an `analytic
steady state'.
Once the analytic steady state $k$ is known, we rescale it to find an
analytic steady state $h_{\rm ss}$ of period $1$ and area $1$. Specifically,
if $P$ and $A$ are the period and area of $k$, then we use
the rescalings \cite[(3.3),(5.1)]{LP1} to determine
$$
\mathcal{B} = A^{q-1} P^{3-q} \quad \mbox{and} \quad
h_{\rm ss} = \frac{P}{A} k.
$$
We note that $k$, and hence $A=A(\alpha)$ and $P=P(\alpha)$, are determined by the
initial value $\alpha$ in (\ref{k_ss_eqn}). Varying $\alpha$, one
finds a family of $1$-periodic steady states of area $1$ each with a
different amplitude and satisfying the steady state equation with a
different value of $\mathcal{B}$. This was used in constructing the
bifurcation diagrams in Figure~\ref{bifurcs}.
To find the finite-difference steady state $\tilde{h}$ close to $h_{\rm ss}$
we need to solve the $N$ equations (\ref{fd_ss}), which we write as
\[
\vec{F}(\tilde{h}) := M \tilde{h} + \vec{V}({\tilde{h}}) = \vec{0}
\]
where $M$ is a tetradiagonal periodic matrix and
$\vec{V}(\tilde{h})_i$ is a nonlinear function of $\tilde{h}_i$ and
$\tilde{h}_{i+1}$. The Newton--Raphson iteration is
\[
\tilde{h}^{new} = \tilde{h}^{old}
- (D \vec{F}(\tilde{h}^{old}))^{-1} \vec{F}(\tilde{h}^{old}).
\]
To iterate, one has to solve $D \vec{F}\vec{x} =
\vec{F}(\tilde{h}^{old})$. We find that $D \vec{F}$ is a singular
matrix of rank $N-1$. We solve $D \vec{F} \vec{x} =
\vec{F}(\tilde{h}^{old})$, using the singular value decomposition of
$D\vec{F}$ obtained using LAPACK's `dgesvd.f' to solve for $\vec{x}$.
The iteration is then started with the initial guess $h_{\rm ss}$ and stopped
when the largest error in the $N$ equations (\ref{fd_ss}) is less than
$10^{-14}$.
\vskip 6pt
{\it Note:} For $q=1$, the finite--difference steady states satisfy a
linear problem:
\[
\tilde{h}_{i+2}
- 3 \tilde{h}_{i+1}
+ 3 \tilde{h}_{i}
- \tilde{h}_{i-1}
+ \Delta x^{\: 2} \tilde{\mathcal{B}}
\left( \tilde{h}_{i+1} - \tilde{h}_{i} \right) = 0 ,
\quad i=1\dots N.
\]
Analytic steady states are $a + b \cos(\sqrt{\mathcal{B}} x + \phi)$ and are
$1$-periodic if $\sqrt{\mathcal{B}}$ is $n 2 \pi$ for some integer $n$.
Sampling such a steady state on a uniform mesh gives
\[
\tilde{h}_j = a + b \cos(\sqrt{\mathcal{B}} \; j \Delta x + \phi ),
\]
and one can check that $\tilde{h}$ is a finite--difference steady state provided
\[
\tilde{\mathcal{B}} = 2 \; \frac{1-\cos(\sqrt{\mathcal{B}} \Delta x)}{\Delta x^2}
= \mathcal{B} - O(\Delta x^2).
\]
That is, there are nontrivial analytic steady states for a countable
collection of Bond numbers, $\sqrt{\mathcal{B}} \in \{ 2 \pi,4 \pi,6 \pi,
\ldots \}$, and nontrivial finite difference steady states for a
nearby countable set of Bond numbers $\tilde{\mathcal{B}}$.
\subsection{Timestepping and accuracy}
\label{acc}
The adaptive timestepping controls the accuracy as follows. An error
tolerance is set, $\epsilon = 10^{-11}$. At each time step, we first use
the Crank-Nicolson scheme to compute $h_1$, an approximation of the
solution at time $t + \Delta t$. We then take two timesteps with $\Delta t/2$
to compute $h_2$, another approximation of the solution at time
$t+\Delta t$. For some constant $C$, the error is
bounded \cite[\S5.2]{Iserles} by the difference of $h_1$ and $h_2$:
\[
\|h(\cdot,t+\Delta t) - h_2 \|_\infty \leq C \| h_1-h_2 \|_\infty.
\]
If $\|h_1-h_2\|_\infty> \epsilon$ then we replace $\Delta t$ with $\Delta t/2$ and try
again (without advancing in time). If $\|h_1-h_2\|_\infty < \epsilon/10$
then we replace $\Delta t$ with $2 \Delta t$ and try again. If
$\|h_1-h_2\|_\infty$ lies between $\epsilon$ and $\epsilon/10$ then we just take the
solution at time $t+\Delta t$ to be $h_2$. Admittedly, since we do not
know the constant $C$ this error control is valid only as long as $C$
is not large.
% Adaptive timestepping takes at least three times longer than using the
% Crank-Nicolson scheme with a fixed timestep. For this reason, we
% performed all of our exploratory studies using a fixed time-step.
% Once we found phenomena of interest, we re-ran using adaptive
% time--stepping. Since the first time--step has $O(\Delta t^2)$ local
% truncation error, rather than $O(\Delta t^3)$, the adaptive time--stepper
% initially refines $\Delta t$ to meet the tolerance. Also, most runs had an
% initial fast transient (see \S\ref{qm3pp}) which required early
% refinement of the timestep.
In practice, we find that the timestep is initially reduced to resolve
the initial fast transient and then coarsens. After this, the
timestep is rarely reduced, except near times when the run has to be
stopped anyway in order to point--double to resolve the growing
derivatives of the solution.
\subsection{Stopping criteria and issues for singularities} \label{stop}
%
To test whether to stop the code, we compute the minimum value of
$h^\ell$ at each time-step. If this minimum is less than or equal to
zero, then we stop the code.
% As discussed in \S\ref{mobility}--\ref{simulations}, we find that the
% stopping criterion is often met.
This is a natural stopping criterion since for any equation in flux
form $h_t + (h U)_x = 0$, if $h_{min}(t) \downarrow 0$ then $U_x
\to \infty$ at the location of $h_{min}$. Because $U$ is
determined by $h$ and its derivatives, this blow-up signals a loss of
smoothness in $h$.
Stopping of the code thus suggests a finite-time singularity, but we
emphasize that the code was designed to preserve the periodic steady
states and not to carefully resolve singularities. Hence in reality
the singularity may occur in infinite time, with the stopping
criterion being met in finite time because of oscillations in the
profile (possibly caused by loss of spectral resolution).
%The stopping criterion was always met after the solution became spectrally
%unresolved.
Also, since the code has no local mesh refinement, we have to
over-resolve much of the solution in order to resolve it near the
singular points, where it is tending to zero. This over--resolution
slows the computation significantly. To examine fine details of the
temporal and spatial scales of the singularities, then, one should
implement a code that has an adaptive spatial mesh.
% Loss of positivity immediately stops our code --- a positivity
% preserving code would keep running, having to be stopped when the
% solution is no longer spectrally resolved. Also, since the code has
% no local mesh refinement, we have to over-resolve much of the solution
% in order to resolve the solution where it is tending to zero. This
% over--resolution away from the singular points slows the computation
% significantly. To find fine details of the temporal and spatial
% scales of the singularities, one should implement a code that
% preserves positivity and has an adaptive spatial mesh.
If a smooth solution touches down in finite time, it may continue to
evolve as a weak solution. Weak solutions also arise when the initial
data has compact support; the degeneracy of the equation can provide
finite-speed of propagation of the support. Computing the evolution
of weak solutions is truly nontrivial. We refer interested readers to
\cite{barrett-blowey,Gruen99,LiyaAndrea}.
\subsection{Point-doubling and keeping solutions spectrally resolved}
\label{point_double}
In \S\ref{qm3pp}, we present a simulation with initial data $h_{\rm ss} -
.0001 h_{\rm ss}''$ for the evolution equation with $q=-3$ and $n=3$.
There, we find that the solution appears to touch-down in finite time.
As the minimum height decreases, the curvature increases, requiring
that after some time the number of meshpoints be increased to keep the
solution spectrally resolved. We do this as follows. We compute the
solution with $2,\!048$ meshpoints until the computation stops because
positivity is lost. We look at the power spectrum of the solution and
choose a time right before the active part of the power-spectrum is
reaching the Nyquist frequency. That is, we find the last time at
which the $1,\!023$rd Fourier amplitude of the solution is at the
level of round--off. We take the solution at this time and compute
its Fourier coefficients, defined for wave numbers $-\text{\scriptsize N}/2+1 \leq k
\leq \text{\scriptsize N}/2-1$ where $\text{\scriptsize N} = 2,\!048$. We pad by zeros, extending the
Fourier coefficients to be defined for wave numbers $-\text{\scriptsize N}+1 \leq k \leq
\text{\scriptsize N}-1$, and then compute the inverse Fourier transform. This yields a
function on $2\text{\scriptsize N}$ meshpoints that is indistinguishable from the
solution at that time, to the level of round--off. Using this
function as initial data, we continue the computation on $2 \text{\scriptsize N}$
meshpoints, repeating this point--doubling process whenever the
solution becomes unresolved.
\subsection{The simulations of Sections~\ref{mobility} and
\ref{simulations}}
\subsubsection{Notes on the $q=5/2$ simulations in \S\ref{mobility}}
\label{mobility_q2.5_params}
% alpha=0.21447368421053
% A = 5.66300993049705, P = 6.8688146976484, B = 1, E = 35.31928352164392
%
% A = 4.33452514815264, P = 2*pi, B = 1.56138376859101, E = 35.31928352164456
We computed a $1$-periodic steady state $h_{\rm ss}$ on $256$ meshpoints with
Bond number $\mathcal{B} = 35.32$ and area $1$, by rescaling $k_\alpha$ with
$\alpha=0.2145$, $P = 6.869$, and $A = 5.663$. All of the
simulations shown in Figures \ref{find_nc}-\ref{mobility_q2.5_n1} were
run until the $8,\!192$ meshpoint simulation lost resolution, except
for the $n=1.6625$ simulation which required $65,\!536$ meshpoints to
resolve the solution when $h_{min}$ was at its smallest. All the
profiles shown in Figure~\ref{one_two} are the final resolved solution
with $8,\!192$ meshpoints.
We now prove our earlier claim that for $q=5/2$ and $\mathcal{B} = 35.32$,
there cannot be two disjoint zero contact angle steady droplets in an
interval of length $1$, if the total area of the droplets is $1$.
There does exist a {\em single} zero angle droplet steady state with
that area and with length less than $1$ --- its length is $P =
\left( E_0(5/2)/\mathcal{B} \right)^2 = .8414 < 1$, where $E_0(5/2) = 32.40$ by
\cite[\S 3.1.2]{LP1}.
% P = 0.84143702864238
% E_0 = 32.39833523925486
Hence there is a zero contact angle droplet to which the solution on
the right of Figure~\ref{mobility_q2.5_n1} could relax. But if there
were a pair of zero contact angle steady state droplets, with areas
$A_1 = \lambda $ and $A_2 = (1-\lambda) $, then the combined length of
the two droplets would be
\[
P_1 + P_2 = \left( \lambda^{-3} + (1-\lambda)^{-3} \right)
\left( E_0(5/2)/\mathcal{B} \right)^2.
\]
The righthand side is a convex function of $\lambda$, achieving its
minimum value $13.46$ at $\lambda = 1/2$. This minimum value is
% 13.46299245827810
greater than $1$, and so one cannot fit the two droplet steady states
in an interval of length $1$. Thus the simulations described earlier
cannot be converging to a pair of steady zero contact angle droplets.
Indeed, our computations show the proto-droplet is slowly draining.
\subsubsection{Notes on the $q=1/2$ simulations in \S\ref{mobility}}
\label{mobility_q.5_params}
% alpha=0.21447368421053
% A = 6.77924328557084, P = 6.1684107099798, B = 1, E = 36.29464864663328
%
% A = 7.16474894335479, P = 2*pi, B = 0.981733055514273, E = 36.29464864663330
We computed a $1$-periodic steady state $h_{\rm ss}$ on $256$ meshpoints with
Bond number $\mathcal{B} = 36.29$ and area $1$, by rescaling $k_\alpha$ with
$\alpha=0.2145$, $P = 6.168$, and $A = 6.779$.
We now prove our earlier claim that for $q=1/2$ and $\mathcal{B} = 36.29$ there
can be two disjoint zero contact--angle droplets in an interval of
length $1$, if their total area is $1$. A single zero contact angle
steady state would satisfy \cite[\S 3.1.2]{LP1}
\[
P = A^{1/5} \left( E_0(1/2)/\mathcal{B} \right)^{2/5},
\quad \quad \mbox{with $A = 1$}
\]
where $E_0(.5) = 32.86$.
%E_0 = 32.86335345030997
We find $P = .9611 < 1$.
%P = 0.96105382707448
Thus a single zero contact angle steady state is a potential
long--time limit. For two zero contact angle steady states, we find
\[
P_1 + P_2 = \left(\lambda^{1/5}+(1-\lambda)^{1/5}\right)
\left( E_0(1/2)/\mathcal{B} \right)^{2/5}.
\]
This is a concave function with its maximum at $\lambda = 1/2$. We
find that $P_1+P_2 < 1$ if $\lambda \leq 10^{-7}$ (approx.) and so
one {\em can} have two zero contact angle droplets --- but one of them
must be fairly small. For example, if $\lambda = 10^{-7}$ then the
length of the smaller droplet is $P_1 = 0.03826$.
%P_1 = 0.03826024198462
In our $n=2$ simulation, the distance between the local minima was
$.03869$ when we stopped the code. This is close to $.02826$; thus a
two--droplet long--time limit is at least a possibility.
\subsubsection{Notes on the $q=-3$ simulations in \S\ref{qm3pp}}
\label{qm3pp_params}
% alpha =0.21447368421053
% A = 120.5845435244233, P = 16.318151714565, B = 1, E = 0.08930159682344
%
%A=6.88364316193449, P = 2*pi, B = 0.00325876775744170, E = 0.08930159682345
We computed a $1$-periodic steady state $h_{\rm ss}$ on $2048$ meshpoints with
Bond number $\mathcal{B} = .08930$ and area $1$, by rescaling $k_\alpha$ with
$\alpha=0.2145$, $P = 16.32$, and $A = 120.6$.
The steady state $h_{\rm ss}$ has large curvature at its local minima, and so
we need a large number of meshpoints to resolve the initial data $h_{\rm ss} +
.0001 h_{\rm ss}''$ with spectral accuracy. We find that for a solution on
$[0,1)$ we need $2,\!048$ meshpoints (following \S\ref{point_double}).
For the initial data $h_{\rm ss} - .0001 h_{\rm ss}''$, the resulting solution
appears to touch down in finite time, requiring point--doubling as the
solution's curvature increases (see \S\ref{point_double}.) In this
way, we computed a resolved solution to time $t = 56.0019$
at which time the $32,\!768$ meshpoint solution became spectrally
unresolved.
%
% I had a P = 2-pi periodic solution with area A = 6.88364316193449.
% it was for n=3 and became singular at time t = .0472496406249
%
% I rescaled this to a 1-periodic solution with area A = 1. so the new time
% t' = P^(n+4) A^(-n) t = 56.00187154016988
%
%
% for the odd perturbations, h + .0001 h', started with 2,\!048
% meshpoints and stopped with 32768 meshpoints at time
% t = .079286 (note: didn't use adaptive timestepping)
%
For the simulation with initial data $1 - .0002 \cos(2 \pi x)$ we
started with $256$ meshpoints and stopped at time $t =
31050.2$ with $16,\!384$ meshpoints.
%
% I had a P = 2-pi periodic solution with area A = .5*2*pi.
% it was for n=3 and became singular at time t = 2.49031918
%
% I rescaled it to a 1-periodic solution with area A = 1. so the new time
% t' = P^(n+4) A^(-n) t = 3.105020514666781e+04
%
\subsubsection{Notes on the $q=.5$ simulations in \S\ref{q.5}}
\label{q.5_params}
% alpha = 0.03000
% A = 7.12187265294875, P = 6.0493965410253, B = 1, E = 33.72741312406922
%
% A = 7.9799036525720, P = 2*pi, B = 0.962791362227189, E = 33.72741312406918
We computed a $1$-periodic steady state $h_{\rm ss}$ on $512$ meshpoints with
Bond number $\mathcal{B} = 33.73$ and area $1$, by rescaling $k_\alpha$ with
$\alpha=0.03000$, $P = 6.049$, and $A = 7.122$. The simulation
with initial data $h_{\rm ss} - .0001 h_{\rm ss}''$ stopped at
$t= 9525.89$ with $8,\!192$ meshpoints.
%
% I had a P = 2-pi periodic solution with area A = 7.9799036525720
% it was for n=1 and became singular at time t = 7.76254062500298
%
% I rescaled it to a 1-periodic solution with area A = 1. so the new time
% t' = P^(n+4) A^(-n) t = 9525.890391143095
%
\subsubsection{Notes on the $q=1$ simulations in \S\ref{q1}}
\label{q1_params}
We found a $256$-meshpoint finite--difference steady state for $\mathcal{B}
\sim 4 \pi^2$ (see \S\ref{compute_id}). We used perturbations of this
steady state for all three simulations in Figure~\ref{q1_fig}. The
simulation in Figure~\ref{q1_fig}b stopped at time $t= 1320.07$, with
$32,\!768$ meshpoints.
%
% I had a P = 2-pi periodic solution with area A = .7*2*pi
% it was for n=1 and became singular at time t = .59288999999979
%
% I rescaled it to a 1-periodic solution with area A = 1. so the new time
% t' = P^(n+4) A^(-n) t = 1320.065736757240
\subsubsection{Notes on the $q=3/2$ simulations in \S\ref{q1.5}}
\label{q1.5_params}
% alpha = 0.21447368421053
% A = 5.97468495749611, P = 6.4527831774677, B = 1, E = 40.06619622956337
%
% A = 5.51586207184591, P = 2*pi, B = 1.08318244359157, E = 40.06619622956358
We computed a $1$-periodic steady state $h_{\rm ss}$ on $256$ meshpoints with
Bond number $\mathcal{B} = 40.07$ and area $1$, by rescaling $k_\alpha$ with
$\alpha=0.2145$, $P = 6.453$, and $A = 5.975$.
The simulation with
initial data $h_{\rm ss} - .0001 \cos(x/2)$ stopped at
$t = 109594$ with $4,\!096$ meshpoints.
%
% I had a P = 2-pi periodic solution with area A = 5.51586207184591
% it was for n=1 and became singular at time t = 61.73049999999795
%
% I rescaled it to a 1-periodic solution with area A = 1. so the new time
% t' = P^(n+4) A^(-n) t = 109593.7376566941
% perturbing the constant: all started with 256 meshpoints. The
% perturbation of 2 was computed to 32768 meshpoints and stopped
% at time 12.6173825000052.
\subsubsection{Notes on the $q=1.768$ simulations in \S\ref{q1.768}}
\label{q1.768_params}
% alpha = .05
% A = 5.65959725101074, P = 6.7034114602923, B = 1, E = 39.458066992079015
% A = 4.66055351848993, P = 2*pi, B = 1.25724005342159, E = 39.458066992079015
%
% alpha = 0.50691271237156
% A = 6.11463784547698, P = 6.3879225342236, B = 1, E = 39.4580669920919859
% A = 5.81877310087267, P = 2*pi, B = 1.06019969639105, E = 39.4580669920919859
We now explain how we found the two steady states $h_{{\rm ss}1}$ and
$h_{{\rm ss}2}$ having period $1$ and area $1$ that are both steady
states for $\mathcal{B} = 39.46$. Let $E(\alpha) = A(\alpha)^{q-1}
P(\alpha)^{3-q}$. First we plot $E(\alpha)$ for $\alpha \in
(0,1)$. We seek $\alpha_1 < \alpha_2$ such that $E(\alpha_1) =
E(\alpha_2)$ with $E'(\alpha_1) < 0$ and $E'(\alpha_2) > 0$. Choosing
$\alpha_1 = .05$, we find that the desired $\alpha_2$ exists (see
\cite[Figure~5]{LP3}). To determine $\alpha_2$, we first compute
$k_{\alpha_1}$ and find its period $P(\alpha_1)= 6.703$, area
$A(\alpha_1) = 5.660$, and $E(\alpha_1)=39.46$. From the graph of
$E(\alpha)$, we make a rough estimate of $\alpha_2$ and choose six
values of $\alpha$ such that $E(\alpha) < E(\alpha_2) = E(\alpha_1)$
at the first three values and $E(\alpha) > E(\alpha_2)$ at the last
three values: then $\alpha_2$ is between the third and fourth values.
We interpolate $E$ at the six values of $\alpha$ with a quintic
polynomial and use Newton--Raphson iteration to find $\alpha_2=
0.5069$ satisfying $E(\alpha_1)=E(\alpha_2)$. We find $k_{\alpha_2}$
has period $P(\alpha_2) = 6.388$ and area $A(\alpha_2) = 6.115$. We
then rescale $k_{\alpha_1}$ to $h_{{\rm ss}1}$ and $k_{\alpha_2}$ to
$h_{{\rm ss}2}$. By construction, they are distinct positive periodic
steady states for $\mathcal{B} = 39.46$. All of the simulations were done with
$256$ meshpoints.
%
%alpha_1 = 0.05
%P_1 = 6.7034114602923
%A_1 = 5.65959725101074
%E_1 = 39.458066992079015
%alpha_2 = 0.50691271237156
%P_2 = 6.3879225342236
%A_2 = 6.11463784547698
%E_2 = 39.4580669920919859
%D_2 = 0.80095123794963329900
\subsubsection{Notes on the $q=5/2$ simulations in \S\ref{q2.5}}
\label{q2.5_params}
% alpha = 0.21447368421053
% A = 5.66300993049712, P = 6.8688146976484, B = 1, E = 35.31928352164458
%
% A = 4.33452514815265, P = 2*pi, B = 1.56138376859101, E = 35.31928352164468
We computed a $1$-periodic steady state $h_{\rm ss}$ on $256$ meshpoints with
Bond number $\mathcal{B} = 35.32$ and area $1$, by rescaling $k_\alpha$ with
$\alpha=0.2145$, $P = 6.869$, and $A = 5.663$. The simulation
with initial data $h_{\rm ss} - .0001 h_{\rm ss}''$ stopped at
$t = 101368$ with $2,\!048$ meshpoints.
%
% I had a P = 2-pi periodic solution with area A = 4.33452514815265
% it was for n=1 and became singular at time t = 44.8685
%
% I rescaled it to a 1-periodic solution with area A = 1. so the new time
% t' = P^(n+4) A^(-n) t = 101367.6470292231
\subsubsection{Notes on the $q=4$ simulations in \S\ref{q4}}
\label{q4_params}
% alpha = 0.21447368421053
% A = 5.54273628800614, P = 7.5356823758703, B = 1, E = 22.59696244299183
%
% A = 3.2128915511588, P = 2*pi, B = 4.28096554684008, E = 22.59696244299193
We computed a $1$-periodic steady state $h_{\rm ss}$ on $256$ meshpoints with
Bond number $\mathcal{B} = 22.60$ and area $1$, by rescaling $k_\alpha$ with
$\alpha=0.2145$, $P = 7.536$, and $A = 5.543$. The energy
$\mathcal{E}(h_0) = .3471$ is positive.
% E = 1.4446096218033D-02
%
% Now, I know that E_{1-periodic} = (2pi)^3/A^2 E_{2pi-periodic}
% and A = 3.2128915511588, E_{2pi-periodic}= .014446096218033 so
% E_{1-periodic} = 0.34713465746149
The simulation with initial data $h_{\rm ss} - .001 h_{\rm ss}''$ stopped $t =
30565.4$ with $4,\!096$ meshpoints. At that time, $h_{max}$ had
increased by a factor of $36.57$.
%
% since the solution was multiplied by P/A = 1.95561699084223, that
% factor of 18.7 becomes 36.57003772874964
%
% we had a P = 2-pi periodic solution with area A = 3.2128915511588
% that became singular at time 10.0282848199749 with n = 1
% We rescaled this to a 1-periodic solution with area A = 1. so the new time
% t' = P^(n+4) A^(-n) t = 30565.38956941903
% min, period,E 0.21447368421053 7.5356823758703, 1.4446096218033D-02
% area = 3.2128915511588
% also tried:
% min, period,E 0.31921052631579 7.2873679452764, -2.7767499031755D-02
% area = 3.6852468649675
% min, period,E 0.42394736842105 7.0499429260269, -7.5626160681070D-02
% area = 4.1879955407977
% min, period,E 0.52868421052632 6.8309057742999, -0.12836039016800
% area = 4.7046544295994
% min, period,E 0.63342105263158 6.6376361095339, -0.18368606958692
% area = 5.2100166883659
% min, period,E 0.73815789473684 6.4775703956385, -0.23729015529975
% area = 5.6690395361573
% min, period,E 0.84289473684211 6.3589681865284, -0.28246947824016
% area = 6.0357751862326
% min, period,E 0.94763157894737 6.2924018367795, -0.31015696322974
% area = 6.2524373141618
%
% and they all blew up...
\subsection{Finding random perturbations}
\label{rand_perts}
A function on $2,\!048$ meshpoints can resolve $1,\!023$ frequencies,
so we choose $\{a_k \}$ and $\{\phi_k \}$ to be two sets of $1,\!023$
uniformly distributed random numbers in $[0,1]$ and define the
perturbation
\[
\phi(x) = \sum_{k=1}^{1,\!023} a_k \exp(-.036 k) \cos(k 2 \pi x + \phi_k ).
\]
This perturbation has zero mean, and has random amplitudes and phases
at each wave--number. The decay rate $0.036$ is chosen so that the
amplitudes at wave numbers $k=900$ and higher are at the level of
round--off error, meaning the initial data is spectrally resolved. We
then normalized $\phi$ to have $L^\infty$ norm equal to 1. (In fact,
we normalized all perturbations in the article in this way.) Random
perturbations on other mesh-sizes are constructed analogously.
% \bibliographystyle{plain}
% \bibliography{films,hele}
\begin{thebibliography}{10}
\bibitem{barrett-blowey}
J.~W. Barrett, J.~F. Blowey, and H.~Garcke.
\newblock Finite element approximation of a fourth order nonlinear degenerate
parabolic equation.
\newblock {\em Numer Math}, 80(4):525--556, 1998.
\bibitem{BBW}
A.~J. Bernoff, A.~L. Bertozzi, and T.~P. Witelski.
\newblock Axisymmetric surface diffusion: dynamics and stability of
self-similar pinchoff.
\newblock {\em J. Statist. Phys.}, 93(3-4):725--776, 1998.
\bibitem{similar}
A.~L. Bertozzi, M.~P. Brenner, T.~F. Dupont, and L.~P. Kadanoff.
\newblock Singularities and similarities in interface flow.
\newblock In L.~Sirovich, editor, {\em Trends and Perspectives in Applied
Mathematics}, volume 100 of {\em Applied Mathematical Sciences}, pages
155--208. Springer--Verlag, New York, 1994.
\bibitem{BPLW}
A.~L. Bertozzi and M.~C. Pugh.
\newblock Long--wave instabilities and saturation in thin film equations.
\newblock {\em Commun Pur Appl Math}, 51:625--661, 1998.
\bibitem{BPFS}
A.~L. Bertozzi and M.~C. Pugh.
\newblock Finite-time blow-up of solutions of some long-wave unstable thin film
equations.
\newblock {\em Indiana Univ. Math. J.}, 49(4):1323--1366, 2000.
\bibitem{Beyn02}
Wolf-J{\"u}rgen Beyn, Alan Champneys, Eusebius Doedel, Willy Govaerts, Yuri~A.
Kuznetsov, and Bj{\"o}rn Sandstede.
\newblock Numerical continuation, and computation of normal forms.
\newblock In {\em Handbook of dynamical systems, Vol. 2}, pages 149--219.
North-Holland, Amsterdam, 2002.
\bibitem{bricmont}
J.~Bricmont, A.~Kupiainen, and J.~Taskinen.
\newblock Stability of {C}ahn-{H}illiard fronts.
\newblock {\em Commun Pur Appl Math}, 52(7):839--871, 1999.
\bibitem{carlen}
E.~A. Carlen, M.~C. Carvalho, and E.~Orlandi.
\newblock A simple proof of stability of fronts for the {C}ahn-{H}illiard
equation.
\newblock {\em Commun. Math. Phys.}, 224(1):323--340, 2001.
\bibitem{First}
P.~Constantin, T.~F. Dupont, R.~E. Goldstein, L.~P. Kadanoff, M.~J. Shelley,
and S.-M. Zhou.
\newblock Droplet breakup in a model of the \mbox{Hele--Shaw} cell.
\newblock {\em Phys Rev E}, 47(6):4169--4181, June 1993.
\bibitem{Doedel97}
E.~J. Doedel, M.~J. Friedman, and B.~I. Kunin.
\newblock Successive continuation for locating connecting orbits.
\newblock {\em Numer. Algorithms}, 14(1-3):103--124, 1997.
\newblock Dynamical numerical analysis (Atlanta, GA, 1995).
\bibitem{FifeGradOrig}
Paul~C. Fife.
\newblock Models for phase separation and their mathematics.
\newblock {\em Electron. J. Differential Equations}, pages No. 48, 26 pp.
(electronic), 2000.
\newblock Originally written in June 1991.
\bibitem{FifeReview}
Paul~C. Fife.
\newblock Pattern formation in gradient systems.
\newblock In {\em Handbook of dynamical systems, Vol. 2}, pages 677--722.
North-Holland, Amsterdam, 2002.
\bibitem{GPS97}
R.~Goldstein, A.~Pesci, and M.~Shelley.
\newblock Instabilities and singularities in {H}ele--{S}haw flow.
\newblock {\em Phys Fluids}, 10(11):2701--2723, 1998.
\bibitem{Gruen00}
G.~Gr\"un and M.~Rumpf.
\newblock Simulation of singularities and instabilities arising in thin film
flow.
\newblock {\em Eur J Appl Math}, 12(3) 3:293--320, 2001.
\bibitem{Gruen99}
G{\"u}nther Gr{\"u}n and Martin Rumpf.
\newblock Nonnegativity preserving convergent schemes for the thin film
equation.
\newblock {\em Numer. Math.}, 87(1):113--152, 2000.
\bibitem{Iserles}
A.~Iserles.
\newblock {\em A first course in the numerical analysis of differential
equations}.
\newblock Cambridge University Press, Cambridge, 1996.
\bibitem{LP1}
R.~S. Laugesen and M.~C. Pugh.
\newblock Properties of steady states for thin film equations.
\newblock {\em Eur J Appl Math}, 11(3):293--351, 2000.
\bibitem{LP2}
R.~S. Laugesen and M.~C. Pugh.
\newblock Linear stability of steady states for thin film and
{C}ahn--{H}illiard type equations.
\newblock {\em Arch Ration Mech An}, 154(1):3--51, 2000.
\bibitem{LP3}
R.~S. Laugesen and M.~C. Pugh.
\newblock Energy levels of steady states for thin film type equations.
\newblock {\em J Diff Eq}, 182(2):377--415, 2002.
\bibitem{Champneys99}
G.~J. Lord, A.~R. Champneys, and G.~W. Hunt.
\newblock Computation of homoclinic orbits in partial differential equations:
an application to cylindrical shell buckling.
\newblock {\em SIAM J. Sci. Comput.}, 21(2):591--619 (electronic), 1999.
\bibitem{Manneville}
P.~Manneville.
\newblock {\em Dissipative structures and weak turbulence}.
\newblock Academic Press Inc., Boston, MA, 1990.
\bibitem{myers}
T.~G. Myers.
\newblock Thin films with high surface tension.
\newblock {\em SIAM Rev.}, 40(3):441--462 (electronic), 1998.
\bibitem{OB}
A.~Oron and S.~G. Bankoff.
\newblock Dewetting of a heated surface by an evaporating liquid film under
conjoining/disjoining pressures.
\newblock {\em J Colloid Interface Sci}, 218(1):152--166, 1999.
\bibitem{oron}
A.~Oron, S.~H. Davis, and S.~G. Bankoff.
\newblock Long--scale evolution of thin liquid films.
\newblock {\em Rev Mod Phys}, 69(3):931--980, Dec 1997.
\bibitem{oron92}
A.~Oron and P.~Rosenau.
\newblock Formation of patterns induced by thermocapillarity and gravity.
\newblock {\em J Phys II}, 2:131--146, 1992.
\bibitem{Otto}
F.~Otto.
\newblock Lubrication approximation with prescribed non--zero contact angle: an
existence result.
\newblock {\em Commun Part Diff Eq}, 23:2077--2164, 1998.
\bibitem{TC}
J.~E. Taylor and J.~W. Cahn.
\newblock Linking anisotropic sharp and diffuse surface motion laws via
gradient flows.
\newblock {\em J. Statist. Phys.}, 77(1-2):183--197, 1994.
\bibitem{WD82}
M.~B. Williams and S.~H. Davis.
\newblock Nonlinear theory of film rupture.
\newblock {\em J Colloid Interf Sci}, 90(1):220--228, 1982.
\bibitem{WB99}
Thomas~P. Witelski and Andrew~J. Bernoff.
\newblock Stability of self-similar solutions for van der {W}aals driven thin
film rupture.
\newblock {\em Phys. Fluids}, 11(9):2443--2445, 1999.
\bibitem{WB00}
Thomas~P. Witelski and Andrew~J. Bernoff.
\newblock Dynamics of three-dimensional thin film rupture.
\newblock {\em Phys. D}, 147(1-2):155--176, 2000.
\bibitem{ZL99}
W.~W. Zhang and J.~R. Lister.
\newblock Similarity solutions for van der {W}aals rupture of a thin film on a
solid substrate.
\newblock {\em Phys Fluids}, 11(9):2454--2462, 1999.
\bibitem{LiyaAndrea}
L.~Zhornitskaya and A.~L. Bertozzi.
\newblock Positivity-preserving numerical schemes for lubrication-type
equations.
\newblock {\em SIAM J. Numer. Anal.}, 37(2):523--555 (electronic), 2000.
\end{thebibliography}
\end{document}