Consider the following differential equation:

$$\begin{align*}&\rho C_p\left(\frac{\partial T}{\partial t}\right)=k\left[\frac{\partial^2 T}{\partial x^2}\right]+\dot{q}\\ &\text{at }x=0,\;\frac{\partial T}{\partial x}=0\\ &\text{at }x=1,\frac{\partial T}{\partial x}=C_1(T(t,1)-C_2)\\ &\text{at }t=0,T(0,x)=C_3 \end{align*}$$

c1 = -10;
c2 = 10;
c3 = 20;    
q[t, x] = 100000;
    heat = 
 NDSolve[{1591920 D[u[t, x], t] == .87 D[u[t, x], x, x] + q[t, x], 
       (D[u[t, x], x] /. x -> 0) == 0, (D[u[t, x], x] /. x -> 1) == 
    c1 (u[t, 1] - c2), u[0, x] == c3}, u, {t, 0, 600}, {x, 0, 1}]

The solver doesn't work and gives the warning:

NDSolve::ibcinc: Warning: boundary and initial conditions are inconsistent. >>

It is clear where the problem lies:

D[u[0, x], x] /. {x -> 1} /. heat
u[0.0, 1] - c2) /. heat



Causing (D[u[t, x], x] /. x -> 1) == c1 (u[t, 1] - c2) to evaluate as false

And rather than approach c2, the heat goes up: Note that editing c1 does affect the graph, while editing c2 does not.

Plot3D[Evaluate[u[t, x] /. heat], {t, 0, 600}, {x, 0, 1}, 
PlotRange -> All, AxesLabel -> {"t(s)", "x(m)"}, 
ColorFunction -> "TemperatureMap"]


Why does this boundary directive refuse to be set, with the whole domain of D[u[0, x], x] returning 0?


  • 1
    $\begingroup$ I've rewritten my answer below, aiming at making it a canonical one, have a look. Any suggestions are welcomed. $\endgroup$
    – xzczd
    Commented Oct 26, 2017 at 18:36

2 Answers 2


Changed the Derivative boundaries to NeumannValue and FEA needed to be specified.

c1 = -10;
c2 = 10;
c3 = 20;
q[t_, x_] := 100000;

heat = NDSolveValue[{
   1591920 D[u[t, x], t] - ((87/100) D[u[t, x], x, x] + q[t, x]) == 
    NeumannValue[0, x == 0] + NeumannValue[c1 (u[t, x] - c2), x == 1],
   u[0, x] == c3},
  u, {t, 0, 600}, {x, 0, 1}, 
  Method -> {"MethodOfLines",  "SpatialDiscretization" -> {"FiniteElement", 
             "MeshOptions" -> {"MaxCellMeasure" -> {"Length" -> 0.001}}}}]

Plot3D[Evaluate[heat[t, x]], {t, 0, 600}, {x, 0, 1}, PlotRange -> All,
  AxesLabel -> {"t(s)", "x(m)"}, ColorFunction -> "TemperatureMap"]

enter image description here


Adding appropriate units

rc = (8950) (385); (*density, Kg/m3 x heat capacity, J/Kg K*)
k = 385; (*thermal conductivity, W/m K*)

h = 10;(*heat transfer coefficient, W/m^2 K*)

ambT = 10; (*T, K*)
iniT = 20; (*T, K*)

q[t_, x_] := 0;  (*W/m^3*)

heat = NDSolveValue[{
   rc D[u[t, x], t] - (k D[u[t, x], x, x] + q[t, x]) == 
    NeumannValue[0, x == 0] + NeumannValue[-h (u[t, x] - ambT), x == 1],
   u[0, x] == iniT},
  u, {t, 0, 600}, {x, 0, 1}, 
  Method -> {"MethodOfLines", "TemporalVariable" -> t, 
    "SpatialDiscretization" -> {"FiniteElement", 
      "MeshOptions" -> {"MaxCellMeasure" -> {"Length" -> 0.001}}}}]

Plot3D[Evaluate[heat[t, x]], {t, 0, 600}, {x, 0, 1}, 
  PlotRange -> All, AxesLabel -> {"t(s)", "x(m)"}, 
 ColorFunction -> Function[{x, y, z}, ColorData["TemperatureMap"][z/40]], 
 ColorFunctionScaling -> False]

q = 0: Looks correct with slight cooling to the ambient temperature of 10.

