Position of discontinuous coefficient influences the solution of PDE

This issue is raised in the discussion under this post about heat flux continuity and I think it's better to start a new question to state it in a clearer way. Just consider the following example:

Lmid = 1; L = 2; tend = 1;
m[x_] = If[x < Lmid, 1, 2];
eq1 = m[x] D[u[x, t], t] == D[u[x, t], x, x];
eq2 = D[u[x, t], t] == D[u[x, t], x, x]/m[x];

Clearly, eq1 and eq2 is mathematically the same, the only difference between them is the position of the discontinuous coefficient m[x]. Nevertheless, the solution of NDSolve will be influenced by this trivial difference, if "FiniteElement" is chosen as the method for "SpatialDiscretization":

opts = Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"FiniteElement",
"MeshOptions" -> {"MaxCellMeasure" -> 0.01}}};

ndsolve[eq_] := NDSolveValue[{eq, u[x, 0] == Exp[x]}, u, {x, 0, L}, {t, 0, tend}, opts];

{sol1, sol2} = ndsolve /@ {eq1, eq2};
Plot[{sol1[x, tend], sol2[x, tend]}, {x, 0, L}] Apparently sol2 is a weak solution that's just 0th order continuous in x direction.

Further check shows that, sol1 is 1st order continuous in x direction, while D[sol2[x, tend]/m[x], x] is continous:

Plot[D[{sol1[x, tend], sol2[x, tend]/m[x]}, x] // Evaluate, {x, 0, L}] To make this post a question, I'd like to ask:

1. Is this behavior of NDSolve intended, or kind of a mistake?

2. Is this behavior controlable? I mean, can we predict what's continuous in the solution, just from the form of the equation?

• @andre Well, maybe they don't like pictures, I tried to improve the format. I didn't make further changes because I'm not sure if I can make it more attractive. (As we can see this question of mine doesn't generate much enthusiasm, too…) – xzczd Nov 26 '16 at 6:10
• I am just realising that the initial condition x[x,0]=Exp[x] is not compatible with flux continuity. Same probleme here and here !! – andre314 Dec 18 '16 at 14:12
• @andre I think it's not a big deal. We can simply consider the inconsistency as an approximation, or a weak solution problem (here is another example). Also, this inconsistency isn't hard to remove. – xzczd Dec 18 '16 at 15:25
• @xzczd, thanks for the bounty! – user21 Dec 19 '16 at 13:15
• @user21 You deserve it :) . – xzczd Dec 19 '16 at 13:27

Here is an explanation of what happens. Let's setup the problem once more.

Lmid = 1; L = 2; tend = 1;
m[x_] = If[x < Lmid, 1, 2];
(*m[x_]=2;*)
eq1 = m[x] D[u[x, t], t] == D[u[x, t], x, x];
eq2 = D[u[x, t], t] == D[u[x, t], x, x]/m[x];
opts = Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"FiniteElement",
"MeshOptions" -> {"MaxCellMeasure" -> 0.01}}};
ndsolve[eq_] :=
NDSolveValue[{eq, u[x, 0] == Exp[x]}, u, {x, 0, L}, {t, 0, tend},
opts];

Equation 1 and 2 are mathematically the same, however, when we evaluate them we get different results as shown here:

sol1 = ndsolve[eq1];
Plot[sol1[x, tend], {x, 0, L}] sol2 = ndsolve[eq2];
Plot[sol2[x, tend], {x, 0, L}] What happens? Let's look at how the PDE gets parsed.

ClearAll[getEquations]
getEquations[eq_] := Block[{temp},
temp = NDSolveProcessEquations[{eq, u[x, 0] == Exp[x]},
u, {x, 0, L}, {t, 0, tend}, opts][];
temp = temp["FiniteElementData"];
temp = temp["PDECoefficientData"];
(# -> temp[#]) & /@ {"DampingCoefficients", "DiffusionCoefficients",
"ConvectionCoefficients"}
]

