# Solution distorsion after long time evolution for 2D+1 PDE [closed]

Continue the discussion 2D+1 PDE problem.

Why the solution eventually becomes messy and show those warning messages.

Here is the code:

a = 1;
T = 10;
ωcb = -50;
ωct = 50;
ωb = -5;
ωt = 5;
A = 10;
γ = 0.1;
kT = 0.1;
φ=π/4;
mol[n_Integer, o_:"Pseudospectral"] := {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid", "MaxPoints" -> n,
"MinPoints" -> n, "DifferenceOrder" -> o}}

With[{u = u[t, θ, ω]},
eq = D[u, t] == -D[ω u, θ] - D[-A Sin[θ] u, ω] - γ kT D[u, ω] + 1/10 D[ω u, ω];

ic = u == EllipticTheta[3, (θ - φ)/2, E^(-a^2/2)]*
E^(-ω^2/(2 a^2))/(2  π)^(3/2) /. t -> 0];

ufun = NDSolveValue[{eq, ic, u[t, -π, ω] == u[t, π, ω], u[t, θ, ωcb] == u[t, θ, ωct]},
u, {t, 0, T}, {θ, -π, π}, {ω, ωcb, ωct},
Method -> mol]; // AbsoluteTiming
(* {24.137131, Null} *)

plots = Table[
Plot3D[Abs[ufun[t, θ, ω]], {θ, -π, π}, {ω, ωb, ωt},
PlotRange -> All, AxesLabel -> Automatic, PlotPoints -> 50,
BoxRatios -> {Pi, ωb, 1}], {t, 0, T, .2}];
ListAnimate[plots] It gives the following message:

NDSolveValue::mxst: Maximum number of 10000 steps reached at the point t == 6.689127286306131.
InterpolatingFunction::dmval: Input value {6.8,-3.14146,-4.9998} lies outside the range of data in the interpolating function. Extrapolation will be used.
InterpolatingFunction::dmval: Input value {7.,-3.14146,-4.9998} lies outside the range of data in the interpolating function. Extrapolation will be used.
InterpolatingFunction::dmval: Input value {7.2,-3.14146,-4.9998} lies outside the range of data in the interpolating function. Extrapolation will be used.
General::stop: Further output of InterpolatingFunction::dmval will be suppressed during this calculation.


How can this happen? I try to enlarge the MinPoints/MaxPoints, but it makes it worse. Any obvious mistake?

Would it related to the Courant condition? If the time step is too large, is there any thumb rule to better control it?

• Your code doesn't run when pasted into a fresh session. Could you check what's missing? Anyhow, the messages give big clues to the problem. The first one says NDSolve gave up at t==6.689 because MaxSteps was reached. The other messages indicate that evaluating the interpolating function after that point is based on extrapolation, hence the weird results. Possible simple solution: add MaxSteps->Infinity to your NDSolve. – Chris K Aug 17 '18 at 17:30
• Thanks! I missed a parameter φ=π/4; – Bob Lin Aug 17 '18 at 17:34
• I am wondering why mathematica can't treat this 2D+1 problem so well just like other problems. – Bob Lin Aug 17 '18 at 17:41
• Numerically solving PDEs is hard. There are entire books, courses, and careers devoted to the subject. Having done it manually, I'm amazed that Mathematica can do it at all. Anyhow, it will sometimes need some guidance from the user as to appropriate numerical methods. – Chris K Aug 17 '18 at 18:04

Adding MaxSteps->Infinity to NDSolveValue fixes the problem.

Plot3D[Abs[ufun[T, θ, ω]], {θ, -π, π}, {ω, ωb, ωt}, PlotRange -> All,
AxesLabel -> Automatic, PlotPoints -> 50, BoxRatios -> {Pi, ωb, 1}] • Sorry there is a typo in my code. The function EllipticTheta[3, (0 - φ)/2, E^(-a^2/2)] should be instead EllipticTheta[3, (θ - φ)/2, E^(-a^2/2)]. But your suggestion is perfect. – Bob Lin Aug 17 '18 at 18:24
• One question about the MaxStep. Does that mean my computation steps has already exceeded 10000? – Bob Lin Aug 17 '18 at 18:26
• Yes, that's what the NDSolveValue::mxst is saying. – Chris K Aug 17 '18 at 18:29
• And the program stop running in my original code, or? – Bob Lin Aug 17 '18 at 18:30
• Yes, NDSolveValue stops at that point. Look at the domain of ufun and you'll see it goes only to t==6.689. When you Plot3D at larger times, the warning InterpolatingFunction::dmval tells you that the t` value is outside the domain, so it can't be trusted. – Chris K Aug 17 '18 at 18:35