Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

Is it possible to accurately solve the 1D Euler equations in Mathematica using NDSolve?

For example, let us consider the problem given here:

Using the notation $(r,v,e)$ for $(\rho,v,E)$ we can formulate the problem in Mathematica as:

eqs = {
    D[r[x, t], t] + D[r[x, t] v[x, t], x] == 0,
    D[r[x, t] v[x, t], t] + D[r[x, t] v[x, t]^2 + p[x, t], x] == 0,
    D[e[x, t], t] + D[v[x, t] (e[x, t] + p[x, t]), x] == 0

with $g:=\gamma$,

p[x_, t_] := (g - 1) (e[x, t] - r[x, t] v[x, t]^2/2);

and initial conditions

r0[x_] := 1.0 Boole[0 < x <= 0.5] + 0.25 Boole[0.5 < x <= 1.0];
v0[x_] := 0.0;
p0[x_] := 1.0 Boole[0 < x <= 0.5] + 0.1 Boole[0.5 < x <= 1.0];

Setting Dirichlet boundary conditions and throwing it into NDSolve:

ppR = 401;
ndsol = NDSolve[Join[eqs, {
    r[x, 0] == r0[x], r[0, t] == r0[0], r[1, t] == r0[1],
    v[x, 0] == v0[x], v[0, t] == v0[0], v[1, t] == v0[1],
    p[x, 0] == p0[x], p[0, t] == p0[0], p[1, t] == p0[1]}],
 {r, v, e}, {x, 0, 1}, {t, 0, 0.1},
 MaxSteps -> 10^3, PrecisionGoal -> 4, AccuracyGoal -> 4, 
 Method -> {"MethodOfLines", "Method" -> "StiffnessSwitching", 
 "SpatialDiscretization" -> {"TensorProductGrid", 
  "DifferenceOrder" -> Pseudospectral, "MaxPoints" -> {ppR}, 
  "MinPoints" -> {ppR}}}]

Taking the solution and plotting

DensityPlot[Evaluate[First[v[x, t] /. ndsol]], {x, 0, 1}, {t, 0, 0.1}, ColorFunction -> Hue]

Solution plot

We see some nasty numerical errors. Is there anyway to get an accurate solution using NDSolve?

Edit: This paper suggests that the Method of Lines should apply:

share|improve this question
Your link isn't working here – Dr. belisarius Oct 8 '12 at 22:24
Strange... Is it something I did wrong? – Dale Roberts Oct 8 '12 at 22:31
perhaps this answer helps? but blind experimentation is not very satisfactory – acl Oct 8 '12 at 22:55
I hammered on this for a while, and got nowhere. BTW, I think in the r0[x_]:= its 0.125 (least of your worries) If you go through and make everything a fraction e.g. 1.4 -> 14/10 then it doesn't even crank, crashing out with a 1/0 ComplexInfinity error... – Eric Brown Oct 9 '12 at 1:00
Both links work fine for me. Perhaps the server was down for maintenance. – stevenvh Oct 9 '12 at 7:29

This is a classical shock-tube problem in which a initially diaphragm separates a hi-pressure, high-density region from one of lower pressure and density. The classical exact solution has multiple discontinuities, a shock wave and a contact-surface (density discontinuity) that propagate to the right, and a continuous rarefaction wave traveling into the high-pressure gas.

The discontinuities cannot be captured by a pseudo spectral scheme,which produces the disastrous oscillations. This also occurs with many other spatial differencing schemes that have too little inherent dissipation. See the book "Difference Methods for Initial Value Problems" by Richtmeyer and Morton (1967), which discusses methods for adding "artificial viscosity" to the Euler equations. This enables the discontinuities to be successfully smeared out over Multiple grid points. There are many more modern methods that accomplish the same purpose. Try Googling the terms "Roe Scheme, ENO Scheme, Flux-Difference Splitting."

In short, Mathematica's NDSolve command is not yet able to handle directly the Euler equation flows that have discontinuous solutions, and this includes all mixed subsonic/supersonic flows such as your case. Theoretically, it should be able to handle subsonic flows, which are always continuous. Unfortunately, I have found to my dismay that NDSolve is not yet sophisticated enough to permit one to use the type of Characteristics-based inflow and outflow boundary conditions that are an absolute necessity for generating stable, accurate subsonic flow solutions. I posted a Question on this site two days ago, but the question was Closed Out (removed) by other more senior users who deemed it inappropriate. That question was titled "NDSolve rejects a common non-reflective subsonic outflow boundary condition for the 1-D Euler equations."

share|improve this answer
I have edited your closed question. You are raising correct issues here but you need to present the problem at hand in a right manner. Best of luck in future! – PlatoManiac May 2 '13 at 23:44

Your Answer


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.