# How to propel the integration of time a little bit further? Numerical solution can not evolve to the max time

I try to solve a nonlinear partial differential equation. I obtain a numerical solution which can not continue to the max time I set, I always receive message

NDSolve::ndcf: Repeated convergence test failure at t ==...,

and

NDSolve::eerr: Warning: scaled local spatial error estimate of ...at t =.... in the direction of independent variable x is much greater than the prescribed error tolerance.

This is my naive code.

ini = 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",(*"DiscretizedMonitorVariables"->True,*)
Method -> {"FixedStep",(*"StepSize"->0.0001,*)
Method -> {"ImplicitRungeKutta", "DifferenceOrder" -> 2,
"StepSizeControlParameters" -> {3/10,
0},(*"StepSizeSafetyFactors"->{8/10,9/10},*)
"ImplicitSolver" -> {"Newton", AccuracyGoal -> 8,
"IterationSafetyFactor" -> 1}}},
"SpatialDiscretization" -> {"TensorProductGrid",
"MinPoints" -> 180,
"MaxPoints" ->
180}},
AccuracyGoal -> 10, WorkingPrecision -> 32, StepMonitor :> Print[t]]


I have tried "ExplicitRungeKutta" and "Adams", and tried to adjust AccuracyGoal, but it does not work.

• Do you have reason to believe there is not a singularity at or near the stopping point? – Michael E2 Sep 15 '14 at 13:25

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];


NDSolveIterate::ndsz: At t == 4.0847545410125389609196036650468602942927835049548442055509298131500274581688440., 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.0847545410125389609196036650468602942927835049548442055509298131500274581688440. *)


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.

• @Michael.Thank you for your devoting to my problem. I'm so sorry for my late reply due to my cursoriness. Your answer is very helpful! Especially, you concluded The value of t where h is zero depends on the accuracy of computation. With time-derivative analysis combined time-integration along axis of symmetry, the suggestion AccuracyGoal -> Infinity is wonderful! Actually, I'm a MMA beginner. I do have some questions which could be easy to you:). 1. As you mentioned "FixedStep should be discarded", I used StepSizeControlParameters in my code, should it be varied time step method? – Enter Sep 21 '14 at 13:44
• Let me continue :P. 2. You employed "DifferenceOrder"-> "Pseudospectral" to evaluate the spatial derivative, do you think there are any other method would work better, such as, "DifferenceOrder" -> 2 or so. 3. Is the ** Newton's method** you mentioned the method of time integration, and how can I get the error estimate for time step from Newton's Method, as in your post -3.62148*10^-37 -- error estimate for t (Newton's method). 4. Please guide me to understand [[1, -10 ;;, 2]] // Differences, which you used to get the number of time steps and the last 9 step sizes. Thanks again! – Enter Sep 21 '14 at 14:07
• @Ixy, hi and thanks! :) I haven't used "StepSizeControlParameters". The description in the manual makes it seem like it might help. "DifferenceOrder"-> "Pseudospectral" seemed to produce smaller spatial error estimates. Each choice has advantages and disadvantages, which are discussed in the "The Numerical Method of Lines" tutorial. The Newton's Method estimate is a local estimate given by Newton's (tangent line) formula Δt = -f[t]/f'[t]. When the function is positive and concave down, this is a theoretical upper bound (but it ignores the accumulated global error of the integration)... – Michael E2 Sep 21 '14 at 14:32
• @Ixy:...In other words, from the computed solution, it is an upper limit on how far we can advance t. It shows that we have determined where the computed approximation is zero highly accurately; however, changing the spatial grid shows that the global error is significant. The notation [[1, -10 ;;, 2]] selects some of the values from the {x, t} grid used by NDSolve. The 1 gets the first line of integration; the -10 ;; code represents the last ten steps (look up Span); the 2 get the t part -- all in all, it is the last ten time steps. Differences gives the differences. – Michael E2 Sep 21 '14 at 14:42

I found another combination of option that doesn't suffer the NDSolve::eerr warning and produces almost the same result as that in Michael E2 but only takes about 0.6 s on my old laptop to finish. The key point here is to use a high enough "DifferenceOrder" and choose a odd number of grid points:

(* Define a auxiliary function *)
mol[n_, o_: "Pseudospectral"] := {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid", "MaxPoints" -> n,
"MinPoints" -> n, "DifferenceOrder" -> o}}

AbsoluteTiming[
sol = h /. First[NDSolve[{D[h[x, t], t] + D[D[h[x, t], x]/h[x, t],
x] + D[h[x, t]^3 D[h[x, t], {x, 3}], x] == 0,
h[(-(1/2)) (3 Sqrt[2] Pi), t] == h[(Sqrt[2] Pi)/2, t],
h[x, 0] == 1 + (1/10) Sin[x/Sqrt[2]]}, h,
{x, (-(1/2)) (3 Sqrt[2] Pi), (Sqrt[2] Pi)/2}, {t, 0, 6},
Method -> mol[75, 11]]]]

{{xl, xr}, {tl, tr}} = sol["Domain"];

Plot3D[sol[x, t], {x, xl, xr}, {t, tl, tr}, PlotRange -> All]
`