getEquations[eq1]
{"DampingCoefficients" -> {{If[x < 1, 1, 2]}},
"DiffusionCoefficients" -> {{{{-1}}}},
"ConvectionCoefficients" -> {{{{0}}}}}

This looks good.

getEquations[eq2]
{"DampingCoefficients" -> {{1}},
"DiffusionCoefficients" -> {{{{-(1/If[x < 1, 1, 2])}}}},
"ConvectionCoefficients" -> {{{{-(If[x < 1, 0, 0]/
If[x < 1, 1, 2]^2)}}}}}

For the second eqn. we get a convection coefficient term. Why is that? The key is to understand that the FEM can only solve this type equation:

$d\frac{\partial }{\partial t}u+\nabla \cdot (-c \nabla u-\alpha u+\gamma ) +\beta \cdot \nabla u+ a u -f=0$

Note, that there is no coefficient in front of the $\nabla \cdot (-c \nabla u-\alpha u+\gamma)$ term. To get things like $h(x) \nabla \cdot (-c \nabla u-\alpha u+\gamma)$ to work, $c$ is set to $h$ and $\beta$ is adjusted to get rid of the derivative caused by $\nabla \cdot (-c \nabla u)$

Here is an example:

c = h[x];
β = -Div[{{h[x]}}, {x}];
(* h[x]*Derivative[u][x] *)

In the case at hand that leads to:

Div[{{1/m[x]}}.Grad[u[x], {x}], {x}] -
Div[{{1/m[x]}}, {x}] // Simplify

(* {Piecewise[{{Derivative[u][x]/2, x >= 1}}, Derivative[u][x]]} *)

But that is the same as specifying:

eq3 = D[u[x, t], t] ==
Inactive[
Div][{{1/If[x < 1, 1, 2]}}.Inactive[Grad][u[x, t], {x}], {x}];

sol3 = ndsolve[eq3];
(* Plot[sol2[x, tend] - sol3[x, tend], {x, 0, L}] *)

I have checked that flexPDE (another FEM tool) gives exactly the same solutions in all three cases. So this issue is not uncommon. In principal a message could be generated but how would one detect when to trigger that message? If you have suggestions about this, let me know in the comments. I think it were also good to add this example to the documentation - if there are no objections. I hope this clarifies the unexpected behavior a bit.

• So the answer for my first question is "this behavior is intended"? "In principal a message could be generated but how would one detect when to trigger that message?" Maybe by detecting if the coefficient is continuous or not, or by detecting if the coefficient involves Piecewise etc.? – xzczd Dec 8 '16 at 11:28
• It's not exactly intended, but as I tried to explained it is a consequence of how the FEM works right now. Concerning the coefficient detection, what should happen in the case of f[x_?Numeric] or an interpolating function that has a discontinuity. Your "etc" is not as trivial as it may seem. I'll think a bit more about this and may I can come up with something. In the meanwhile I added this example to the documentation, pending approval, which I think it will get. – user21 Dec 12 '16 at 9:13

This is not a answer, only a comment. It is related to the continuity problem (see end of this comment).

The equations given by xzczd are the heat equation along a rod that has a thermal (volumic) capacity that double at point x=1. There are no boundaries conditions, so NDSolve[..., "FiniteElement"...] will take the Neuman boundaries conditions =0 (this is equivalent to thermal flux = 0, ie adiabatic boundaries). In this case the total heat quantity in the rod should stay constant over time. This quantity is very easy to calculate :

at t=0 :

NIntegrate[sol1[x, 0], {x, 0, 1}] +
2 NIntegrate[sol1[x, 0], {x, 1, 2}]

11.0598

at t=tend:

sol1 :

NIntegrate[sol1[x, 1], {x, 0, 1}] +
2 NIntegrate[sol1[x, 1], {x, 1, 2}]

11.0598

OK

sol2 :

NIntegrate[sol2[x, 1], {x, 0, 1}] +
2 NIntegrate[sol2[x, 1], {x, 1, 2}]

