1

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[50]]; // 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]

Time evolution

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?

Bob Lin
  • 445
  • 2
  • 8
  • 2
    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
  • 1
    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

1 Answers1

4

Adding MaxSteps->Infinity to NDSolveValue fixes the problem.

Plot3D[Abs[ufun[T, θ, ω]], {θ, -π, π}, {ω, ωb, ωt}, PlotRange -> All,
  AxesLabel -> Automatic, PlotPoints -> 50, BoxRatios -> {Pi, ωb, 1}]

Mathematica graphics

Chris K
  • 20,207
  • 3
  • 39
  • 74
  • 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
  • I am still confused why as the time evolves, the computational effort will also increase and eventually kill the program. Is this a normal issue in PDE? – Bob Lin Aug 17 '18 at 18:33
  • 1
    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
  • If you want more solution, NDSolve will generally take more steps. Just put MaxSteps->Infinity whenever you run into this problem. – Chris K Aug 17 '18 at 18:36
  • Could you replace the EllipticTheta[3, (0 - φ)/2, E^(-a^2/2)] to EllipticTheta[3, (θ - φ)/2, E^(-a^2/2)] , will you manage to obtain the solution? – Bob Lin Aug 17 '18 at 18:37
  • I also found that the solution's amplitude is increasing. That's strange. Hmm... – Bob Lin Aug 17 '18 at 18:43