I attempted to use NDSolve for the 1-D isentropic unsteady flow equations with low subsonic inflow velocity and prescribed inflow total enthalpy; along with a non-reflective subsonic outflow boundary condition (BC) based on the method of characteristics. The latter BC is commonly used in Computational Fluid Dynamics (CFD) and in Computational Aero-Acoustics. NDSolve rejects this BC outright in Mathematica 9.0.1.


Auxiliary Eq's and parameters.

gamma = 1.4;
sigma = 0.25;
uInfty = 0.02;
M$ref = uInfty/Sqrt[gamma - 1];
dhInfty = (-(1/2)) uInfty^2;
h[x_, t_] := 1 + dh[x, t];
c[x_, t_] := ((gamma - 1)*h[x, t])^(1/2);
k = sigma*(1 - M$ref^2);
rule = {
 residu -> D[dh[x, t] + (1/2)*u[x, t]^2, x],
 residh -> u[x, t]*D[dh[x, t], x] + (gamma - 1)*D[u[x, t], x]};

PDE's for dependent variables u and dh vs. x and t.

AllEqs = {
 ueqn -> D[u[x, t], t] + residu == 0,
 dheqn -> D[dh[x, t], t] + residh == 0} /. rule;

BC and IC.

BCICRules = {
 bcInFlow -> {(*Inflow Boundary x=0 (subsonic)*)
  u[0, t] == uInfty,
  dh[0, t] + (1/2)*u[0, t]^2 == 0},
 bcOutFlow -> {(*Outflow Boundary x=1 (subsonic, non-reflective)*)
  Derivative[0, 1][dh][1, t] - c[1, t]*Derivative[0, 1][u][1, t] + 
   k*c[1, t]*(dh[1, t] - dhInfty) == 0},
 ic -> {(*IC*)
   u[x, 0] == (1 - x)*uInfty,
   dh[x, 0] == (1 - x)*dhInfty}


Solve the equations numerically to steady state

NDSolve[Flatten@({ueqn, dheqn, bcInFlow, bcOutFlow, ic} /. AllEqs /. BCICRules), 
{dh[x, t], u[x, t]}, {t, 0, 3.}, {x, 0, 1.}]

Following error occurs!

NDSolve::bdord: Boundary condition 0.157956 (0.0002 +dh[1,t]) Sqrt[1+dh[1,t]]+(dh^(0,1))[1,t]-0.632456 Sqrt[1+dh[1,t]] (u^(0,1))[1,t] should have derivatives of order lower than the differential order of the partial differential equation. >>


  • The reason for failure is bcOutFlow specification. Changing it to the following by vanishing the time derivatives results into a solution without any complain.

New bcOutFlow.

bcOutFlow -> {
  Derivative[0, 0][dh][1, t] - c[1, t]*Derivative[0, 0][u][1, t] + 
  k*c[1, t]*(dh[1, t] - dhInfty) == 0}

Here is the plot.

Plot3D[Evaluate[{#[x, t]} /. soln0], {t, 0, 3}, {x, 0, 1}, 
Mesh -> None, ColorFunction -> "DarkRainbow", 
PlotStyle -> Opacity[.7], PlotPoints -> 40, ImageSize -> 400, 
PlotLabel -> #] & /@ {dh, u} // Row

enter image description here

  • Above error states that BC is not allowed to contain time derivatives of the same order as the time derivatives in the PDE's (namely, first-order time derivatives). THIS CLAIM IS PATENTLY INCORRECT FOR THE COMPRESSIBLE FLOW EQUATIONS. The non-reflective boundary condition containing time derivatives is widely used, as anyone can verify from the CFD literature of the past two decades.
  • I had hoped that, within Mathematica's Method of Lines, the non-reflective outflow boundary condition should just give rise to another ODE in time at the outflow-boundary grid point $x=1$. That ODE should be coupled automatically to those at the interior grid points through the stencil of the built-in default spatial differencing scheme.
  • The pre-test has caused Mathematica to give the NDSolve::bdord: error and reject the problem outright. If the pre-test were bypassed, then there might be hope that NDSolve would translate the time-derivative outflow boundary condition into a time-dependent ODE for the dependent variables at the boundary grid point $x=1$. It could then go forward with the numerical integration of the ODE system.


  1. Does anyone out there know of a way to circumvent this difficulty with the NDSolve implementation in Mathematica 9.0.1?
  2. Specifically, does anyone know if it is possible to bypass NDSolve's pre-test for derivatives in the boundary conditions?
  3. Any suggestion how to accommodate reflective BC for a PDE like the above one in Mathematica?
  • 4
    $\begingroup$ I find two problems with your post. First, your code is unreadable. Please format it in the normal Mathematica input form. Second, it isn't clear what kind of help you are expecting. $\endgroup$
    – m_goldberg
    May 1, 2013 at 3:15
  • 2
    $\begingroup$ You should copy the InputForm version of your code to make it better readable. The warning message you get from Mathematica is probably not about what is theoretically possible, but about BC conditions that Mathematica can handle. $\endgroup$ May 1, 2013 at 5:22
  • $\begingroup$ From the Author: I improved code readability by converting Mathematica StandardForm commands to InputForm. I also tried to be more specific about what kind of help I hope someone can give. Would you please consider re-opening the Question? I hope my editing will answer the objections posed by commentors @m_goldberg and Sjoerd C. de Vries $\endgroup$ May 1, 2013 at 20:05

1 Answer 1


I managed to solve your equation set by

  1. adding a very small viscosity μ to the equation to circumvent the computational difficulty caused by the discontinuities of spatial derivative i.e. D[u[x, t], x] and D[dh[x, t], x] in the solution.

  2. discretizing the PDEs to a set of ODEs to circumvent the bdord warning, via pdetoode.

The following is the code:

gamma = 14/10; sigma = 25/100; uInfty = 2/100;
M$ref = uInfty/Sqrt[gamma - 1];
dhInfty = (-(1/2)) uInfty^2;
h[x_, t_] = 1 + dh[x, t];
c[x_, t_] = ((gamma - 1) h[x, t])^(1/2);
k = sigma (1 - M$ref^2);
residu = D[dh[x, t] + (1/2) u[x, t]^2, x]; 
residh = u[x, t] D[dh[x, t], x] + (gamma - 1) D[u[x, t], x];

With[{μ = 5 10^-4}, {ueqn, dheqn} = {D[u[x, t], t] + residu == μ D[u[x, t], x, x], 
    D[dh[x, t], t] + residh == 0}];
bcInFlow = {u[0, t] == uInfty, dh[0, t] + (1/2) u[0, t]^2 == 0}; 
bcOutFlow = {Derivative[0, 1][dh][1, t] - c[1, t] Derivative[0, 1][u][1, t] + 
    k c[1, t] (dh[1, t] - dhInfty) == 0};
ic = {u[x, 0] == (1 - x) uInfty, dh[x, 0] == (1 - x) dhInfty};

domain = {L, R} = {0, 1}; tend = 3;
points = 500; xdifforder = 8;
grid = Array[# &, points, domain];
var = {u /@ grid, dh /@ grid};
(* Definition of pdetoode isn't included in this code piece,
   please find it in the link above. *)
ptoo = pdetoode[{u, dh}[x, t], t, grid, xdifforder];
odeu = Delete[#, {{1}, {-1}}] &@ptoo@ueqn;
odedh = Delete[#, {{1}(*, {-1}*)}] &@ptoo@dheqn;
odeic = ptoo@ic;
{odebcin, odebcout} = ptoo@{bcInFlow, bcOutFlow};
sollst = NDSolveValue[{odeu, odedh, odeic, odebcin, odebcout}, var, {t, 0, tend}, 
   MaxSteps -> Infinity, SolveDelayed -> True];
{solu, soldh} = rebuild[#, grid, -1] & /@ sollst;
Manipulate[GraphicsRow[Plot[#[x, t], {x, L, R}, 
     PlotRange -> #2] & @@@ {{solu, {0, 0.03}}, {soldh, {0, 0.016}}}, 
  ImageSize -> 800], {t, 0, tend}]

enter image description here


  1. In principle SolveDelayed -> True is equivalent to Method -> {"EquationSimplification" -> "Solve"}, and the latter is recommended, but actually, when added as a counter measure to the warning ndsdtc, the latter will make NDSolve return unevaluated, at least in v9.0.1, I think it's a bug.

  2. The definition of odeu and odedh can also be odeu = Delete[#, {{1}(*,{-1}*)}] &@ptoo@ueqn; odedh = Delete[#, {{1}, {-1}}] &@ptoo@dheqn;, but this definition takes more time to solve.

  3. It seems that points should be large enough, and xdifforder should be a large enough even number.

  4. μ should be small enough, but not too small, or NDSolve may spit out the mconly warning and fails.


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.