# Unstable solution of 2D+1 time PDE with periodic boundary condition

Now I am trying to solve the following 2D+1 type of PDE:

$\partial_t u(t,x,y)=-y\partial_{x}u+\partial_{y}\left[a y+b sin(x)u+c\partial_{y}u\right]$

with $u(0,x,y)=\frac{1}{2\pi}e^{-((x-\pi/4)^2+y^2)/2}$

and periodic boundary condition:

$u(t,-\pi,y)=u(t,\pi,y)$

In $y$-direction, it is unbounded.

Here is the code:

a = 1;
T = 1;
ωcb = -50;
ωct = 50;
ωb = -5;
ωt = 5;
A = 10;
γ = 0.1;
kT = 0.1;
ufun = NDSolveValue[{D[u[t, θ, ω], t] == -D[ω u[t, θ, ω], θ] -
D[-A Sin[θ] u[t, θ, ω] - γ kT D[u[t, θ, ω], ω], ω] + 0.1 D[ω u[t, θ, ω], ω],
u[0, θ, ω] ==
1/(2 a π)
E^(-((θ - π/4)^2/(2 a^2)) - (ω)^2/(2 a^2)),
u[t, -π, ω] == u[t, π, ω],
u[t, θ, ωcb] == u[t, θ, ωct]},
u, {t, 0,
T}, {θ, -π, π}, {ω, ωcb, ωct},
Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid",
"MinPoints" -> 200, "MaxPoints" -> 1000}}];
plots = Table[
Plot3D[Abs[
ufun[t, θ, ω]], {θ, -π, π}, {ω, ωb, ωt}, PlotRange -> All,
ColorFunction -> "LakeColors"], {t, 0, T, .1}];
ListAnimate[plots]


The solution is a mess. I believe there are some problem on mesh? Also there are some message

NDSolveValue::mxsst: Using maximum number of grid points 1000 allowed by the MaxPoints or MinStepSize options for independent variable $\theta$.

NDSolveValue::ibcinc: Warning: boundary and initial conditions are inconsistent.

NDSolveValue::eerr: Warning: scaled local spatial error estimate of 62.663959713694915at t = 1. in the direction of independent variable θ is much greater than the prescribed error tolerance. Grid spacing with 1001 points may be too large to achieve the desired accuracy or precision. A singularity may have formed or a smaller grid spacing can be specified using the MaxStepSize or MinPoints method options.

But when I try

Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid",
"MinPoints" -> 200, "MaxPoints" -> 10000}}


It cannot manage to finish the program. What's going on ? Thanks for any suggestion!

• Is this equation taken from literature or derived by yourself? Can you add some background information? Aug 16, 2018 at 9:06
• Yes, I derive it by myself. It's the Klein-Kramers equation or 2D Fokker-Planck equation in rotational coordinate with sinusoidal potential. Aug 16, 2018 at 10:09
• The initial data is not symmetric, so the periodic boundary conditions for the variable θ are inconsistent with the initial data. Aug 16, 2018 at 11:25
• We must pay attention to the code or the equations that are declared by the author? The equations in the code and in the beginning of the topic do not coincide. Aug 16, 2018 at 12:00
• @AlexTrounev I guess OP's equation in the code is a variant of the one in the reference. (He mentioned that "I derive it by myself" in his first comment. ) Of course this is something that OP should clarify. Aug 16, 2018 at 12:12

You need the magic of "Pseudospectral":

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 == E^(-((θ - π/4)^2/(2 a^2)) - ω^2/(2 a^2))/(2 a π) /. 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, .04}];
ListAnimate[plots] # Update: A faster approach

Your equation turns out to be another example where the somewhat strange difference strategy discussed in this post causes trouble, so the problem can be solved With the fix function in that post:

(* The fix function isn't included in this post,
ufun = fix[T, 4]@
NDSolveValue[{eq, ic, u[t, -π, ω] == u[t, π, ω],
u[t, θ, ωcb] == u[t, θ, ωct]},
u, {t, 0, T}, {θ, -π, π}, {ω, ωcb, ωct},
Method -> mol]; // AbsoluteTiming
(* {13.415239, Null} *)


The resulting animation is similar so I'd like to omit it here.

• The method is good, but the problem is not correct, so we have a message on the output NDSolveValue::ibcinc: Warning: boundary and initial conditions are inconsistent. Aug 16, 2018 at 11:34
• @Alex ibcinc warning pops up because the i.c. isn't periodic technically, but it's close to constant near $\theta=\pm \pi$ i.e. approximately periodic, so the warning isn't a big deal, I think. One can eliminate the warning by e.g. modifying i.c. to With[{u = u[t, θ, ω]}, ic = u == Piecewise[{{E^(-((θ - π/4)^2/(2 a^2)) - ω^2/(2 a^2))/( 2 a π), Abs@θ < 0.9 Pi}}] /. t -> 0], and the solution won't change much. Of course when the b.c. involves derivative, ibcinc often means fatal error, see this post for more information. Aug 16, 2018 at 11:57
• Note that the equations in the code and in the beginning of the topic do not coincide. Aug 16, 2018 at 12:08
• Thanks for pointing out the problem. Indeed, they looks different but essentially the same. Aug 16, 2018 at 13:41
• Could someone explain what's going on? I am confused with the method. What's the problem with my original code? Anyway, thanks! Aug 16, 2018 at 13:43

I took the equations and the statement of the problem at the beginning of the topic and wrote my code. I got the expected rotation of the wave as a result.

{a, b, c, k, L, T} = {1, 10, 1, 1, 5, 2};
u0[x_, y_] := Exp[-(x^2 + y^2)/2]/(2*Pi)
eq = D[u[t, x, y], t] == -k*y*D[u[t, x, y], x] +
D[a*y + b*Sin[x]*u[t, x, y] + c*D[u[t, x, y], y], y];
ic = u[0, x, y] == u0[x, y];
bc = {u[t, x, -L] == u[t, x, L], u[t, -Pi, y] == u[t, Pi, y]};
sol = NDSolveValue[{eq, ic, bc},
u, {t, 0, T}, {x, -Pi, Pi}, {y, -L, L},
Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid",
"MinPoints" -> 40, "MaxPoints" -> 100,
"DifferenceOrder" -> "Pseudospectral"}}, MaxSteps -> 10^6]
plots = Table[
Plot3D[Abs[sol[t, x, y]], {x, -Pi, Pi}, {y, -L, L},
PlotRange -> All, ColorFunction -> "LakeColors",
Mesh -> None], {t, 0, T, .05*T}];
ListAnimate[plots] • The solution is perfect! Many thanks. Do you have any opinion on the mesh you use? I am not familiar with specifying these things. Aug 16, 2018 at 14:38
• I used a very simple option Mesh -> None Aug 16, 2018 at 14:46
• Sorry I must misunderstand. What I want to ask is Method -> {"MethodOfLines", "SpatialDiscretization" -> {"TensorProductGrid", "MinPoints" -> 40, "MaxPoints" -> 100, "DifferenceOrder" -> "Pseudospectral"}. How did you find those setup? Aug 16, 2018 at 14:50
• I solved a similar equation, only non-linear in the complex plane. I picked up the parameters experimentally. Aug 16, 2018 at 15:19
• Sorry.. I copy your code but I cannot reproduce your solution. If your ic is centered at origin, I guess it won't rotate, right? Aug 16, 2018 at 17:33