\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}