# Instability, Courant Condition and Robustness about solving 2D+1 PDE

After several discussions, I would like to focus on the robustness of solving 2D+1 PDE by considering all suggested methods from @xzczd (see here) I found that the Ratio between the convection term and diffusion is crucial. Here is the code.

Clear["Global*"]
Clear[tosameorder, fix]
tosameorder[state_NDSolveStateData, order_] :=
state /. a_NDSolveFiniteDifferenceDerivativeFunction :>
NDSolveFiniteDifferenceDerivative[a@"DerivativeOrder",
a@"Coordinates", "DifferenceOrder" -> order,
PeriodicInterpolation -> a@"PeriodicInterpolation"]

fix[endtime_, order_] :=
Function[{ndsolve},
Module[{state =
First[NDSolveProcessEquations @@ Unevaluated@ndsolve],
newstate}, newstate = tosameorder[state, order];
NDSolveIterate[newstate, endtime];
Unevaluated[ndsolve][[2]] /. NDSolveProcessSolutions@newstate],
HoldAll]
a = 1;
T = 1;
ωb = -15; ωt = 15;
A = 6.5;
γ = .1;
kT = 0.1;
φ = 0;

With[{u = u[t,θ, ω]},
eq = D[u, t] == -D[ω u,θ] - D[-A Sin[3θ] u, ω] + γ (1 + Sin[3θ])  kT  D[
u, {ω, 2}] + γ  (1 + Sin[3θ]) D[ω u, ω];
ic = u == E^(-((ω^2 +(θ+Pi/4)^2)/(2 a^2))) 1/(2 π a) /. t -> 0];
ufun = NDSolveValue[{eq, ic, u[t, -π, ω] == u[t, π, ω],
u[t,θ, ωb] == 0, u[t,θ, ωt] == 0}, u, {t, 0, T}, {θ, -π, π}, {ω, ωb, ωt},
Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid",
"MaxPoints" -> {81, 51},
"MinPoints" -> {41, 31}}}]; // AbsoluteTiming
plots = Table[
Plot3D[Abs[ufun[t,θ, ω]], {θ, -π, π}, {ω, ωb, ωt}, AxesLabel -> Automatic,
PlotPoints -> 30, BoxRatios -> {Pi, ωb, 1},
ColorFunction -> "LakeColors", PlotRange -> All], {t, 0, T,
T/10}]; // AbsoluteTiming
ListAnimate[plots]


One can see that the coefficient of diffusion term (2nd-order term) is really small ($0.1*0.1*sin(3\theta)$) especially when $sin(3\theta)=-1$, how could it possible satisfy the Courant condition with such vanishing diffusion term?

The below is the result:

What I expect is the something similar to following (obtain by adding artificial diffusion)

The main question is the robustness about solving this partially convection-dominated problem in a efficient and stable way. Many thanks.

Max[γ[θ], 0.3] D[ u, {ω, 2}]


This is how I put the artificial diffusion. But for the problem, the angle-dependent diffusion and convection are different.

Full code is here, it's a little bit messy:

Clear["Global*"]
(*//////////////////////////////////////         parameters         \
/////////////////////////////////////////////////////////////////////////////////////////////////////////////\
*)

n = 3;
ϕ = π/2; vg0 = 5;
vg[t_] := vg[t] = vg0;
(*vg[t_]:=vg[t]=2vg0*1/(\[ExponentialE]^(k(t-Tp1))+1)-vg0;*)
(*vg[t_]:=vg[t]=2vg0*1/(\[ExponentialE]^(k(t-Tp1))+1)-2vg0*1/(\
\[ExponentialE]^(k(t-Tp2))+1)+vg0;*)

τi = 3; Γ = 10;

Vb = 0; μ[α_] := (-1)^(α + 1) Vb/2  ;

XTicks1 = Table[2 π*j, {j, -10, 10}];
XTicks2 = Table[π/6*j, {j, -10, 10}];
YTicks = Table[2 π*j, {j, -10, 10}];

(*//////////////////////////////////////         Karrasch poles and \
coefficients          \
/////////////////////////////////////////////////////////////////////////////////////////////////////////////\
*)

