Is it possible to accurately solve the 1D Euler equations in Mathematica using NDSolve
?
For example, let us consider the Sod shock tube problem. Introduction to this problem can be found in the following links:
http://www.csun.edu/~jb715473/examples/euler1d.htm
https://en.wikipedia.org/wiki/Sod_shock_tube
Using the notation {r, v, e}
for $(\rho,v,E)$ we can formulate the problem in Mathematica as:
g = 1.4;
p[x_, t_] := (g - 1) (e[x, t] - r[x, t] v[x, t]^2/2);
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
};
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]
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: http://math.lanl.gov/~mac/papers/numerics/H79.pdf