EDIT #2
My error was useful. It brought me to the conclusion that the difficulties in solving the PDE of the OP are due to the drift term
$$\frac{\partial (x u(x,t))}{\partial x}$$
If the drift term is included, many boundary problems are ill defined. It turns out that there are cases where mathematically there is only a trivial solution u = 0 but numerically you obtain nice non trivial pictures.
Example (notice the linear dispersion term)
$$\frac{\partial u(x,t)}{\partial t}=\frac{\partial ^2u(x,t)}{\partial x^2}+\frac{\partial (x u(x,t))}{\partial x}$$
$$u(x=\pm 1,t)=0$$
$$u(x,0)=\cos ^2\left(\frac{\pi x}{2}\right)$$
I shall describe the details soon.
EDIT
After some days of consideration I found that with a drastic simplification we can find a satisfactory solution.
The simplification consists in dropping the "drift" term in the PDE.
The "energy leak paradoxon" which persists here, was solved.
Further discussion and revision of the original post are still to be done.
With an intial condition quit close to that of the OP
h[x_] := Exp[-18 x^2]
the solution is (notice the missing "drift" term)
uu[x_, t_] =
u[x, t] /.
NDSolve[D[u[x, t], t] == D[u[x, t] D[u[x, t], x], x] &&
u[-1, t] == 0 && u[+1, t] == 0 && u[x, 0] == h[x],
u[x, t], {x, -1, 1}, {t, 0, 10}][[1]];
and a 3D plot of the function
Plot3D[uu[x, t], {x, -1, 1}, {t, 0, 1}, PlotRange -> All,
AxesLabel -> {"x", "t", "u"}, LabelStyle -> {12, Bold},
ViewPoint -> 10/22 {Pi/2, -3 Pi/2 , \[Pi]},
PlotLabel ->
Style["Numerical solution of PDE\n\
\!\(\*SubscriptBox[\(\[PartialD]\), \(t\)]\)u = \!\(\*SubscriptBox[\(\
\[PartialD]\), \(x\)]\)(u \!\(\*SubscriptBox[\(\[PartialD]\), \
\(x\)]\)u)\nu(x=\[PlusMinus]1,t) = 0\nu(x,t = 0) = Exp[-(6 \
x\!\(\*SuperscriptBox[\()\), \(2\)]\)/2]\n", 14, Bold]]
There is no numeric instability seen here anymore.
This goes also for initial conditions of the form
h1[x_] = Cos[ Pi/2 x ]^k (* k=1,2,3,4 *)
The stationary state reached for large t is the trivial one, u = 0.
I took some days off in order to clarify the "mystery" of the energy leak which is in fact specific to this non linear problem because the function appears as a factor before the derivative. This is in contrast to the situation with the diffusion equation. The solution - in brief - is that u(x) u'(x) can remain finite for x->1 if the derivative diverges. As an example take f(x) = Sqrt(1-x^2) which gives f(x)f'(x) = - x. The numerical results for large enough t show large (divergent) spatial derivatives at the boundaries. Hence the "energy" will leak out indefinitely until it is zero.
The "energy" of the inital state is
w0 = Integrate[h[x], {x, -1, 1}]
(* Out[13]= 1/3 Sqrt[\[Pi]/2] Erf[3 Sqrt[2]] *)
The normalized "energy" as a function of time is then given by
w[t_] := Integrate[uu[x, t], {x, -1, 1}]/w0
Plots over different regions of time are
Plot[w[t], {t, 0, 1}, PlotRange -> {0, 1.1},
PlotLabel ->
Style["normalized 'energy' as a function of time\ninitial condition u(x,t = \
0) = Exp[-(6 x\!\(\*SuperscriptBox[\()\), \(2\)]\)/2]" , 14, Bold],
AxesLabel -> {"t", "w"}, LabelStyle -> {12, Bold}]
Original post
The following is the attempt of an in depth study of this interesting problem. Due to a lack of time here are only the first steps stating some of the specifiy problems encountered. I shall complete it subsequently.
The simplified problem
To begin with, let us simplify the problem without losing the main features.
We consider the interval -1 < x < 1 instead of the infinite interval, and take as inital condition a power of (1-x^2).
The PDE is
eq = D[f[x, t], t] == D[x f[x, t], x] + D[f[x, t] D[f[x, t], x], x]
(*
Out[6]=
Derivative[0, 1][f][x, t] ==
f[x, t] + x*Derivative[1, 0][f][x, t] + Derivative[1, 0][f][x, t]^2 +
f[x, t]*Derivative[2, 0][f][x, t]
*)
The initial distribution is taken as
f0[x_, k_] = (1 - x^2)^k;
The boundary conditions are
bound = {f[1, t] == 0 && f[-1, t] == 0};
The numerical solution is then (with k = 4, say)
sol = NDSolve[eq && f[x, 0] == f0[x, 4] && bound,
f[x, t], {x, -1, 1}, {t, 0, 20}][[1]];
ff[x_, t_] = f[x, t] /. sol;
The solution approaches the steady state solution in a few time steps:
pt = Table[
Plot[ff[x, t], {x, -1, 1}, PlotRange -> {{-1, 1}, {-1, 1}}], {t, 0, 4, 2}];
Show[pt]
The steady state solution g[x] is given by a solution of the ODE
eqs = x g[x] + g[x] g'[x] == const
(*
Out[23]= x g[x] + g[x] Derivative[1][g][x] == const
*)
The constant can be calculated from the boundary conditions or from the condition of symmetry g'[0]==0 to give const = 0.
Hence we find g[x] = 0 or, applying the boundary condition,
g[x_] := (1 - x^2)/2
Transistion to the steady state
The transisiton from the initial state to the steady state requires some discussion, as it is not as smooth as it might seem from the plots.
Indeed, a 3D plot with "continuous" time steps shows some instability at t ~= 0.15
Plot3D[ff[x, t], {t, 0, 1}, {x, -1, 1}]
Notice that this is the region in t where the solution makes the transisiton from positive second spatial derivative to the negative one of the steady state solution. This will be investiged later in more detail.
Conservation law
An interesting quantity is the integral w over the solution over the complete x-interval. This could be interpreted as the total thermal energy in the system
w [t_] := Integrate[f[x, t], {x, -1, 1}]
From the original PDE it might seem that this quantity must not change in time.
Indeed, integrating the PDE over x from -1 to +1 and changing the order of the temporal derivative and the spatial integral gives
D[w, t] = q[x, t]/x -> +1 - q[x, t]/x -> -1
where
q[x_, t_] = x f[x, t] + f[x, t] D[f[x, t], x];
But because ot the boundary conditions, q is zero.
On the other hand w assumes different values for t = 0 and t->[Infinity], so it must have changed in the meantime. This strange behaviour will also be examined in more detail, of course.