Np = 25; M = 2 Np;
B = Normal[
SparseArray[{Band[{2, 1}] ->
Table[N[1/(2 Sqrt[(2 n - 1) (2 n + 1)])], {n, 1, M - 1}],
Band[{1, 2}] ->
Table[N[1/(2 Sqrt[(2 n - 1) (2 n + 1)])], {n, 1, M - 1}]}, M]];
{bvals, bvecs} = Eigensystem[B];
Zp = Table[Abs[N[1/bvals[[2 p]]]], {p, 1, Np}];
Rp = Table[
N[(Normalize[bvecs[[2 p]]][[1]]/(2 bvals[[2 p]]))^2], {p, 1, Np}];

(*/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////\
*)

σ0[θ_, V_, τ0_, Γ_,
Vg_, ϕ_] := σ0[θ, V, τ0, Γ,
Vg, ϕ] =
1/2 - I/(
4 π) (PolyGamma[
1/2 - I/(
2 π) (τ0 (Sin[n θ + ϕ] -
Sin[n θ]) - V/2 + Vg) + Γ/(
4 π)] -
PolyGamma[
1/2 + I/(
2 π) (τ0 (Sin[n θ + ϕ] -
Sin[n θ]) - V/2 + Vg) + Γ/(
4 π)] +
PolyGamma[
1/2 - I/(
2 π) (τ0 (Sin[n θ + ϕ] -
Sin[n θ]) + V/2 + Vg) + Γ/(
4 π)] -
PolyGamma[
1/2 + I/(
2 π) (τ0 (Sin[n θ + ϕ] -
Sin[n θ]) + V/2 + Vg) + Γ/(
4 π)]);
τ[θ_, V_, τ0_, Γ_,
Vg_, ϕ_] := τ[θ, V, τ0, Γ,
Vg, ϕ] = -τ0 n (Cos[n θ + ϕ] -
Cos[n  θ]) σ0[θ,
V, τ0, Γ, Vg, ϕ];
σ1[θ_, V_, τ0_, Γ_,
Vg_, ϕ_] := σ1[θ, V, τ0, Γ,
Vg, ϕ] =
1/(8 π^2 Γ)*(PolyGamma[1,
1/2 - I/(
2 π) (τ0 (Sin[n θ + ϕ] -
Sin[n θ]) - V/2 + Vg) + Γ/(
4 π)] +
PolyGamma[1,
1/2 + I/(
2 π) (τ0 (Sin[n θ + ϕ] -
Sin[n θ]) - V/2 + Vg) + Γ/(
4 π)] +
PolyGamma[1,
1/2 - I/(
2 π) (τ0 (Sin[n θ + ϕ] -
Sin[n θ]) + V/2 + Vg) + Γ/(
4 π)] +
PolyGamma[1,
1/2 + I/(
2 π) (τ0 (Sin[n θ + ϕ] -
Sin[n θ]) + V/2 + Vg) + Γ/(
4 π)]) - \!$$\*UnderoverscriptBox[\(∑$$, $$p = 1$$, $$Np$$]$$Rp[\([$$$$p$$$$]$$] $$( \*FractionBox[\(3 \*SuperscriptBox[\((τ0 \((Sin[n\ θ + ϕ] - \ Sin[n\ θ])$$ - V/2 +
Vg)\), $$2$$] $$(Γ/2 + Zp[\([$$$$p$$$$]$$])\) -
\*SuperscriptBox[$$(Γ/2 + Zp[\([$$$$p$$$$]$$])\), $$3$$]\),
SuperscriptBox[$$( \*SuperscriptBox[\((\ τ0 \((Sin[n\ θ + ϕ] - \ Sin[n\ θ])$$ - V/2 + Vg)\), $$2$$] +
\*SuperscriptBox[$$(Γ/2 + Zp[\([$$$$p$$$$]$$])\), $$2$$])\), $$3$$]] +
\*FractionBox[$$3 \*SuperscriptBox[\((τ0 \((Sin[n\ θ + ϕ] - \ Sin[n\ θ])$$ + V/2 +
Vg)\), $$2$$] $$(Γ/2 + Zp[\([$$$$p$$$$]$$])\) -
\*SuperscriptBox[$$(Γ/2 + Zp[\([$$$$p$$$$]$$])\), $$3$$]\),
SuperscriptBox[$$( \*SuperscriptBox[\((\ τ0 \((Sin[n\ θ + ϕ] - \ Sin[n\ θ])$$ + V/2 + Vg)\), $$2$$] +
\*SuperscriptBox[$$(Γ/2 + Zp[\([$$$$p$$$$]$$])\), $$2$$])\), $$3$$]])\)\)\);
Uprime[θ_, τ0_, ϕ_] := τ0 n Cos[n θ];

