Bug introduced in or after 10.3, persisting through 11.2.

I'm trying to solve following PDE (heat equation):

$$ \begin{cases} u_t = a \, u_{xx} \\ u(x,0) = 0\\ \lim_{x\to \infty}u(x,t) =0\\ \alpha\, (\theta_0-u(0,t))+\dot{q}_0=-\lambda u_x(0,t) \end{cases}$$

Where basically I have an initial temperature of $0$ everywhere, a constant heat flux at the beginning of the rod, and convection between air and the rod at its beginning ($\theta_0$ is the air temperature which is assumed to be constant).

I found following analytical solution for the problem:

$$u(x,t) = \frac{\dot{q}_0+\alpha \, \theta_0}{\alpha}\left[ \mathrm{erfc}\, \left(\frac{x}{2\sqrt{a \, t}} \right) -\mathrm{exp}\,\left(\frac{\alpha}{\lambda}x+a \frac{\alpha^2}{\lambda^2}t \right)\mathrm{erfc}\,\left(\frac{x}{2\sqrt{a\, t}} + \frac{\alpha}{\lambda} \sqrt{a\, t} \right) \right] $$

which is physically meaningful. With mathematica, however, I get some meaningless results, probably due to my not that good mathematica skills.

This is what I tried to do:

λ = 0.8; c = 880; ρ = 1950; a = λ/(c ρ)
α = 15; θair = 0;

heqn = D[u[x, t], t] == a D[u[x, t], {x, 2}];
ic1 = u[x, 0] == 0;
bc1 = α (θair - u[0, t]) + 650 == -λ Derivative[1, 0][u][0, t];

sol = DSolveValue[{heqn, ic1, bc1}, u[x, t], {x, t}][[1, 1, 1]]

Which leads me to a complex solution (plotted below)

DensityPlot[sol, {t, 0, 2*3600}, {x, 0, 0.1}, 
 PlotLegends -> Automatic, FrameLabel -> {t[s], x[m]}]

