# What is going wrong in this simulation using NDSolve

I am trying to simulate the non-linear KdV using NDSolve with initial condition -2Sech^2(x) and periodic boundary conditions. Here is my code:

a = 30;

eq = {D[u[t, x], {t, 1}] + D[u[t, x], {x, 3}] ==
6 u[t, x] D[u[t, x], {x, 1}], u[0, x] == -2 Sech[x] ^2,
u[t, -a] == u[t, a]};

sol = NDSolve[eq, u, {t, 0, 30}, {x, -a, a}, MaxStepSize -> 0.14];

data = Flatten[Table[{t, x, u[t, x]} /. sol, {t, 0, 30}, {x, -a, a}],
2];

ListPlot3D[data, Mesh -> None, ColorFunction -> "Rainbow",
PlotRange -> {{0, 30}, {-30, 30}, {-2, 2}}, Lighting -> "Automatic",
AxesLabel -> {"t", "x"}]


The output is correct at first but then something unexpected happens:

What is going wrong here? Normally the peaked-shaped object would continue to propagate.

Looks like it's just a glitch with too big a MaxStepSize:

sol = NDSolve[eq, u, {t, 0, 30}, {x, -a, a}, MaxStepSize -> .07];

yields:

• I compared this method with some alternatives in my answer, in case you're interested. (+1) – Michael E2 Aug 8 '17 at 21:49

Another approach is to pick a better method for approximating the derivatives. Setting "DifferenceOrder" to either 6 or "Pseudospectral" solves the problem.

sol6 = NDSolve[eq, u, {t, 0, 30}, {x, -a, a},
Method -> {"MethodOfLines",
"SpatialDiscretization" -> {
"TensorProductGrid", "DifferenceOrder" -> 6},
"TemporalVariable" -> t}]; // AbsoluteTiming
(*  {2.66469, Null}  *)

solSP = NDSolve[eq, u, {t, 0, 30}, {x, -a, a},
Method -> {"MethodOfLines",
"SpatialDiscretization" -> {
"TensorProductGrid", "DifferenceOrder" -> "Pseudospectral"},
"TemporalVariable" -> t}]; // AbsoluteTiming
(*  {3.22224, Null}  *)


Note that the spatial discretization corresponding to "MaxStepSize" -> 0.07 is equivalent to "MinPoints" -> 858. However, the option MaxStepSize -> 0.07 used by itself also limits the time steps.

sol07s = NDSolve[eq, u, {t, 0, 30}, {x, -a, a},
Method -> {"MethodOfLines",
"SpatialDiscretization" -> {
"TensorProductGrid", "MaxStepSize" -> 0.07},
"TemporalVariable" -> t}]; // AbsoluteTiming
(*  {3.95492, Null}  *)

sol858 = NDSolve[eq, u, {t, 0, 30}, {x, -a, a},
Method -> {"MethodOfLines",
"SpatialDiscretization" -> {
"TensorProductGrid", "MinPoints" -> 858},
"TemporalVariable" -> t}]; // AbsoluteTiming
(*  {4.10221, Null}  *)

sol07 = NDSolve[eq, u, {t, 0, 30}, {x, -a, a},  (* John Joseph M. Carrasco *)
MaxStepSize -> .07]; // AbsoluteTiming
(*  {6.75514, Null}  *)


The two with only a limitation on spatial discretization are the same:

sol858 === sol07s
(*  True  *)


But the last way sol07, which limits both the spatial and time step size, is different:

Length /@ First[u["Coordinates"] /. sol858]
Length /@ First[u["Coordinates"] /. sol07]
(*
{173, 859}
{432, 859}
*)


By comparison, the "DifferenceOrder" adjustment yields a smaller solution (with respect to the number of steps and memory):

Length /@ First[u["Coordinates"] /. sol6]
Length /@ First[u["Coordinates"] /. solSP]
(*
{173, 491}
{175, 239}
*)


• Nice! Very clear analysis, and educational! – John Joseph M. Carrasco Aug 8 '17 at 22:00