mol[m_Integer, n_Integer] := {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid", "MaxPoints" -> m,
"MinPoints" -> n}}
τeff =
FunctionInterpolation[-τi n Cos[n θ] + τ[θ,
Vb, τi, Γ,
vg0, ϕ], {θ, -π, π}];
γ =
FunctionInterpolation[(τi n (Cos[n θ + ϕ] -
Cos[n  θ]))^2 σ1[θ,
Vb, τi, Γ,
vg0, ϕ], {θ, -π, π}, AccuracyGoal -> 5,
InterpolationPrecision -> MachinePrecision];
Plot[{τeff[θ], γ[θ]}, {θ, -π, \
π}, PlotRange -> All]
a = 0.7;
T = 20;
ωb = -8;
ωt = 8;
θ0 = -π/4;
L = 1;
With[{u = u[t, θ, ω]},
eq = D[u, t] == -ω D[ u, θ] -
1/L*τeff[θ] D[u, ω] +
Re[γ[θ]] D[ ω u, ω] +
Max[γ[θ], 0.3] D[ u, {ω, 2}];
ic = u ==
E^(-(ω^2 + (θ + π/
4)^2)/(2 a^2))/(2  π a^2) /. t -> 0];

ufun = NDSolveValue[{eq, ic,
u[t, -π, ω] == u[t, π, ω],
u[t, θ, ωb] == 0, u[t, θ, ωt] == 0},
u, {t, 0,
T}, {θ, -π, π}, {ω, ωb, ωt},
Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid",
"MaxPoints" -> {41, 61}, "MinPoints" -> {41, 61},
"DifferenceOrder" -> 4}}]; // AbsoluteTiming
(*Plot3D[ufun[T,θ,ω],{θ,-π,π},{ω,\
ωb,ωt},PlotRange\[Rule]All,AxesLabel\[Rule]Automatic,\
PlotPoints\[Rule]30,ColorFunction\[Rule]"LakeColors"]*)
plots = Table[
Plot3D[Abs[
ufun[t, θ, ω]], {θ, -π, π}, {\
ω, ωb, ωt}, PlotRange -> All,
AxesLabel -> Automatic, PlotPoints -> 30,
ColorFunction -> "LakeColors"], {t, 0, T,
T/50}]; // AbsoluteTiming
ListAnimate[plots] // AbsoluteTiming


Stable Julia code for the same problem: Matrix of 2D+1 PDE

function F_eff(x, Gamma, Delta, Q, EpOm, A)
return -A*sin(x + Delta)
end

# effective friction
function gamma_eff(x, Gamma, Delta, Q, EpOm, A)
return Gamma
end

# effective diffusion
function D_eff(x, Gamma, Delta, Q, EpOm, A, kBT)
return 2.0*kBT*Gamma
end

function make_FPE_matrix(xi,vj,Gamma, Delta, Q, EpOm, A, kBT)

Nx = length(xi)
Nv = length(vj)
hx = xi[2]-xi[1]
hv = vj[2]-vj[1]

mat = zeros(Float64, Nx*Nv+1, Nx*Nv)

for i=0:(Nx-1)
for j=0:(Nv-1)

mat[i*(Nv)+j+1,i*(Nv)+j+1] = -D_eff( xi[i+1], Gamma, Delta, Q, EpOm, A, kBT)/(4*hv^2)

# -d/dx (v P)
if i == 0
mat[i*(Nv)+j+1, (Nx-1)*(Nv)+j+1] = vj[j+1]/(2*hx) # PBC
end
if i > 0
mat[i*(Nv)+j+1, (i-1)*(Nv)+j+1] = vj[j+1]/(2*hx)
end
if i < Nx-1
mat[i*(Nv)+j+1, (i+1)*(Nv)+j+1] = -vj[j+1]/(2*hx)
end
if i == Nx-1
mat[i*(Nv)+j+1, (0)*(Nv)+j+1] = -vj[j+1]/(2*hx)   # PBC
end

