The fluid here has been assumed as single component perfect gas i.e. it obeys the equation $p=ρ R T$, the thermal conductivity is assumed as a constant, so the equation set is:

NDSolve[{D[ρg[t, x], t] + D[ρg[t, x] u[t, x], x] == 0,

         ρg[t, x] D[u[t, x], t] + ρg[t, x] u[t, x] D[u[t, x],x] 
          == -D[ρg[t, x] te[t, x], x],

         ρg[t, x] D[te[t, x], t] + ρg[t, x] u[t, x] D[te[t, x],x]
          == -ρg[t, x] te[t, x] D[u[t, x], x] + D[te[t, x], x, x],

         te[0, x] == 298, te[t, -0.5] == 298, te[t, 0.5] == 298,

         ρg[0, x] == (1 - x^2),

         u[0, x] == 0, u[t, -0.5] == 0, u[t, 0.5] == 0},

        {ρg[t, x], te[t, x], u[t, x]}, {t, 0, 1}, {x, -0.5, 0.5}]

After I ran the code, I only got the warning message NDSolve::ndsz and NDSolve::eerr, I've checked the the equations for times and I think they are correct, and the initial and boundary conditions are also simple and seems to be reasonable (at least from the perspective of physics). So…What's wrong with it?…Well, to tell you the truth, what I really want to ask is, does NDSolve lack the ability to solve system of partial differential equations?

