How to use NDSolve with moving boundary conditions?

Question

So I am trying to solve the movement in space and time of a spreading gravity current. The interface satisfies the following PDE:

$ \frac{\partial h}{\partial t} = \frac{\partial}{\partial x}\left(h^3 \frac{\partial h}{\partial x}\right) $.

The nose, $h(x,t) = 0$, of the current can move forward in space. If we say the initial condition is $h(x,0) = 1-x$, then the nose, $x_N$ is initially at $x = 1$. The base of the current can be fixed at $h(0,t) = 1$. So I have one IC and one BC, I need another BC to close the system. Now, as time moves forward the nose (fixed at h = 0) propagates according to the expression:

$\int_0^{x_N} h \, dx = \frac{1}{2} + t$

or kinematically $\dot{x_N} = -h^2 h_x$.

Clearly the second boundary condition in $x$ needs to come from the nose condition, but I have no idea how to input an integral boundary condition or the kinematic condition. Secondly, I have no idea how to deal with xmax in NDSolve as clearly my xmax is moving forward with each time step.

Here is my attempt, which gives a "There are fewer dependent variables error"

NDSolve[{D[h[x, t], t] == D[h[x, t]^3 D[h[x, t], x], x], h[0, t] == 1,
h[x, 0] == 1 - x (D[h[x, t], t] /. x -> 1) == -h[x, t]^2 (D[h[x, t], x] /. 
  x -> 1)}, h, {x, 0, 1}, {t, 0, 1}]

EDIT

So I solved this problem by amending the code in this Stefan problem post. I modified the statement of the problem slightly to work more closely to @ybeltukov excellent solution, and actually solved the following system:

$ \frac{\partial h}{\partial t} = \frac{\partial}{\partial x}\left(h^3 \frac{\partial h}{\partial x}\right) \\ \dot{s} = -h^2 h_x \\ s(0) = 0 \\ h(x,0) = \begin{cases}1 \quad \text{for } x = 0 \\ 0 \quad \text{otherwise}\end{cases} \\ h(s(t),t) = 0 \\ h(0,t) = 1 $.

I made the same change to a normalised variable for my PDE and then amended finite difference code:

n = 100;
\[Delta]\[Xi] = 1./n;

ClearAll[dv, t];
dv[v_List] := 
With[{s = First@v, u = Rest@v}, 
With[{ds = u[[-1]]^3/(s \[Delta]\[Xi]), \[Xi] = N@Range[n - 1]/n, 
 d1 = ListCorrelate[{-0.5, 0, 0.5}/\[Delta]\[Xi], #] &, 
 d2 = ListCorrelate[{1, -2, 1}/\[Delta]\[Xi]^2, #] &}, 
Prepend[3 u^2 d1[#]^2/s^2 + u^3 d2[#]/s^2 + \[Xi] ds d1[#]/s &@
  Join[{1}, u, {0.}], ds]]];
 s0 = 0.001;
 v0 = Flatten@Prepend[ConstantArray[0.001, n - 2], {s0, 1.}];
 sol = NDSolve[{v'[t] == dv[v[t]], v[0] == v0}, v, {t, 0, 50}][[1, 1, 
2]];

Returning from the normalised variable is identical code.

You can then make a nice time series with this code:

ListAnimate[
Table[Plot[u[t, x], {x, 0, 10}, PlotRange -> {{0, 10}, {0, 3}}], {t, 
 0, 50, 0.1}]]

Which returns a spreading current as desired.

Answer 1

The PDE can be easily solved numerically if a simple additional boundary condition is imposed.

Adding a second bundary condition h[1,t] = 0 we find numerically, without any error message

hh[x_, t_] = 
  h[x, t] /. 
   NDSolve[{D[h[x, t], t] == D[h[x, t]^3 D[h[x, t], x], x], h[0, t] == 1, 
      h[1, t] == 0, h[x, 0] == (1 - x)}, h[x, t], {x, 0, 1}, {t, 0, 5}][[1]];

Plot3D[hh[x, t], {x, 0, 1}, {t, 0, 3}, PlotRange -> {0, 1}, 
 PlotLabel -> "Solution of a PDE", AxesLabel -> {"x", "t", "h[x,t]"}]

Executing NDSolve without a secod boundary condition leads to about the same graph and the error message that a second boundary condition is missing.

The time dependence at various positions is given by the graph

Plot[{hh[0.2, t], hh[0.5, t], hh[0.8, t], hh[0.95, t]}, {t, 0, 3}, 
 PlotRange -> {0, 1}, 
 PlotLabel -> 
  "Solution of a PDE\nTime dependence for various locations x = 1 - h[x,0]", 
 AxesLabel -> {"t", "h[x,t]"}]

The asymptotic function for large t is given by

ha[x_] = (1 - x)^(1/4)

which solves the steady state ODL derived from the original PDE.

We could study the approch of the solution to the equlibirum by linearizing the equation, i.e. letting

h = hL[x] + ha[x]

and dropping all non linear terms in hL.

Partial solutions are given by Bessel functions. The interested reader will surely be able to follow this path on his own.

Moving structures

Let us have a look at moving structures allowed by the PDE.

Letting

h[x,t] = g[x - v t]

The PDE reduces to an ODE which can be solved

gg[x_] = g[x] /. 
   DSolve[g'[x] v == 3 g[x]^2 g'[x]^2 + g[x]^3 g''[x], g[x], x][[1]] /. 
  C[2] -> 0

(* Out[341]= InverseFunction[-((C[1]^3 Log[C[1] + v #1])/v^4) + (C[1]^2 #1)/v^3 - (
    C[1] #1^2)/(2 v^2) + #1^3/(3 v) &][x] *)

Defining

h1[s_] = (gg[x][[0, 1]] /. {C[1] -> -10, v -> 1} // FullSimplify)@s

(* Out[393]= 100 s + 5 s^2 + s^3/3 + 1000 Log[-10 + s] *)

Gives

Plot[Re[h1[s]], {s, 0, 10}, PlotRange -> {0, 3540}, 
 PlotLabel -> "'Nose' profile", AxesLabel -> {"x - t", "g"}]

This section is just a sketch which must be worked out more thoroughly.