if j > 0
mat[i*(Nv)+j+1, (i)*(Nv)+(j-1)+1] = F_eff(xi[i+1],Gamma, Delta, Q, EpOm, A)/(2*hv) -
gamma_eff(xi[i+1],Gamma, Delta, Q, EpOm, A)*vj[j+1-1]/(2*hv)
if j> 1
mat[i*(Nv)+j+1, (i)*(Nv)+(j-2) + 1] = 0.5*D_eff(xi[i+1], Gamma, Delta, Q, EpOm, A, kBT)/(4*hv^2)
end
end
if j < Nv-1
mat[i*(Nv)+j+1, (i)*(Nv)+(j+1)+1] = -F_eff(xi[i+1],Gamma, Delta, Q, EpOm, A)/(2*hv) +
gamma_eff(xi[i+1], Gamma, Delta, Q, EpOm, A)*vj[j+1+1]/(2*hv)
if j < Nv-2
mat[i*(Nv)+j+1, (i)*(Nv)+(j+2)+1] = 0.5*D_eff(xi[i+1], Gamma, Delta, Q, EpOm, A, kBT)/(4*hv^2)
end
end

end
end

for i = 1:Nx*Nv
mat[end,i] = hx*hv
end

return mat
end


Update (9/4) I try to run the compiled code (with VS2017 community). Here is my code:

Clear["Global*"]
(*//////////////////////////////////////         parameters         \
/////////////////////////////////////////////////////////////////////////////////////////////////////////////\
*)

n = 3;
ϕ = π/2; vg0 = 5;
vg[t_] := vg[t] = vg0;

τi = 3; Γ = 10;

Vb = 0; μ[α_] := (-1)^(α + 1) Vb/2  ;

XTicks1 = Table[2 π*j, {j, -10, 10}];
XTicks2 = Table[π/6*j, {j, -10, 10}];
YTicks = Table[2 π*j, {j, -10, 10}];

(*//////////////////////////////////////         Karrasch poles and \
coefficients          \
/////////////////////////////////////////////////////////////////////////////////////////////////////////////\
*)

Np = 25; M = 2 Np;
B = Normal[
SparseArray[{Band[{2, 1}] ->
Table[N[1/(2 Sqrt[(2 n - 1) (2 n + 1)])], {n, 1, M - 1}],
Band[{1, 2}] ->
Table[N[1/(2 Sqrt[(2 n - 1) (2 n + 1)])], {n, 1, M - 1}]}, M]];
{bvals, bvecs} = Eigensystem[B];
Zp = Table[Abs[N[1/bvals[[2 p]]]], {p, 1, Np}];
Rp = Table[
N[(Normalize[bvecs[[2 p]]][[1]]/(2 bvals[[2 p]]))^2], {p, 1, Np}];

(*/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////\
*)

