4
$\begingroup$

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<l\\k_2^{},x\ge l\end{array}\right.\,,$

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<l\\\frac{a}{k_1^{}}l+\frac{a}{k_2^{}}(x-l),\,x\ge l\end{array}\right. \,,$

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.

$\endgroup$
  • 1
    $\begingroup$ 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. $\endgroup$ – Alexei Boulbitch Mar 18 at 12:11
  • $\begingroup$ 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. $\endgroup$ – Fabio Mar 18 at 15:10
4
$\begingroup$

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

Mathematica graphics

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 

Mathematica graphics

$\endgroup$
4
$\begingroup$

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]

enter image description here

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]  

enter image description here

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

just added for completeness :

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]  

enter image description here

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]  

enter image description here

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]

enter image description here

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.

$\endgroup$
  • 1
    $\begingroup$ (+1) Great summary. $\endgroup$ – user21 Mar 19 at 5:33
  • $\begingroup$ How about adding a solution using continuous approximation of UnitStep? $\endgroup$ – xzczd Mar 19 at 6:31
  • $\begingroup$ 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. $\endgroup$ – Fabio Mar 19 at 11:51

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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