# How to demonstrate lack of stability with advection equation

I am trying to that using a coarse grid with an explicit method for, say, the advection equation leads to an unstable solution. The trouble is Mathematica avoids unstable solutions with good Methods or I am unable to set a coarse spatial grid or temporal grid spacing.

My code is as follows:

sol = NDSolve[{D[u[t, x], t] ==
0.5 D[u[t, x], x, x] + u[t, x] D[u[t, x], x],
u[t, -Pi] == u[t, Pi] == 0 , u[0, x] == Sin[x]},
u, {t, 0, 2}, {x, -Pi, Pi}, MaxStepFraction -> 1/10]


And a plot of sol is:

Plot3D[Evaluate[u[t, x] /. First[sol]], {t, 0, 2}, {x, -Pi, Pi},
PlotRange -> All


Should I not be using MaxStepFraction to demonstrate instability? What is the most coarse explicit Method I could use? I have tried Method->{"ExplicitEuler"} but this was good enough to resolve instabilities and I found a stable solution without kinks at larger times.

-

Forcing large enough constant step size (in time) for NDSolve method lead to instabilitity.

 sol = NDSolve[{D[u[t, x], t] == 0.5 D[u[t, x], x, x] + u[t, x] D[u[t, x], x],
u[t, -Pi] == u[t, Pi] == 0, u[0, x] == Sin[x]}, u, {t, 0, 2}, {x, -Pi, Pi},
Method -> {"FixedStep", "Method" -> "ExplicitEuler", "StepSize" -> 1/10}];

Plot3D[Evaluate[u[t, x] /. First[sol]], {t, 0, 2}, {x, -Pi, Pi},
PlotRange -> All]


Max step fraction you used "StartingStepSize" -> 1/10 gives maximal step size of order $1/20$ for this particular integration time interval, which is small enough and gives stable evolution for explicit Euler algorithm. Of course stability depends on the ratio of spatial and temporal grid spacings. For this set of options Mathematica sets (and adapts) spatial grid resolution automatically. To force specific spatial grid, with discrete differentiation method, use the "MethodOfLines" option

 sol = NDSolve[{D[u[t, x], t] == 0.5 D[u[t, x], x, x] + u[t, x] D[u[t, x], x],
u[t, -Pi] == u[t, Pi] == 0, u[0, x] == Sin[x]}, u, {t, 0, 2}, {x, -Pi, Pi},
Method -> {"MethodOfLines",
Method -> {"FixedStep", "Method" -> "ExplicitEuler", "StepSize" -> 1/4},
"SpatialDiscretization" -> {"TensorProductGrid", "MinPoints" -> 16,
"MaxPoints" -> 16, "DifferenceOrder" -> 2}}
]


With the same values of "MinPoints" and "MaxPoints" we force constant grid over the evolution. Manipulating this and the time step size you can demonstrate instability of the used method. For more on numerical methods for PDEs have a look at the tutorial.

-
Thank you. So is it difficult to back-out a "CFL" number out of this example when one uses mathematica since mma always "sets/adapts spatial grid resolution"? –  drN Aug 7 '13 at 17:12
@drN Yes of course, see the edit to my question. –  mmal Aug 7 '13 at 17:18
Why when I increase the MaxPoints while holding MinPoints constant, does the error increase? I would assume the reverse should happen... –  drN Aug 7 '13 at 17:30
@drN Yes it should, but it also depends on the steps size you set... –  mmal Aug 7 '13 at 18:02