enter image description here

q = 100000: A convection coefficient of 10 can't dissipate this source so the temperature rises.

enter image description here

q[t_, x_] := If[0 < t < 100, 100000, 0]

enter image description here

 Plot[Evaluate[heat[t, x]], {x, 0, 1}, PlotRange -> {{0, 1}, {0, 40}},
   PlotLabel -> t, AxesLabel -> {"x(m)", "T(K)"}, 
  ColorFunction -> Function[{x, y}, ColorData["TemperatureMap"][y/40]], 
  ColorFunctionScaling -> False], {t, 0, 600}]

Reference my answer here: Couple a PDE and ODE in NDSolve

  • $\begingroup$ Temp still seems to be increasing rather than decreasing when running Manipulate[ Plot[Evaluate[heat[t, x]], {x, 0, 1}, PlotRange -> All, AxesLabel -> {"t(s)", "x(m)"}, ColorFunction -> "TemperatureMap"], {t, 0, 600}] $\endgroup$
    – Feyre
    Commented Aug 5, 2016 at 21:12
  • $\begingroup$ However, the questioned problem does seem to be answered, it seems now a point around x=0.95 is playing havoc. $\endgroup$
    – Feyre
    Commented Aug 5, 2016 at 21:14
  • $\begingroup$ @Feyre I added some updates that may help. $\endgroup$
    – Young
    Commented Aug 5, 2016 at 21:53
  • $\begingroup$ Excellent work. The initial values weren't mine and I suspected them part of the problem. If I could I'd accept again. $\endgroup$
    – Feyre
    Commented Aug 6, 2016 at 10:03

Short Answer


Method -> {"MethodOfLines", 
         "DifferentiateBoundaryConditions" -> {True, "ScaleFactor" -> 1}}

inside NDSolve will resolve the problem. It's not necessary to set "ScaleFactor" to 1, it just needs to be a not-that-small positive number.

Long Answer

The answer for this problem is hidden in this obscure tutorial. I'll try my best to retell it in an easier-to-understand way.

Let's consider the following simpler initial-boundary value problem (IBVP) of heat equation that suffers from the same issue:

$$\frac{\partial u}{\partial t}=\frac{\partial^2 u}{\partial x^2}$$ $$u(0,x)=x(1-x)$$ $$u(t,0)=0,\ \frac{\partial u}{\partial x}\bigg|_{x=1}=0$$ $$t>0,\ 0\leq x\leq 1$$

Clearly $u(0,x)=x(1-x)$ and $\frac{\partial u}{\partial x}\bigg|_{x=1}=0$ is inconsistent. When you solve it with NDSolve / NDSolveValue, ibcinc warning will be spit out:

tend = 1; xl = 0; xr = 1;

With[{u = u[t, x]},
 eq = D[u, t] == D[u, x, x];
 ic = u == (x - xl) (xr - x) /. t -> 0;
 bc = {u == 0 /. x -> xl, D[u, x] == 0 /. x -> xr};]

sol = NDSolveValue[{eq, ic, bc}, u, {t, 0, tend}, {x, xl, xr}]


and further check shows the boundary condition (b.c.) $\frac{\partial u}{\partial x}\bigg|_{x=1}=0$ isn't satisfied at all:

Plot[D[sol[t, x], x] /. x -> xr // Evaluate, {t, 0, tend}, PlotRange -> All]

Mathematica graphics

Why does this happen?

The best way for explaination is re-implementing the method used by NDSolve in this case i.e. method of lines, all by ourselves.

As mentioned in the document, method of lines is a numeric method that discretizes the partial diffential equation (PDE) in all but one dimension and then integrating the semi-discrete problem as a system of ordinary differential equations (ODEs) or differential algebraic equations (DAEs). Here I discretize the PDE with 2nd order centered difference formula:

$$f'' (x_i)\simeq\frac{f (x_{i}-h)-2 f (x_i)+f (x_{i}+h)}{h^2}$$

formula = eq /. {D[u[t, x], t] -> u[x]'[t], 
   D[u[t, x], x, x] -> (u[x - dx][t] - 2 u[x][t] + u[x + dx][t])/dx^2}

points = 5;
dx = (xr - xl)/(points - 1);
ode = Table[formula, {x, xl + dx, xr - dx, dx}]
{u[1/4]'[t] == 16 (u[0][t]   - 2 u[1/4][t] + u[1/2][t]), 
 u[1/2]'[t] == 16 (u[1/4][t] - 2 u[1/2][t] + u[3/4][t]), 
 u[3/4]'[t] == 16 (u[1/2][t] - 2 u[3/4][t] + u[1][t])}

I've chosen a very coarse grid for better illustration. The initial condition (i.c.) should also be discretized:

odeic = Table[ic /. u[t_, x_] :> u[x][t] // Evaluate, {x, xl, xr, dx}]
{u[0][0] == 0, u[1/4][0] == 3/16, u[1/2][0] == 1/4, u[3/4][0] == 3/16, u[1][0] == 0}

We still need to deal with b.c.s. The Dirichlet b.c. doesn't need discretization:

bcnew1 = bc[[1]] /. u[t, x_] :> u[x][t]
u[0][t] == 0

The Neumann b.c. contains the derivative of $x$ so we need to discretize it with the one-sided difference formula:

$$f' (x_n)\simeq \frac{f (x_{n}-2h)-4 f (x_{n}-h)+3 f (x_n)}{2 h}$$

bcnew2 = bc[[2]] /. 
  D[u[t, x_], x_] :> (u[x - 2 dx][t] - 4 u[x - dx][t] + 3 u[x][t])/(2 dx)
2 (u[1/2][t] - 4 u[3/4][t] + 3 u[1][t]) == 0

"OK, 5 unknowns, 5 equations, we can now solve the system with any ODE solver! Just as NDSolve does!" Sadly you were wrong if you thought this statement is correct, because:

  1. Though {ode, odeic, bcnew1, bcnew2} is already a solvable system, it's not a set of ODEs, but DAEs. Notice here ODE refers to explicit ODE i.e. the coefficient of derivative term can't be $0$. Clearly, bcnew1 and bcnew2 don't explicitly contain the derivative of $t$.

  2. Though NDSolve is able to handle the DAE system {ode, odeic, bcnew1, bcnew2} directly, it doesn't solve the PDE in this way by default. Instead, it'll try to transform the DAE system to an explicit ODE system, probably because its ODE solver is generally stronger than the DAE solver (at least now).


To be precise, internally NDSolve doesn't "transform the DAE system to an explicit ODE system", it directly discretize the PDE to an explicit ODE system that can be viewed as a transformation of {ode, odeic, bcnew1, bcnew2}, i.e. {ode, odeic, bcnew1, bcnew2} isn't really generated, but this truth doesn't influence the discussion in this answer.

So, how does NDSolve transform the DAE system to ODE system? "That's simple! Just eliminate some of the variables in ode and odeic with bcnew1 and bcnew2! " Yeah this is a possible method, but not the one implemented in NDSolve. NDSolve has chosen a method that may be rather unusual at first glance. It mixes the original b.c. with its 1st order derivative respects to $t$. For our specific problem, the b.c. becomes:

odebc1 = D[#, t] + scalefactor1 # & /@ bcnew1
scalefactor1 u[0][t] + u[0]'[t] == 0
odebc2 = D[#, t] + scalefactor2 # & /@ bcnew2
2 scalefactor2 (u[1/2][t] - 4 u[3/4][t] + 3 u[1][t]) + 
   2 (u[1/2]'[t] - 4 u[3/4]'[t] + 3 u[1]'[t]) == 0

Where scalefactor1 and scalefactor2 are properly chosen coefficients.

It's not hard to notice this approach is systematic and easy to implement, and I guess that's the reason why NDSolve chooses it for transforming algebraic equations to ODEs. Nevertheless, this method has its disadvantage. The generated b.c. is equivalent to the original b.c., only if the original b.c. is continuous and i.c. is consistent with the b.c..

Let's use odebc1 as an example. In our case, bc[[1]] is continuous, and it's consistent with ic, so it can be easily rebuilt from odebc1:

DSolve[{odebc1, odeic[[1]]}, u[0][t], t]
{{u[0][t] -> 0}}

However, if the i.c. is something that isn't consistent with bc[[1]], for example u[0][0] == 1, the b.c. rebuilt from odebc1 will become:

DSolve[{odebc1, u[0][0] == 1}, u[0][t], t]
{{u[0][t] -> E^(-scalefactor1 t)}}

It's no longer equivalent to bc[[1]], but when scalefactor1 is a large positive number, this b.c. will converge to the original one.

Now here comes the key point. As stated in the document:

With the default "ScaleFactor" value of Automatic, a scaling factor of 1 is used for Dirichlet boundary conditions and a scaling factor of 0 is used otherwise.

i.e. scalefactor2 will be set to 0. Guess what b.c. will be rebuilt in this way?:

With[{scalefactor2 = 0}, 
 DSolve[{D[#, t] + scalefactor2 # & /@ bc[[2]], D[ic, x] /. x -> xr}, 
  D[u[t, x], x] /. x -> xr, t]]
{{Derivative[0, 1][u][t, 1] -> -1}}

It's a completely different b.c..

Back to the problem mentioned in the question, we can analyse its b.c. with the same method as above:

bcInQuestion = (D[u[t, x], x] /. x -> 1) == c1 (u[t, 1] - c2);
icInQuestion = u[0, x] == c3;
With[{sf = 0}, 
  DSolve[{D[#, t] + sf # & /@ bcInQuestion, D[icInQuestion, x] /. x -> 1}, 
   D[u[t, x], x] /. x -> 1, t]] /. u[0, _] :> c3
{{Derivative[0, 1][u][t, 1] -> -c1 c3 + c1 u[t, 1]}}

We see c1 is still in the rebuilt b.c., while c2 is completely killed by the derivation, that's why OP found "editing c1 does affect the graph, while editing c2 does not".

OK, then why does NDSolve choose such a strange setting for scaling factor? The document explains as follows:

There are two reasons that the scaling factor to multiply the original boundary condition is zero for boundary conditions with spatial derivatives. First, imposing the condition for the discretized equation is only a spatial approximation, so it does not always make sense to enforce it as exactly as possible for all time. Second, particularly with higher-order spatial derivatives, the large coefficients from one-sided finite differencing can be a potential source of instability when the condition is included. …

but personally I think this design is just too lazy: why not make NDSolve choose a non-zero scaling factor at least when ibcinc warning pops up and the order of spatial derivatives isn't too high (this is usually the case, the differential order of most of PDEs in practise is no higher than 2 )?

Anyway, now we know how to fix the issue. Just choose a positive scaling factor:

c1 = -10;
c3 = 20;
q[t, x] = 100000;
heat = ParametricNDSolveValue[{1591920 D[u[t, x], t] == .87 D[u[t, x], x, x] + 
      q[t, x], (D[u[t, x], x] /. x -> 0) == 0, (D[u[t, x], x] /. x -> 1) == 
     c1 (u[t, 1] - c2), u[0, x] == c3}, u, {t, 0, 600}, {x, 0, 1}, c2, 
   Method -> {"MethodOfLines", 
     "DifferentiateBoundaryConditions" -> {True, "ScaleFactor" -> 1}}];
Plot[heat[#][300, x] & /@ Range[10, 50, 10] // Evaluate, {x, 0.9, 1}, PlotRange -> All, 
 AxesLabel -> {"x(m)", "T"}]

Mathematica graphics

Now c2 influences the solution.

Young has also solved the problem with the new-in-v10 "FiniteElement" method, I guess it's probably because b.c.s are imposed in a completely different way when "FiniteElement" method is chosen, but I'd like not to talk too much about it given I'm still in v9 and haven't looked into "FiniteElement".

  • 3
    $\begingroup$ Thanks for the clear and in-depth explanation! $\endgroup$
    – lxy
    Commented Oct 23, 2018 at 6:30
  • $\begingroup$ This problem persists in v11.3 at least, although Plot[D[sol[t, x], x] /. x -> xr // Evaluate, {t, 0, tend}, PlotRange -> All] gives a different wrong curve $\endgroup$
    – lxy
    Commented Mar 18, 2023 at 9:37
  • $\begingroup$ I have read this post once more. Btw, on the LHS of that three FD formulae, a prime is missing in each one. $\endgroup$
    – lxy
    Commented Mar 18, 2023 at 10:02
  • $\begingroup$ please check your answer, where it states {u[1/4][t] == ..., {u[1/2][t] ==, ... {u[3/4][t] ==...} $\endgroup$
    – lxy
    Commented Mar 18, 2023 at 12:07
  • $\begingroup$ @jsxs Oops, I see. Edited. Thx for pointing out. For such mistakes it's OK to directly edit the post without asking me :) . $\endgroup$
    – xzczd
    Commented Mar 18, 2023 at 12:25

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.