Oh, someone may feel strange that there's no boundary condition for ρg[t, x], that's because, I found that only four boundary conditions are necessary for the solving of the equations though I don't know the exact reason (I found it in times of trial when I set ρg[0, x] as a constant ).

  • 1
    $\begingroup$ Did you check carefully all signs for all terms? Can this be just a typo with sign leading to sort of exp growth? $\endgroup$ Aug 10, 2012 at 5:31
  • $\begingroup$ @Vitaliy Kaurov After I saw your comment I checked all signs for another time, and…I still don't find any mistake… $\endgroup$
    – xzczd
    Aug 10, 2012 at 5:47
  • $\begingroup$ @Vitaliy Kaurov Er…what about this? See the 5.17, 5,31, 5.77. (Of course we still need some deformation to get the form for the one-dimensional cases.) $\endgroup$
    – xzczd
    Aug 10, 2012 at 6:17
  • 1
    $\begingroup$ 5.77 seems to differ quite a bit from your implementation. Sure you're correct there? $\endgroup$ Aug 10, 2012 at 6:37
  • $\begingroup$ @SjoerdC.deVries I think you mean the expression involved τ (shear stress), right? Since the fluid here is ideal gas (oh, maybe I should say perfect gas for it's in the field of Fluid Mechanics), it's zero. $\endgroup$
    – xzczd
    Aug 10, 2012 at 6:47

2 Answers 2


Perhaps setting the difference order to "DifferenceOrder" -> "Pseudospectral" is what you are looking for:

showStatus[status_] := 
clearStatus[] := showStatus[""];
nxy = 33;
sol = NDSolve[{D[ρg[t, x], t] + D[ρg[t, x] u[t, x], x] == 
   0, ρg[t, x] D[u[t, x], t] + ρg[t, x] u[t, x] D[u[t, x], 
      x] == -D[ρg[t, x] te[t, x], 
     x], ρg[t, x] D[te[t, x], t] + ρg[t, x] u[t, x] D[
      te[t, x], x] == -ρg[t, x] te[t, x] D[u[t, x], x] + 
    D[te[t, x], x, x], te[0, x] == 298, te[t, -0.5] == 298, 
  te[t, 0.5] == 298, ρg[0, x] == (1 - x^2), u[0, x] == 0, 
  u[t, -0.5] == 0, u[t, 0.5] == 0}, {ρg[t, x], te[t, x], 
  u[t, x]}, {t, 0, 1}, {x, -0.5, 0.5},
 Method -> {"MethodOfLines", 
   "SpatialDiscretization" -> {"TensorProductGrid", 
     "MaxPoints" -> nxy, "MinPoints" -> nxy, 
     "DifferenceOrder" -> "Pseudospectral"}, Method -> "Adams"},
 MaxSteps -> Infinity, 
 EvaluationMonitor :> showStatus["t = " <> ToString[CForm[t]]]];

The Method->Adams is not necessary.

Depicted below is how the solutions look. To avoid scales incongruousness they are visualized on different plots.

Plot3D[Evaluate[{#[t, x]} /. sol], {t, 0, .2}, {x, -0.5, 0.5}, 
    Mesh -> True, MeshStyle -> Opacity[.2], 
    ColorFunction -> "DarkRainbow", PlotStyle -> Opacity[.7], 
    PlotRange -> All, PlotPoints -> 40, ImageSize -> 200, 
    PlotLabel -> #] & /@ {ρg, u, te} // Row

enter image description here

  • $\begingroup$ any chance of making it faster? Is it the optimal speed for NDSolve? $\endgroup$ Aug 10, 2012 at 19:26
  • $\begingroup$ @ruebenko Big +1 - I am in awe you got this to work. I was almost sure there is an intrinsic divergence. I was also trying "Pseudospectral" with cyclic BC but to no avail. I am still puzzled why this works an nothing else does. I took a liberty to plot the beautiful solutions. Please feel free to remove my edit if you feel it's inappropriate. $\endgroup$ Aug 11, 2012 at 9:15
  • 5
    $\begingroup$ @ruebenko: can you explain how you found that exactly this combination of options does work? Do you have a strategy how to find a working set of parameters or are you using a clever combination of experience+guessing+trying? If you could explain, that would probably make the answer even more valuable, you eared a +1 anyway... $\endgroup$ Aug 11, 2012 at 10:15
  • 4
    $\begingroup$ @AlbertRetey, in this case the error message NDSolve::ndsz made me try Method->Adams (which actually is not needed) as an alternative to the automatic stiffness switching. NDSolve::eerr: made me choose my own grid. The DifferenceOrder was because OP mentioned something about the boundary conditions in the post and that made me think this might be worth a try. So, yes, some experience and a try. I do have a note where I collect these type of questions and solutions that are found for them and then I consult this when needed. This might be something for a community wiki. $\endgroup$
    – user21
    Aug 13, 2012 at 16:19
  • 2
    $\begingroup$ As late as this is, changing the Method from Adams to LSODA made is really fast to run. $\endgroup$
    – dearN
    Feb 16, 2018 at 17:29

With the experience gained in past 6 years, I manage to find out a more efficient solution for this problem :D .

This problem turns out to be another example on which the difference scheme implemented in NDSolve doesn't work well. The issue has been discussed in this post. Equipped with the fix function therein, the ndsz warning no longer pops up. Though the eerr warning remains, the estimated error is small, and the result is consistent with the one in ruebenko (now user21)'s answer.

xL = -1/2; xR = 1/2;

With[{T = T[t, x], ρ = ρ[t, x], u = u[t, x]},

 eq = {D[ρ, t] + D[ρ u, x] == 0,
       D[u, t] + u D[u, x] == -(1/ρ) D[ρ T, x],
       D[T, t] + u D[T, x] == -T D[u, x] + 1/ρ D[T, x, x]};

 ic = {T == 298, ρ == (1 - x^2), u == 0} /. t -> 0;

 bc = {{T == 298, u == 0} /. x -> xL, {T == 298, u == 0} /. x -> xR};]

endtime = 0.2; difforder = 2; points = 300;

mol[n : _Integer | {_Integer ..}, o_: "Pseudospectral"] := {"MethodOfLines", 
    "SpatialDiscretization" -> {"TensorProductGrid", "MaxPoints" -> n, 
        "MinPoints" -> n, "DifferenceOrder" -> o}}
mol[tf : False | True, sf_: Automatic] := {"MethodOfLines",
  "DifferentiateBoundaryConditions" -> {tf, "ScaleFactor" -> sf}}

(* Definition of fix isn't included in this post, 
   please find it in the link above. *)
mysol = fix[endtime, difforder]@
    NDSolveValue[{eq, ic, bc}, {ρ, T, u}, {t, 0, endtime}, {x, xL, xR}, 
     Method -> Union[mol[False], mol[points, difforder]], 
     MaxStepFraction -> {1/10^4, Infinity}, MaxSteps -> Infinity]; // AbsoluteTiming
(* {8.111344, Null} *)

(* Please find definition of sol in user21's answer. *)
  Plot[{sol[[1, i, -1]], mysol[[i]][t, x]} // Evaluate, {x, xL, xR}, 
    PlotRange -> All, PlotStyle -> {Automatic, Directive[{Thick, Dashed}]}], 
   {i, 3}], {t, 0, endtime}]

enter image description here

Just for comparison, in v9.0.1, ruebenko's solution takes about 35 seconds to finish computing with nxy = 33, and about 170 seconds with nxy = 43. (Yes, the speed drops dramatically when nxy increases. )


  1. I choose T instead of te when coding equation in this answer, because it's more straightforward and I believe T is unlikely to be introduced as a built-in symbol in the future.

  2. MaxStepFraction option can be taken away, and NDSolveValue will solve the system in less than 3 seconds then, but the solution will be slightly noisy in some region, for example:

    Plot[{sol[[1, 3, -1]], mysol[[3]][t, x]} /. t -> 0.187 // Evaluate, {x, xL, xR}, 
     PlotRange -> All, PlotStyle -> {Automatic, Directive[{Thick, Dashed}]}]

    Mathematica graphics

  3. mol[False] i.e. "DifferentiateBoundaryConditions" -> False can be taken away, but NDSolveValue will be slower without it.

  4. The difference between sol and mysol is relatively obvious in certain moment e.g. t == 0.187, but further check by adjusting points shows mysol seems to be the reliable one.

  5. Though I can't figure it out at the moment, I suspect there exists even more suitable way to discretize the system, given the precedent here.

  • 7
    $\begingroup$ It's always nice when people come back to questions they previously asked, armed with better understanding. :) $\endgroup$ Oct 4, 2018 at 16:09

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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