8.64626

KO

This problem is related to the continuity problem because if the continuity of flux=(conductivity*D[u,x]) fails at x=1 (conductivity = 1 here), then the global heat quantity is not conserved.

• "The equations given by xzczd are the heat equation along a rod that has a thermal (volumic) capacity that double at point $x=1$." Well, not necessarily. Actually m[x] represents $C/k$ where $C$ is thermal capacity and $k$ is thermal conductivity. So the equation in my question can also represent "a rod with a thermal conductivity that reduces by half at $x=1$", in this case the total heat quantity is NIntegrate[NIntegrate[temperature[time, x], {x, 0, 2}] and sol2 is (almost) the solution that obeys energy conservation. – xzczd Nov 19 '16 at 9:24
• Disagree : I'm agree that m[x] D[u[x, t], t] == D[u[x, t], x, x] and D[u[x, t], t] == D[u[x, t], x, x]/m[x] should be equivalent, but I still think that Div k[x] Grad u[x,t] is not equal to k[x] Div Grad u[x,t] even if k[x] is piecewise constant. That's the reason why I have assumed that conductivity=1 along the whole rod (may be you have not seen my last update above, where I have written conductivity=1 at the end of the text) – andre314 Nov 19 '16 at 12:00
• If you mean, the discontinuity in $k$ should be represented by HeavisideTheta rather than UnitStep so Div[k[x] Grad[u[t, x], {x, y}], {x, y}] actually involves DiracDelta, then, though I don't know how to prove it analytically, numeric experiments indeed show that DiracDelta has no contribution in the solution. Here's the code I used: appro = With[{coe = 1000}, ArcTan[coe #1]/π + 1/2 &]; unitStepExpand = SimplifyPWToUnitStep@PiecewiseExpand@# &;k[x_] = unitStepExpand[If[x < 1, 2, 1]] /. UnitStep -> appro;eq = D[u[t, x], t] == Div[k[x] Grad[u[t, x], {x, y}], {x, y}];…… – xzczd Nov 19 '16 at 12:50
• ……eqrefer = D[u[t, x], t] == k[x] Div[Grad[u[t, x], {x, y}], {x, y}];points = 5000;{sol, solrefer} = NDSolveValue[{#, u[0, x] == (1 - x) x, u[t, 0] == 0, u[t, 1] == 0}, u, {t, 0, 2}, {x, 0, 1}, Method -> {"MethodOfLines", "SpatialDiscretization" -> {"TensorProductGrid", "MaxPoints" -> points, "MinPoints" -> points, "DifferenceOrder" -> 4}}] & /@ {eq, eqrefer}; Plot[{sol[1/10, x], solrefer[1/10, x]}, {x, 0, 1}] You can adjust coe and points and observe, it's reasonable to think that the DiracDelta doesn't play a role in the solution. – xzczd Nov 19 '16 at 12:50
• @xzczd In the code of your 2 latest comments, the discontinuity is at x=1 and the right boundary is also at x=1. Is it a mistake or intentional ? – andre314 Nov 19 '16 at 17:34

The problem mentioned by the OP arises even with a ODE :

Compare :

m[x_]=If[x<0.5,1,2];

f=NDSolveValue[
{m[x] y''[x]==0, y==0, y==1},
y,
{x,0,1},
Method->{"PDEDiscretization"->{"FiniteElement"}}
];

Plot[f[x],{x,0,1}] with :

m[x_]=If[x<0.5,1,2];

f=NDSolveValue[
{y''[x]==0/m[x], y==0, y==1},
y,
{x,0,1},
Method->{"PDEDiscretization"->{"FiniteElement"}}
];

Plot[f[x],{x,0,1}] This time, we have the analytical solutions, which are trivial : y = a x + b, eventually in several segments.

• I think this case is a bit different. 0/m[x] will evaluate to 0 before being handled by NDSolve. – xzczd Jul 7 '17 at 2:55