# Heat equation with variable conductivity

Let us consider the (very simple) problem

$$\qquad \partial_x[k(x)\,\partial_xu]=0,\quad x\in [0,1]$$

$$\qquad u(0)=0,\ u(1)=1,$$ with $$l\in(0,1)$$

$$\qquad k(x)=\left\{\begin{array}{l}k_1^{},x

which has the analytical solution given by

$$\qquad k(x)\,\partial_xu=a\rightarrow u(x)=\left\{\begin{array}{l}\frac{a}{k_1^{}}x,\,x

and for which $$a$$ must be sought as satistfying the boundary condition for $$x=1$$

$$\qquad a=\frac{1}{\frac{l}{k_1^{}}+\frac{1-l}{k_2^{}}}$$.

The following code

k1 = 1;
k2 = 2;
l = 3/4;
a = 1/(l/k1 + (1 - l)/k2)
k[x_] := Piecewise[{{k1, x < l}}, k2]
sol = DSolveValue[D[k[x] D[y[x], x], x] == 0 && y[0] == 0 && y[1] == 1, y[x], x]


yields $$x$$ as its result, which is not the solution for $$k_1\neq k_2$$ since the derivative of $$u$$ must change in passing from one conductivity to the other while the heat flux is constant.

Is this because a misuse of Piecewise? Or what else?

I am using version 10.3.1.

• What about an idea of making a continuous function with a narrow kink? You may use something as k[x_]:=k1+k2*(Tanh[(x-l)/d]+1)/2 or something else suitable. Here I assumed the case k2>k1. Then you can play with d decreasing it and see, if the solution converges. – Alexei Boulbitch Mar 18 '19 at 12:11
• Hi Alexei, thanks for your suggestion. Beyond the fact that it should work, I'd like to know if the behavior I found with my implementation is expected or not. – Fabio Mar 18 '19 at 15:10

You can help Mathematica a little bit to get the solution. Since $$(k(x) u_x)_x=0$$ then this means $$k(x) u_x=c$$ where $$c$$ is constant. So we solve this ODE with the left boundary condition only, then at the end solve for $$c$$ from the second boundary condition

ClearAll[k, y, x, c];
k1     = 1;
k2     = 2;
L      = 3/4;
k[x_]  := Piecewise[{{k1, x < L}, {k2,True}}];
bc     = y[0] == 0;
sol    = DSolveValue [{k[x] D[y[x], x] == c, bc}, y[x], {x, 0, 1}];
foundConstant = c /. First@Solve[(sol[[2]] == 1) /. x -> 1, c];
sol /. c -> foundConstant


Compare to the analytical solution you showed

a = 1/(L/k1 + (1 - L)/k2);
analytical = Piecewise[{{a/k1 x, x < L}, {(a/k1) L + (a/k2) (x - L), x > L}}]//Simplify


Here is a review of differents approachs.

NDSolve with Finite Elements

This is the result the OP expects. Temperature and flux are continuous.
Unfortunetally it is only a numerical solution.

k[x_] := Piecewise[{{1, x < 3/4}}, 2];
sol00 = NDSolveValue[
D[k[x] D[y[x], x], x] == 0 && y[0] == 0 && y[1] == 1,
y[x], {x, 0, 1},
Method -> {"PDEDiscretization" -> {"FiniteElement"}}];
Plot[Evaluate[With[{sol = sol00}, {sol, k[x] D[sol, x]}]], {x, 0, 1},
PlotLegends -> {"temperature", "flux"}, GridLines -> Automatic]


NDSolve with default Method option

The default option is : Method -> {"MethodOfLines", {"SpatialDiscretization" -> {"TensorProductGrid"}}}

This is a bad choice for the OP's problem :

k[x_] := Piecewise[{{1, x < 3/4}}, 2];
sol01 = NDSolveValue[
D[k[x] D[y[x], x], x] == 0 && y[0] == 0 && y[1] == 1,
y[x], {x, 0, 1},
Method -> {"MethodOfLines", {"SpatialDiscretization" -> \
{"TensorProductGrid"}}}];
Plot[Evaluate[With[{sol = sol01}, {sol, k[x] D[sol, x]}]], {x, 0, 1},
PlotLegends -> {"temperature", "flux"}, GridLines -> Automatic]


This behaviour is normal because this method is the Finite Difference method, which assumes that the temperature and its derivative are continuous (if the derivative is continuous then k[x] times the derivative is not continuous).

OP's attempt with NDSolve

k[x_] := Piecewise[{{1, x < 3/4}}, 2];
sol02 = DSolveValue[
D[k[x] D[y[x], x], x] == 0 && y[0] == 0 && y[1] == 1, y[x], x];
Plot[Evaluate[With[{sol = sol02}, {sol, k[x] D[sol, x]}]], {x, 0, 1},
PlotLegends -> {"temperature", "flux"}, GridLines -> Automatic]


Attempt with DSolve and UnitStep :

k[x_] := 1 + UnitStep[x - l];
sol03 = DSolveValue[
D[k[x] D[y[x], x], x] == 0 && y[0] == 0 && y[1] == 1, y[x], x];
Plot[Evaluate[With[{sol = sol03}, {sol, k[x] D[sol, x]}]], {x, 0, 1},
PlotLegends -> {"temperature", "flux"}, GridLines -> Automatic]


Attempt with DSolve and HeavisideTheta :

k[x_] := 1 + HeavisideTheta[x - l];
sol04 = DSolveValue[
D[k[x] D[y[x], x], x] == 0 && y[0] == 0 && y[1] == 1, y[x], x];
Plot[Evaluate[With[{sol = sol04}, {sol, k[x] D[sol, x]}]], {x, 0, 1},
PlotLegends -> {"temperature", "flux"}, GridLines -> Automatic]


Mathematically, the solution is indeed not unique without restrictive assumptions. Physically, it is obvious that the solution should be continuous nearly everywhere. "nearly everywhere" but not "absolutely everywhere", because one can easily imagine real physical problems that are modeled with some discontinuity at material interfaces. So mathematica really can't guess that you are expecting the continuity of flux.

The problem is that these different behaviours are explained nowhere in the documentation. This is probably because the explanations rely on some depth understandings. For example, it seems that the finite element method always respects the continuity of flux (more precisely the term "k[x] Grad[u[x]]" in "Div[ k[x] Grad[u[x]]]"), but I don't see why.

One more thing : One can see a stiff step as a limit case of a serie of smooth steps. One is then tempted to think that the solution of the underspecified initial problem is unique. This is false.

• (+1) Great summary. – user21 Mar 19 '19 at 5:33
• How about adding a solution using continuous approximation of UnitStep? – xzczd Mar 19 '19 at 6:31
• andre314 thanks for the detailed answer... the only behavior I recognize as "real" is the FEM solution. I think that the continuity is ensured by the fact that at x=l the variational formulation constrains uL JL=uR JR. since for x=l there is only one DOF must be uL=uR and as a consequence JL=JR. – Fabio Mar 19 '19 at 11:51