σ0[θ_, V_, τ0_, Γ_,
Vg_, ϕ_] := σ0[θ, V, τ0, Γ,
Vg, ϕ] =
1/2 - I/(
4 π) (PolyGamma[
1/2 - I/(
2 π) (τ0 (Sin[n θ + ϕ] -
Sin[n θ]) - V/2 + Vg) + Γ/(
4 π)] -
PolyGamma[
1/2 + I/(
2 π) (τ0 (Sin[n θ + ϕ] -
Sin[n θ]) - V/2 + Vg) + Γ/(
4 π)] +
PolyGamma[
1/2 - I/(
2 π) (τ0 (Sin[n θ + ϕ] -
Sin[n θ]) + V/2 + Vg) + Γ/(
4 π)] -
PolyGamma[
1/2 + I/(
2 π) (τ0 (Sin[n θ + ϕ] -
Sin[n θ]) + V/2 + Vg) + Γ/(
4 π)]);
τ[θ_, V_, τ0_, Γ_,
Vg_, ϕ_] := τ[θ, V, τ0, Γ,
Vg, ϕ] = -τ0 n (Cos[n θ + ϕ] -
Cos[n  θ]) σ0[θ,
V, τ0, Γ, Vg, ϕ];
σ1[θ_, V_, τ0_, Γ_,
Vg_, ϕ_] := σ1[θ, V, τ0, Γ,
Vg, ϕ] =
1/(8 π^2 Γ)*(PolyGamma[1,
1/2 - I/(
2 π) (τ0 (Sin[n θ + ϕ] -
Sin[n θ]) - V/2 + Vg) + Γ/(
4 π)] +
PolyGamma[1,
1/2 + I/(
2 π) (τ0 (Sin[n θ + ϕ] -
Sin[n θ]) - V/2 + Vg) + Γ/(
4 π)] +
PolyGamma[1,
1/2 - I/(
2 π) (τ0 (Sin[n θ + ϕ] -
Sin[n θ]) + V/2 + Vg) + Γ/(
4 π)] +
PolyGamma[1,
1/2 + I/(
2 π) (τ0 (Sin[n θ + ϕ] -
Sin[n θ]) + V/2 + Vg) + Γ/(
4 π)]) - \!$$\*UnderoverscriptBox[\(∑$$, $$p = 1$$, $$Np$$]$$Rp[\([$$$$p$$$$]$$] $$( \*FractionBox[\(3 \*SuperscriptBox[\((τ0 \((Sin[n\ θ + ϕ] - \ Sin[n\ θ])$$ - V/2 +
Vg)\), $$2$$] $$(Γ/2 + Zp[\([$$$$p$$$$]$$])\) -
\*SuperscriptBox[$$(Γ/2 + Zp[\([$$$$p$$$$]$$])\), $$3$$]\),
SuperscriptBox[$$( \*SuperscriptBox[\((\ τ0 \((Sin[n\ θ + ϕ] - \ Sin[n\ θ])$$ - V/2 + Vg)\), $$2$$] +
\*SuperscriptBox[$$(Γ/2 + Zp[\([$$$$p$$$$]$$])\), $$2$$])\), $$3$$]] +
\*FractionBox[$$3 \*SuperscriptBox[\((τ0 \((Sin[n\ θ + ϕ] - \ Sin[n\ θ])$$ + V/2 +
Vg)\), $$2$$] $$(Γ/2 + Zp[\([$$$$p$$$$]$$])\) -
\*SuperscriptBox[$$(Γ/2 + Zp[\([$$$$p$$$$]$$])\), $$3$$]\),
SuperscriptBox[$$( \*SuperscriptBox[\((\ τ0 \((Sin[n\ θ + ϕ] - \ Sin[n\ θ])$$ + V/2 + Vg)\), $$2$$] +
\*SuperscriptBox[$$(Γ/2 + Zp[\([$$$$p$$$$]$$])\), $$2$$])\), $$3$$]])\)\)\);
Uprime[θ_, τ0_, ϕ_] := τ0 n Cos[n θ];

mol[m_Integer, n_Integer] := {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid", "MaxPoints" -> m,
"MinPoints" -> n}}
τeff =
FunctionInterpolation[-τi n Cos[n θ] + τ[θ,
Vb, τi, Γ,
vg0, ϕ], {θ, -π, π}];
γ =
FunctionInterpolation[(τi n (Cos[n θ + ϕ] -
Cos[n  θ]))^2 σ1[θ,
Vb, τi, Γ,
vg0, ϕ], {θ, -π, π}, AccuracyGoal -> 5,
InterpolationPrecision -> MachinePrecision];
Plot[{τeff[θ], γ[θ]}, {θ, -π, \
π}, PlotRange -> All]
a = 0.7;
T = 20;
ωb = -8;
ωt = 8;
θ0 = -π/4;
L = 1;
Clear[τeff, γ]
points@θ = 100; points@ω = 50;
delta@θ = (Pi + Pi)/(points@θ - 1);
delta@ω = (ωt - ωb)/(points@ω - 1);
τefflst =
Chop@N@Array[
Function[θ, -τi n Cos[n θ] + τ[θ, Vb, τi, Γ, vg0, ϕ]], points@θ, {-π, π}];
γlst =
Chop@Array[
Function[θ, (τi n (Cos[n θ + ϕ] - Cos[n θ]))^2 σ1[θ, Vb, τi, Γ, vg0, ϕ]],
points@θ, {-π, π}];

