Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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

enter image description here

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.

share|improve this question

1 Answer 1

up vote 3 down vote accepted

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]

enter image description here

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.

share|improve this answer
    
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

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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