Given that the plot of h[-(Sqrt[2]*π)/2, t]
looks something like this,

and that the differential equation has a singularity at h == 0
, one might expect that NDSolve
ought not to continue the solution to the max time the OP sets, which is t == 6
. NDSolve
might do it, but that would be because of numerical error when h
is close to zero (insufficient accuracy).
The solution is symmetric with respect to x == -(Sqrt[2]*π)/2
, the derivative $\partial h/\partial x$ along this line is zero. For the time integration along this line, therefore, the differential equation is equivalent to
$${\partial h \over \partial t} + {1 \over h}\; {\partial^2 h \over \partial x^2} + h^3\; {\partial^4 h \over \partial x^4} = 0\,.$$
When $h \approx 0$, assuming ${\partial^4 h / \partial x^4}$ does not grow too fast, we can reduce this to
$${\partial h \over \partial t} \approx - {1 \over h}\; {\partial^2 h \over \partial x^2}\,.$$
Thus the value of the time derivative will be very large and the step NDSolve
will be very small (the equation becomes stiff). If in taking a discrete step, NDSolve
pushes the value of h
across zero, the time derivative changes sign. The solution can then be expected to bounce back and forth across zero. With AccuracyGoal -> 10
, one might expect the solution to oscillate between $\pm 10^{-10}$ or so (see below). The step size at this point turns out to be around $10^{-20}$, which means it would take forever to carry out the pointless task of advancing much farther.
Therefore, one should use the setting AccuracyGoal -> Infinity
.
The other issue one might infer from the discussion above is that a fixed step size will severely limit Mathematica's ability to home in on the singularity. So the next recommendation is that the "FixedStep"
method should be discarded.
Another issue is finding a spatial grid that yields an acceptable error estimate. The OP used "MinPoints" -> 180
. I found "MinPoints" -> 101
was sufficient to avoid a warning, but "MinPoints" -> 201
produced a better looking graph. (At some point I raised the WorkingPrecision
to 40
. I do not think this is strictly necessary. It takes a long time to run a trial and such a difference in precision is relatively unimportant, so I did not experiment with it much.)
Here is the code I used. It took almost half an hour.
sol = NDSolve[
{D[h[x, t], t] + D[h[x, t]^-1*D[h[x, t], x], x] + D[h[x, t]^3*D[h[x, t], {x, 3}], x] == 0,
h[-((3 Sqrt[2] π)/2), t] == h[(Sqrt[2] π)/2, t],
h[x, 0] == 1 + 1/10*Sin[x/Sqrt[2]]},
h,
{x, -((3 Sqrt[2] π)/2), (Sqrt[2] π)/2},
{t, 0, 6},
Method -> {"MethodOfLines", "SpatialDiscretization" -> {"TensorProductGrid",
"MinPoints" -> 201, "MaxPoints" -> 1001, "DifferenceOrder" -> "Pseudospectral"}},
AccuracyGoal -> Infinity, WorkingPrecision -> 40];
NDSolve`Iterate::ndsz
: At t == 4.08475454101253896091960366504686029429278350495484420555092981315002745816884`40.
, step size is effectively zero; singularity or stiff system suspected.
We can check the domain and get the value of t
reported in the warning message:
(h["Domain"] /. First[sol])[[2, 2]]
(* 4.08475454101253896091960366504686029429278350495484420555092981315002745816884`40. *)
We can see that the computed value of h[-(Sqrt[2]*π)/2, t0]
is very small. The derivative is quite large, so Newton's method estimates that the zero is within a little under 10^-36
of the last value calculated by NDSolve
.
critpt = {x -> (-(Sqrt[2]*π)/2), t -> (h["Domain"] /. First[sol])[[2, 2]]};
h[x, t] /. First[sol] /. critpt // N
D[h[x, t], t] /. First[sol] /. critpt // N
%% / %
(*
7.16714*10^-18 -- value of h
-1.97906*10^19 -- t-derivative
-3.62148*10^-37 -- error estimate for t (Newton's method)
*)
We can see that the spatial grid has 201 points and the time integration took almost 250 steps. The last nine steps were very small indeed, around 10^-36
to 10^-37
.
(h["Grid"] /. First[sol]) // Dimensions
(h["Grid"] /. First[sol])[[1, -10 ;;, 2]] // Differences
(*
{201, 246, 2}
{1.419*10^-36, 2.218*10^-36, 1.609*10^-36, 1.504*10^-36, 1.96*10^-37,
1.96*10^-37, 3.07*10^-37, 2.22*10^-37, 2.08*10^-37}
*)
Let's check that the value of h
has not crossed zero. The midpoint of the spatial grid is position 101
and corresponds to x == -(Sqrt[2] π)/2
. The following returns the positions of any negative values along this line; there are none:
Position[(h["ValuesOnGrid"] /. First[sol2])[[101]], _?Negative]
(* {} *)
Here is a look at the solution. One can see oscillations beginning to form at the time front, which suggests the error oscillations might be getting large.
Plot3D[h[x, t] /. First[sol],
{x, -((3 Sqrt[2] π)/2), (Sqrt[2] π)/2},
Evaluate@Prepend[(h["Domain"] /. First[sol])[[2]], t],
PlotRange -> {0, All}, WorkingPrecision -> 40
]

Caveat
The value of t
where h
is zero depends on the accuracy of the computed solution. I carried out the solution with "MinPoints" -> 101
. Mathematica did not complain about the solution, but the stopping point was about 4.08323
, which differs from the above stopping point by about 0.0015
. That is a rather large difference. The value of h
was 6.34605*10^-18
, which is about the same. One would also infer a similar estimate for the error in t
(4*10^-37
) from Newton's Method. Probably there is an accumulation of error throughout the time integration. One might consider doubling the spatial grid until a semblance of convergence occurs; perhaps WorkingPrecision
might need to be increased. I would expect that to take some computation time. It seems feasible, but it is more time than I have to devote to this problem.
Verification of oscillations when h ≈ 0
The following was produce with a spatial grid of 101
points and an AccuracyGoal
of 10
. The first few hundred points (beyond the plot range) track the solution. For the last 9000+ points, the steps oscillate between ±2*10^-8
.
