3
$\begingroup$

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:

enter image description here

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

$\endgroup$
4
$\begingroup$

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:

enter image description here

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

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}
*)

Mathematica graphics

| improve this answer | |
$\endgroup$
  • $\begingroup$ Nice! Very clear analysis, and educational! $\endgroup$ – John Joseph M. Carrasco Aug 8 '17 at 22:00

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.