Solving a PDE gives stiffness error

I'm solving a PDE with

Needs["DifferentialEquationsNDSolveProblems"];
Needs["DifferentialEquationsNDSolveUtilities"];

tmin = -2; tmax = 2;
(* u[x, 0] \[Equal] 4Sech[x]^2*)

ν = 1; ϵ = 0.001;(*0.005*);

xmin = -2; xmax = 2;
sol =
NDSolve[ {D[u[x, t], t] +
u[x, t]*D[u[x, t], x] + ϵ^2*D[u[x, t], {x, 3}] == ν*
D[u[x, t], {x, 2}], u[x, 0] == Cos[x], u[xmin, t] ==  u[xmax, t] },
u, {x, xmin, xmax}, {t, tmin, tmax}, Method -> "Automatic"]


During evaluation of In[40]:= NDSolve::mxsst: Using maximum number of grid points 10000 allowed by the MaxPoints or MinStepSize options for independent variable x.

During evaluation of In[40]:= NDSolve::ndsz: At t == -5.72262*10^-7, step size is effectively zero; singularity or stiff system suspected.

During evaluation of In[40]:= NDSolve::eerr: Warning: scaled local spatial error estimate of 5.058814843192396^10 at t = -5.7226210^-7 in the direction of independent variable x is much greater than the prescribed error tolerance. Grid spacing with 10001 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.

(* Out[40]= {{u -> -7 InterpolatingFunction[{{…, -2., 2., …}, {-5.72262 10 , 2.}}, <>]}} *)

Solve Burgers in form ν ∂u^2/∂^2x= -(∂u/∂t)- u ∂u/∂x-ϵ^2 ∂^3u/∂x^3

Plot3D[u[x, t] /. Flatten[sol], {x, xmin, xmax}, {t, tmin, tmax},
PlotPoints -> 100, PlotRange -> All, ColorFunction -> "Rainbow"]


But can't get an output because

NDSolve::ndsz: At t == -5.72262*10^-7, step size is effectively zero; singularity or stiff system suspected.

Any ideas how to fix?

• There are more than 100 questions on the topic on this site. Have you gone through those? Have you tried switching to a different Method such as Method -> "StiffnessSwitching"? Commented Apr 8, 2022 at 14:44
• Yes, tried that and Runge Kutta. Still didn't resolve it Commented Apr 8, 2022 at 14:46
• With respect to the question title, the message does not just say stiffness. The first problem suggested is a singularity. Did you check that? The NDSolve::mxsst warning tends to make that more plausible. Commented Apr 8, 2022 at 22:36
• Perhaps noteworthy: the initial condition is u[x, 0] == Cos[x] at t=0 yet you want to NDSolve from t=-2. If u[x, -2] == Cos[x] or the NDSolve is from t=0 to 2 Mathematica can solve it. Is this a typo?
– bmf
Commented Apr 8, 2022 at 22:46
• Yep, tyop. How would this new solution look? Commented Apr 9, 2022 at 16:48

Note: in the comment section under the OP we found a typo, which has been taken into consideration in the suggested solution below.

The following runs without any errors and complaints

tmin = -2;
tmax = 2;
ν = 1;
ε = 0.001;
xmin = -2;
xmax = 2;
points = 2000;

sol = NDSolveValue[{D[u[x, t], t] +
u[x, t]*D[u[x, t], x] + ε^2*D[u[x, t], {x, 3}] == ν*
D[u[x, t], {x, 2}], u[x, -2] == Cos[x],
u[xmin, t] == u[xmax, t]},
u[x, t], {x, xmin, xmax}, {t, tmin, tmax},
Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid",
"MaxPoints" -> points, "MinPoints" -> points,
"DifferenceOrder" -> "Pseudospectral"}, Method -> "LSODA"},
MaxSteps -> Infinity]


Visualizing

Plot3D[sol, {x, xmin, xmax}, {t, tmin, tmax}]
`

• ... which was the problem I stated with Commented Apr 8, 2022 at 15:16
• @Roberto_1986 do you agree with the newly suggested solution that has no errors?
– bmf
Commented Apr 26, 2022 at 13:59