Density plot of the resulting equation[scale bar[2]

Plot[Evaluate[Table[sol, {t, 3600, 7200, 3600}]], {x, 0, .1}, 
 PlotRange -> All, Filling -> Axis, Axes -> True, AxesLabel -> {x[m], θ[°C]}]

Plot at t=1h and t=2h

I'm aware of the fact that I didn't consider the condition at infinity. To do this, I tried to follow this answer without success. Also, mathematica finds the solution without this condition as soon as $\alpha = 0$.

This is the plot of the analytical solution I get:

Plot[{u[x, 600], u[x, 3600], u[x, 7200]}, {x, 0, .2}, Filling -> Axis, 
  Axes -> True, AxesLabel -> {x[m], θ[°C]}]

Plot of analytical function

  • $\begingroup$ Do you want to obtain the same analytic solution with Mathematica, or just want to solve the problem numerically or analytically? $\endgroup$
    – xzczd
    Commented Jan 17, 2018 at 11:08
  • $\begingroup$ BTW, something seems to be wrong with the $\LaTeX$ formula, I failed to reproduce the graph with it: sola[t_, x_] = ( q + α θair)/α (Erfc[x/(2 Sqrt[a t])] - Exp[α/λ x + a α^2/λ^2 t] Erfc[ x/(2 Sqrt[α t]) + α/λ Sqrt[a t]]) $\endgroup$
    – xzczd
    Commented Jan 17, 2018 at 11:26
  • $\begingroup$ @xzczd at the moment I cannot check if I did some mistake by copying the analytical solution. I'll check it as soon as possible (probably in 1-2hours). Regarding your second question: I'm preparing a benchmark problem for a FEM simulation. A numerical solution is therefore in principle also ok. The only thing that seemed a bit strange to me is that Mathematica cannot solve a PDE which can be solved analytically (although I used Laplace transformations and this isn't probably the easiest thing for a software) $\endgroup$
    – David
    Commented Jan 17, 2018 at 11:44
  • 1
    $\begingroup$ @xzczd The latex formula is correct, you had a small mistake in your definition of the term Erfc[ x/(2 Sqrt[α t]) + α/λ Sqrt[a t]]. It should be x/(2 Sqrt[a t]) instead of x/(2 Sqrt[α t]). $\endgroup$
    – David
    Commented Jan 17, 2018 at 13:35
  • $\begingroup$ Oh, I see… I suggest adding the Mathematica code for the formula to the question so people can easily test it. BTW, thx for accepting, but you actually don't need to accept that fast, you can wait for 24 hours or so to see if someone will come up with a better answer :) . $\endgroup$
    – xzczd
    Commented Jan 17, 2018 at 13:39

2 Answers 2


Semi-analytic Solution

The approach in the linked post can be used for solving your problem analytically. We just need an extra step i.e. making Laplace transform:

heqn = D[u[x, t], t] == a D[u[x, t], {x, 2}];
ic1 = u[x, 0] == 0;
ic2 = α (θair - u[0, t]) + q == -λ Derivative[1, 0][u][0, t];

teqn = LaplaceTransform[{heqn, ic2}, t, s] /. Rule @@ ic1 /. 
  HoldPattern@LaplaceTransform[a_, __] :> a

tsol = Collect[u[x, t] /. First@DSolve[teqn, u[x, t], x], Exp[_], Simplify]
E^((Sqrt[s] x)/Sqrt[a]) C[1] + (E^(-((Sqrt[s] x)/Sqrt[
   a])) (s^(3/2) λ C[1] + Sqrt[a] (q + α (θair - s C[1]))))/(
 Sqrt[a] s α + s^(3/2) λ)

tsol is the transformed solution, with one constant to be determined. Since the solution should vanish at $\infty$, C[1] can only be 0:

tsolfunc[s_, x_] = tsol /. C[1] -> 0

The last step is to transform back, but sadly InverseLaplaceTransform can't handle it. Nevertheless, we've indeed obtained an analytic solution involving integration: $$ u=\frac{1}{2 \pi i } \int_{\gamma-i \infty }^{\gamma+i \infty } \frac{\sqrt{a} e^{-\frac{\sqrt{s} x}{\sqrt{a}}} (\alpha \theta_0+q)}{\sqrt{a} \alpha s+\lambda s^{3/2}} e^{st} \, ds $$

If we want to get the numeric result, we need to make use of numeric inverse Laplace transform package e.g. this:

λ = 8/10; c = 880; ρ = 1950; a = λ/(c ρ);
α = 15; θair = 0; q = 650;
(* Definition of FT isn't included in this post,
   please find it in the link above. *)
Clear@sol; sol[t_, x_] = FT[tsolfunc[#, x] &, t];

Plot[sol[#, x] & /@ {600, 3600, 7200} // Evaluate, {x, 0, .2}, Filling -> Axis, 
 AxesLabel -> {x[m], θ[°C]}]

Mathematica graphics

Analytic Solution

Based on this answer, the following is one way to find the symbolic Laplace inversion:

coe = (q + α θair);
mid = Apart[LaplaceTransform[tsolfunc[s, x]/coe, x, S] /. s -> sqrts^2] /. 
  sqrts -> Sqrt[s]

mid2 = InverseLaplaceTransform[mid, s, t]; // AbsoluteTiming
(* {33.611455, Null} *)
final = FullSimplify[InverseLaplaceTransform[mid2[[1]] // Expand, S, x], {t > 0, a > 0}]

coe final is the analytic solution:

coe final

$$u=\frac{(\alpha \theta_0+q)}{\alpha } \left(\text{erfc}\left(\frac{x}{2 \sqrt{a t}}\right)-e^{\frac{\alpha (a \alpha t+\lambda x)}{\lambda ^2}} \text{erfc}\left(\frac{2 a \alpha t+\lambda x}{2 \lambda \sqrt{a t}}\right)\right)$$


sola[t_, x_] = (q + α θair)/α (Erfc[x/(2 Sqrt[a t])] - 
    Exp[α/λ x + a α^2/λ^2 t] Erfc[x/(2 Sqrt[a t]) + α/λ Sqrt[a t]])

sola[t, x] == final coe // FullSimplify
(* True *)

Numeric Solution

NDSolve can also be used for solving the problem, of course. We just need to adjust the option a little:

inf = 0.2;

mol[tf : False | True, sf_: Automatic] := {"MethodOfLines",
  "DifferentiateBoundaryConditions" -> {tf, "ScaleFactor" -> sf}}

nsol = NDSolveValue[{heqn, ic1, ic2, u[inf, t] == 0}, u, {x, 0, inf}, {t, 0, 7200}, 
   Method -> mol[True, 100]];

np = Plot[nsol[x, #] & /@ {600, 3600, 7200} // Evaluate, {x, 0, inf}, Filling -> Axis, 
  AxesLabel -> {x[m], θ[°C]}}]

The result looks the same so I'd like to omit it here. To learn more about why the option needs to be adjusted, check this post.

  • $\begingroup$ Thanks for your answer! I accepted it after getting the correct semi-analytical solution. However, now I'm trying the numerical one and I still get a NDSolveValue::ibcinc error (bc and ic are inconsistent), which, according to your other post should be avoided by using the method of lines. $\endgroup$
    – David
    Commented Jan 17, 2018 at 16:34
  • 1
    $\begingroup$ @david Well, the linked post is not about avoiding ibcinc, it's about getting an accurate enough result even if it pops up. Notice ibcinc is merely a warning, it doesn't necessarily mean the result will be unreliable. (You can compare np to the analytic solution. ) $\endgroup$
    – xzczd
    Commented Jan 18, 2018 at 2:52
  • $\begingroup$ Sure, you're right. thanks $\endgroup$
    – David
    Commented Jan 18, 2018 at 7:39

Analytical solution seems to be wrong!

Just try

heqn = Derivative[0, 1][u][x, t] == a Derivative[2, 0][u][x, t];
ic1 = u[x, 0] == 0;
ic2 = \[Alpha] (\[Theta]air - u[0, t]) +qp0 == -\[Lambda] Derivative[1, 0][u][0, t];

sol = DSolve [{heqn, ic1, ic2}, u[x, t], {x, t}][[1]] 
(*{u[x, t] -> -((\[Alpha] (qp0 + \[Alpha] \[Theta]air) ((2 \[Alpha] +x \[Lambda]) Erfc[x/(2 Sqrt[a t])] - 
2 E^((\[Lambda] (x \[Alpha] +a t \[Lambda]))/\[Alpha]^2) \[Alpha] Erfc[(
    x \[Alpha] + 2 a t \[Lambda])/(2 Sqrt[a t] \[Alpha])]))/(2\[Lambda]^3))} *)

sol looks like the expected analytical solution!

u[infinity,t]->0 seems to be fullfilled:

Simplify[Normal[ Series[u[x, t] /. sol, {x, Infinity, 0}] ], {x > 0,t > 0, \[Alpha] > 0, \[Lambda] > 0, a > 0} ]
(*-((E^(-(x^2/(4 a t))) Sqrt[a t] \[Alpha] (qp0 + \[Alpha] \[Theta]air))/(Sqrt[\[Pi]] \[Lambda]^2))*)

because E^(-(x^2/(4 a t))) decreases with x!

The boundary condition u[x,0]==0 can be verified(a>0,x>0):

Simplify[Limit[u[x, t] /. sol[[1]], t -> 0], {x > 0, a > 0}]
(* 0 *)

but the solution doesn't fullfill the pde

(E^(-(x^2/(4 a t))) Sqrt[a t] \[Alpha] (qp0 + \[Alpha] \[Theta]air))/\[Lambda] == 0

Also ic2 isn't true!

  • $\begingroup$ Sadly it's not… you'll see the 2nd graph in the question if you plot it. $\endgroup$
    – xzczd
    Commented Jan 17, 2018 at 12:15
  • $\begingroup$ @xzczd hopefully it does. I edited my answer... $\endgroup$ Commented Jan 17, 2018 at 21:45
  • $\begingroup$ Well, but it doesn't satisfy the b.c. at $x=0$, try test = {u[x, t] -> -((\[Alpha] (qp0 + \[Alpha] \[Theta]air) ((2 \[Alpha] + x \[Lambda]) Erfc[x/( 2 Sqrt[a t])] - 2 E^((\[Lambda] (x \[Alpha] + a t \[Lambda]))/\[Alpha]^2) \[Alpha] Erfc[( x \[Alpha] + 2 a t \[Lambda])/(2 Sqrt[a t] \[Alpha])]))/( 2 \[Lambda]^3))} /. (u[x, t] -> rest_) :> (u -> Function @@ {{x, t}, rest}); Simplify[ic2 /. test, t > 0]. You can also plot the solution and compare it to the graph in the question. $\endgroup$
    – xzczd
    Commented Jan 18, 2018 at 3:06
  • $\begingroup$ @xzczd u[x,0}==0 !I edited my answer... $\endgroup$ Commented Jan 18, 2018 at 8:17
  • $\begingroup$ This is the initial condition… I mean the boundary condition i.e. ic2 in the code. $\endgroup$
    – xzczd
    Commented Jan 18, 2018 at 8:22

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.