According to the documentation about the pseudospectral difference-order:

It says:

enter image description here

Following the discussion here:

I found the messy behavior is always on the artificial boundary in $\omega$-direction ($u(t,\theta,\omega_{cutoff})=0$ because I want $\omega$ to be unbounded.) Perhaps, this is so called Runge phenomenon? In principle, we should not use pseudospectral difference-order for all directions. However, it is not clear how to specify them separately.

Here is code:

a = 1;
T = 1;
ωb = -15; ωt = 15;
A = 8;
γ = .1;
kT = 0.1;
φ = 0;
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[3θ] u, ω] + γ (1 + Sin[3θ])  kT  D[
       u, {ω, 2}] + γ  (1 + Sin[3θ]) D[ω u, ω];
ic = u == E^(-((ω^2 +θ^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 -> mol[61], MaxSteps -> Infinity]; // 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


enter image description here


enter image description here


enter image description here

One can clearly see the large deviation occurs only in $\omega$-direction, which is consistent with the description as above (neither periodic nor Chebyshev).

Is it possible to have something like:

"DifferenceOrder" ->{"Pseudospectral", Automatic}

The above simply doesn't work.

Update: Finally, I figure out the problem is just due to convection-domination. The problem is depending on the ratio of convection term and diffusion term. Artificial diffusion or denser grid points is necessary.

Update(8/25): After using the implicit RungeKutta scheme, the solution is much stable. Now the another problem is the convergency.

What I expect is something similar to the following smooth behavior. enter image description here

But so far their is no such method which can arrive this, or? enter image description here

  • $\begingroup$ Perturbations increase at the borders of the region. It is necessary to experiment with the parameters, to determine what affects stability. I indicated several possible solutions, but there are others. And, of course, you can use other methods of solution. $\endgroup$ Aug 22, 2018 at 17:14
  • $\begingroup$ Your discussion includes, "large deviation occurs only in y-direction", but y is not one of the independent variables in your code. Also, the sentence, "In principle, we should not should pseudospectral difference-order for all direction." is garbled. For better responses by readers, please correct these and any other issues in the question. $\endgroup$
    – bbgodfrey
    Aug 22, 2018 at 18:08
  • 1
    $\begingroup$ Your understanding for the material is wrong. Notice the following sentence: However, there're two instances where this is not the case. Pseudospectral in Mathematica is implemented following these 2 instances. $\endgroup$
    – xzczd
    Aug 23, 2018 at 5:50
  • 1
    $\begingroup$ Once again, your understanding for the material is wrong. The meaning of that paragraph can be summarized as: 1. Pseudospectral derivative may cause Runge phenomenon; 2. But the implementation for Pseudospectral in NDSolve is free from Runge phenomenon. $\endgroup$
    – xzczd
    Aug 23, 2018 at 6:18
  • 1
    $\begingroup$ I didn't post the link as a possible solution… as far as I can tell, there's no universal technique for approximating b.c. at infinity, and I don't expect artificial b.c. designed for completely different PDE to work for your problem. Efficient artificial b.c. for certain PDE usually needs to be specially designed, and sadly that's beyond my reach. $\endgroup$
    – xzczd
    Aug 23, 2018 at 12:55

2 Answers 2


The computation in the question appears to suffer from a Courant instability. To illustrate, repeat the computation with higher plotting resolution and a slightly simpler code.

a = 1; T = 1; n = {61, 61};
ωb = -15; ωt = 15;
A = 8; γ = 1/10; kT = 1/10;
eq = D[u[t, θ, ω], t] == -D[ω u[t, θ, ω], θ] - D[-A Sin[3 θ] u[t, θ, ω], ω] 
  + γ (1 + Sin[3 θ]) kT  D[u[t, θ, ω], {ω, 2}] + γ (1 + Sin[3 θ]) D[ω u[t, θ, ω], ω];
ic = u[t, θ, ω] == E^(-((ω^2 + θ^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" -> n, "MinPoints" -> n, "DifferenceOrder" -> "Pseudospectral"}}]; 
Plot3D[Abs[ufunot[.9, θ, ω]], {θ, -π, π}, {ω, ωb, ωt}, PlotPoints -> 2 n, 
  PlotRange -> All, BoxRatios -> {Pi, ωb, 1}, ColorFunction -> "LakeColors", 
  ImageSize -> Large, AxesLabel -> {θ, ω, u}, LabelStyle -> {Black, Bold, Medium}, 
  Mesh -> None]

enter image description here

The higher plotting resolution displays significant fine structure in the numerical behavior near ω == ωt. A second plot focusing on that region makes the fine structure even more apparent.

Plot3D[Abs[ufun[T, θ, ω]], {θ, -π, π}, {ω, 12, ωt}, PlotPoints -> {122, 40}, 
  PlotRange -> All, ImageSize -> Large, AxesLabel -> {θ, ω, u}, 
  LabelStyle -> {Black, Bold, Medium}, Mesh -> None]

enter image description here

Spatial oscillations with wavelengths on the order of the cell size are the hallmark of the Courant instability. Reducing the number of grid points in ω from 61 to 59 to 57 steadily reduces the instability growth rate, and at 55 the instability disappears. Repeating the computation above with T = 10; n = {61, 55} shows no sign of the Courant instability.

enter image description here

There are, however, two obvious issues. First, waves have reached the boundaries in ω and appear to be reflecting. (The PDE is approximately advective at large Abs[ω].) Second, spatial resolution may no longer be sufficient to accurately represent the short wavelengths evident in the plot. The runtime for this computation was of order four minutes, and doubling the resolution would require over a half-hour of calculation, as well as some experimentation to find the optimal ratio of grid points in the two spatial dimensions. For completeness, here is a plot of the latter solution at T == 5, where the spatial resolution and boundary reflection problems are not yet significant.

enter image description here

  • $\begingroup$ Thanks. I agree with all you said since I actually did the same thing trying to optimize the grid points and boundary. However, this is really frustrating since It takes a lot of simulation time to obtain one experimental result but implies nothing useful still. What surprises me is the same thing can be done in Julia code by solving discretized matrix equation in about a minute. Would you think playing the computational method in time like explicit, implicit, Runge-Kutta methods would help? Again many thanks for your engagement to this tricky problem. $\endgroup$
    – Bob Lin
    Aug 24, 2018 at 19:56
  • $\begingroup$ Now we know this problem is really nontrivial, hope mathematica can solve this problem one day. $\endgroup$
    – Bob Lin
    Aug 24, 2018 at 20:30
  • 1
    $\begingroup$ @BobLin Mathematica can solve the problem today by the method you described in your comment above. Also, determining the largest eigenvalue of the matrix then would give the Courant limit. Expanding the PDE as a Fourier series in θ might also be interesting., decomposing the PDE into a coupled system of ODEs. Obtaining the optimal boundary condition in ω is tricky but might be feasible, starting from the Fourier decomposition method I just mentioned and then estimating analytically the boundary condition that minimizes reflections. In all, this is a very interesting problem. $\endgroup$
    – bbgodfrey
    Aug 24, 2018 at 21:49
  • $\begingroup$ @BobLin As I asked before, what difference scheme is used in that Julia code? $\endgroup$
    – xzczd
    Aug 25, 2018 at 2:14
  • 1
    $\begingroup$ @BobLin No, the time integration is trivial, and the difference scheme plays a critical role in this problem. I've included the solution here. With this solution the problem in this question can be solved in less than 10 seconds in NDSolve with Method -> {"MethodOfLines", "SpatialDiscretization" -> {"TensorProductGrid", "MaxPoints" -> 50, "MinPoints" -> 50, "DifferenceOrder" -> 2}}. $\endgroup$
    – xzczd
    Aug 26, 2018 at 12:53

You might try some different NDSolve methods. E.g.,

mol[n_Integer, o_: "Pseudospectral"] := {"MethodOfLines", "SpatialDiscretization" ->
  {"TensorProductGrid", "MaxPoints" -> n, "MinPoints" -> n, "DifferenceOrder" -> o},
  Method -> {"IDA", "ImplicitSolver" -> {"GMRES"}}}

takes twice as much time as your code, but doesn't have the artifacts at the boundary at t==0.9:

Mathematica graphics

  • 1
    $\begingroup$ Just a side note: this is probably related to 1. DifferentiateBoundaryConditions -> False. When IDA is chosen, the PDE is discretized to a DAE system and the b.c. is followed in a more strict way. The main disadvantage of this method is, we know, DAE solver of NDSolve is weaker than its ODE solver, and in certain cases it even fails for solvable problem. Related: mathematica.stackexchange.com/a/127411/1871 2. The implicit solver, because implicit method is usually unconditionally stable or almost unconditionally stable so is free from the issue discussed in bbgodfrey's answer. $\endgroup$
    – xzczd
    Aug 25, 2018 at 2:21
  • $\begingroup$ For what it's worth, I've always had much faster results solving PDEs by implementing the method of lines using VODE in Fortran, but that's much uglier than Mathematica. $\endgroup$
    – Chris K
    Aug 25, 2018 at 2:25
  • $\begingroup$ Just a side note 2: But one of the disadvantage of implicit solver is, it's usually slower than the explicit solver. $\endgroup$
    – xzczd
    Aug 25, 2018 at 2:37
  • $\begingroup$ @ChrisK Nice approach (+1). I ran your calculation to T = 10. Your solution was a bit noisier at the upper boundary than my solution but also a bit faster computationally. $\endgroup$
    – bbgodfrey
    Aug 25, 2018 at 5:24

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.