With[{u = u[θ, ω]},
rhs2 = -ω ct@D[u, θ] - 1/L τeff[θ] ct@D[u, ω] +
γ[θ] ct@D[ω u, ω] + γ[θ] fw@D[bw@D[u, ω], ω];
iclst2 = Table[E^(-((ω^2 + (θ + π/4)^2)/(2 a^2)))/(2 π a^2),
{θ, -Pi, Pi, delta@θ}, {ω, ωb, ωt, delta@ω}]];

rt = RescalingTransform[{{-Pi, Pi}, {ωb, ωt}}, {{1, points@θ}, {1, points@ω}}];
rttheta = RescalingTransform[{{-Pi, Pi}}, {{1, points@θ}}];

With[{rc = RuleCondition, cg = CompileGetElement},
rhsfunc2 = Hold@
Compile[{{u, _Real, 2}, {τeff, _Real, 1}, {γ, _Real, 1}},
Table[rhs2, {θ, -Pi, Pi, delta@θ}, {ω, ωb, ωt, delta@ω}],
RuntimeOptions -> EvaluateSymbolically -> False, CompilationTarget -> C] /.
OwnValues@rhs2 /.
u[theta_, omega_] :>
rc@(cg[u, Mod[#, points@θ - 1, 1], Mod[#2, points@ω - 1, 1]] & @@
Round@rt@{theta, omega}) /. (coef : τeff | γ)[theta_] :>
rc@(cg[coef, First@Round@rttheta@{theta}]) /. DownValues@delta /.
DownValues@points /. Flatten[OwnValues /@ Unevaluated@{ωb, ωt}] //
ReleaseHold // Last];
T = 30;
ulstfunc2 =
NDSolveValue[{u'[t] == rhsfunc2[u[t], τefflst, γlst], u[0] == iclst2},
u, {t, 0, T}, MaxSteps -> Infinity]; // AbsoluteTiming
(* {36.177260, Null} *)

lst = Table[
ListPlot3D[ulstfunc2[t]\[Transpose], PlotRange -> All,
DataRange -> {{-Pi, Pi}, {ωb, ωt}}], {t, 0, 5, 1/20}];
ListAnimate@lst


There is an error message:

CompiledFunction::rterr: -- Message text not found -- (compiledFunction5) (8)

• See section "Stabilization of Convection-Dominated Equations" in this guide. Note that this is also an ongoing research topic. Aug 29, 2018 at 9:17
• @HenrikSchumacher Thanks for the suggestion. Actually, the second animation is obtained from adding artificial diffusion as suggested in the documentation. In question is "How far away is solution relative to the actual solution" and "Can we actually improve the stability by specifying a proper grids or boundary condition". Of course, the Courant condition is an important issue but hard to control since time-step is adaptive, which is not able to modify manually. Aug 29, 2018 at 10:11
• "What I expect is the something similar to following (obtain by adding artificial diffusion)" How do you add the artificial diffusion? Can you include the specific setting? Aug 29, 2018 at 11:09
• @xzczd I put the code in my question. The idea is to keep some minimal value of diffusion such that the solution is still stable. Aug 29, 2018 at 13:08
• The Julia code is incomplete. Anyway, I managed to find a hint from it. Aug 30, 2018 at 11:26

OK, since neither of the default difference scheme nor the fix function works properly on this PDE, let's discretize the spatial derivative all by ourselves in a way different from that in NDSolve or in fix[…, …]@NDSolve.

The spatial discretization will be done in the following way: \begin{aligned} \frac{\partial^2 f}{\partial x^2}\Biggl|_{x=x_i} \approx & \frac{f(x_i+h)-2f(x_i)+f(x_i-h)}{h^2} \\ \frac{\partial f}{\partial x}\Biggl|_{x=x_i} \approx & \frac{f(x_i+h)-f(x_i-h)}{2h} \\ \frac{\partial (x f)}{\partial x}\Biggl|_{x=x_i} \approx & \frac{(x_i+h) f(x_i+h)-(x_i-h)f(x_i-h)}{2h} \end{aligned}

Notice the difference formula for $\frac{\partial (x f)}{\partial x}\Bigl|_{x=x_i}$ is critical here, because it produces a result different from

$$\frac{\partial (x f)}{\partial x}\Biggl|_{x=x_i}=\left(x\frac{\partial f}{\partial x}+f\frac{\partial x}{\partial x}\right)\Biggl|_{x=x_i}\approx x_i \frac{f(x_i+h)-f(x_i-h)}{2h}+f(x_i)$$

Periodic b.c. will be set in both directions because it's relatively easier to implement. Since the Dirichlet b.c. in $\omega$ direction is just an approximation for b.c. at infinity, this should not have any significant influence on the result.

The remaining work is just coding:

ωb = -5; ωt = 5;
a = 1; A = 6.5; γ = .1; kT = 0.1; φ = 0;

ClearAll[fw, bw, ct]
SetAttributes[#, HoldAll] & /@ {fw, bw, ct};
fw@D[expr_, x_] := Subtract @@ (expr /. {{x -> x + delta@x}, {x -> x}})/delta@x
bw@D[expr_, x_] := Subtract @@ (expr /. {{x -> x}, {x -> x - delta@x}})/delta@x
ct@D[expr_, x_] :=
Subtract @@ (expr /. {{x -> x + delta@x}, {x -> x - delta@x}})/(2 delta@x)

Clear[delta]
delta[a_ + b_] := delta@a + delta@b
delta[k_. delta[_]] := 0

points@θ = 100; points@ω = 50;
delta@θ = (Pi + Pi)/(points@θ - 1);
delta@ω = (ωt - ωb)/(points@ω - 1);
With[{u = u[θ, ω]},
rhs =-ct@D[ω u, θ] - ct@D[-A Sin[3 θ] u, ω] +
γ (1 + Sin[3 θ]) kT bw@D[fw@D[u, ω], ω] + γ (1 + Sin[3 θ]) ct@D[ω u, ω];
iclst = Table[
E^(-((ω^2 + (θ + Pi/4)^2)/(2. a^2))) 1/(2 π a),
{θ, -Pi, Pi, delta@θ}, {ω, ωb, ωt, delta@ω}]
];

rt = RescalingTransform[{{-Pi, Pi}, {ωb, ωt}}, {{1, points@θ}, {1, points@ω}}];

With[{rc = RuleCondition, cg = CompileGetElement},
rhsfunc = Hold@Compile[{{u, _Real, 2}},
Table[rhs, {θ, -Pi, Pi, delta@θ}, {ω, ωb, ωt, delta@ω}],
RuntimeOptions -> EvaluateSymbolically -> False,
CompilationTarget -> C] /. OwnValues@rhs /.
u[theta_, omega_] :>
rc@(cg[u, Mod[#, points@θ - 1, 1], Mod[#2, points@ω - 1, 1]] & @@
Round /@ rt@{theta, omega}) /. DownValues@delta /.
DownValues@points /. Flatten[OwnValues /@ Unevaluated@{ωb, ωt}] //
ReleaseHold // Last];

T = 30;
ulstfunc = NDSolveValue[{u'[t] == rhsfunc[u[t]], u[0] == iclst}, u, {t, 0, T},
MaxSteps -> Infinity]; // AbsoluteTiming
(* {33.1812583, Null} *)

lst = Table[
ListPlot3D[ulstfunc[t]\[Transpose], PlotRange -> All,
DataRange -> {{-Pi, Pi}, {ωb, ωt}}], {t, 0, T, 1}];
ListAnimate@lst


The result seems to be what you're expecting.

If you don't have a C compiler installed, take the CompilationTarget option and // Last away. I do recommend you to install one though.

Some advanced technique has been used in this code piece to make the discretization less tedious. To understand it better, you may want to read the following posts:

When should I, and when should I not, set the HoldAll attribute on a function I define?

# Update: Uncompilable coefficient function treatment

The coefficients of equation solved above is compilable. When they're not, such as in your 2nd code sample, some further modification is needed. The key idea is calculating coefficient values at grid points first and pass the value list to Compile:

Clear[τeff, γ]
points@θ = 100; points@ω = 50;
delta@θ = (Pi + Pi)/(points@θ - 1);
delta@ω = (ωt - ωb)/(points@ω - 1);
τefflst =
Chop@N@Array[
Function[θ, -τi n Cos[n θ] + τ[θ, Vb, τi, Γ, vg0, ϕ]], points@θ, {-π, π}];
γlst =
Chop@Array[
Function[θ, (τi n (Cos[n θ + ϕ] - Cos[n θ]))^2 σ1[θ, Vb, τi, Γ, vg0, ϕ]],
points@θ, {-π, π}];

With[{u = u[θ, ω]},
rhs2 = -ω ct@D[u, θ] - 1/L τeff[θ] ct@D[u, ω] +
γ[θ] ct@D[ω u, ω] + γ[θ] fw@D[bw@D[u, ω], ω];
iclst2 = Table[E^(-((ω^2 + (θ + π/4)^2)/(2 a^2)))/(2 π a^2),
{θ, -Pi, Pi, delta@θ}, {ω, ωb, ωt, delta@ω}]];

rt = RescalingTransform[{{-Pi, Pi}, {ωb, ωt}}, {{1, points@θ}, {1, points@ω}}];
rttheta = RescalingTransform[{{-Pi, Pi}}, {{1, points@θ}}];

With[{rc = RuleCondition, cg = CompileGetElement},
rhsfunc2 = Hold@
Compile[{{u, _Real, 2}, {τeff, _Real, 1}, {γ, _Real, 1}},
Table[rhs2, {θ, -Pi, Pi, delta@θ}, {ω, ωb, ωt, delta@ω}],
RuntimeOptions -> EvaluateSymbolically -> False, CompilationTarget -> C] /.
OwnValues@rhs2 /.
u[theta_, omega_] :>
rc@(cg[u, Mod[#, points@θ - 1, 1], Mod[#2, points@ω - 1, 1]] & @@
Round@rt@{theta, omega}) /. (coef : τeff | γ)[theta_] :>
rc@(cg[coef, First@Round@rttheta@{theta}]) /. DownValues@delta /.
DownValues@points /. Flatten[OwnValues /@ Unevaluated@{ωb, ωt}] //
ReleaseHold // Last];
T = 30;
ulstfunc2 =
NDSolveValue[{u'[t] == rhsfunc2[u[t], τefflst, γlst], u[0] == iclst2},
u, {t, 0, T}, MaxSteps -> Infinity]; // AbsoluteTiming
(* {36.177260, Null} *)

lst = Table[
ListPlot3D[ulstfunc2[t]\[Transpose], PlotRange -> All,
DataRange -> {{-Pi, Pi}, {ωb, ωt}}], {t, 0, 5, 1/20}];
ListAnimate@lst


Notice definitions of σ1 etc. are not included here, please find them in the body of the question.

• I notice that you use the same difference scheme in the d(xf)/dx as in my friend's Julia code. Such a subtle detail can really affect the stability very much? Anyway, thanks a lot! Aug 30, 2018 at 14:14
• @BobLin Yeah, surprising, but it's true. Notice this phenomenon is rather rare, before your question came out, this seems to be the only example in this site. Aug 30, 2018 at 14:28
• @BobLin If you're going to be doing a lot of numerical solution of PDEs for your work, you should get more general background on the topic. One nice starting place is Numerical Recipes by Press et al. Aug 30, 2018 at 21:11
• @xzczd Is it possible to apply your code to the PDE with f, which are interpolating function or user define function? In my second code, there are interpolating function like $τ_{eff}(θ)$ and $γ(θ)$. I want to apply your method to get the exact solution in the second example, but it seems to run quite too long. Aug 31, 2018 at 8:38
• I notice that the functions are always outside the the derivative. Something like With[{u = u[t, θ, ω]}, eq = D[u, t] == -ω D[ u, θ] - 1/L*τeff[θ] D[u, ω] + γ[θ] D[ ω u, ω] + γ[θ]D[ u, {ω, 2}]; ic = u == E^(-(ω^2 + (θ + π/ 4)^2)/(2 a^2))/(2 π a^2) /. t -> 0]; But it still very slow, when I run NDSolve. Aug 31, 2018 at 9:02