# How to handle discontinuity in diffusion coefficient?

I am looking to solve the diffusion equation with a discontinuous jump in the diffusion coefficient. In 1D, the diffusion equation for $$u(t,x)$$ is: $$\partial_t u = \partial_x (D \partial_x u),$$ where $$D(x)$$ is the spatially varying diffusion coefficient. Let's use $$D(x) = 1$$ if $$x < 1$$ and $$D(x) = 3$$ if $$x > 1$$.

Question: Is there a better/smarter way to handle the discontinuity than approximating the jump with a continuous function? Is there a way to solve the equation in a piecewise manner, on $$x \in [0,1)$$ with $$D=1$$, on $$x \in (1,2]$$ with $$D=3$$, and somehow impose the condition that $$D^\text{left} \partial_x u^\text{left} = D^\text{right} \partial_x u^\text{right}$$ at $$x=1$$?

What I am trying to avoid is the following approximation:

We can approximate the discontinuity in $$D$$ by a sharp continuous function.

diffConst[x_] := (1 + 2 LogisticSigmoid[50 (x - 1)])

Plot[diffConst[x], {x, 0, 2}, PlotRange -> {0, 3}] Then we can solve the equation like so, with a sufficiently fine-grained spatial discretization:

fun = NDSolveValue[
{
D[z[x, t], t] == D[diffConst[x] D[z[x, t], x], x],
z[0, t] == 2,
z[2, t] == 1,
z[x, 0] == 2
},
z,
{x, 0, 2}, {t, 0, 20},
Method -> {"PDEDiscretization" -> {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid",
"MinPoints" -> 300}}}
]


Plot the solution:

Animate[
Plot[fun[x, t], {x, 0, 2}],
{t, 0, 5}
] • In principle, FEM would not require any continuity of the diffusion constant (just make sure that there is a mesh node at the jump point). But I am not sure whether one can access this feature from the high-level user interface... Apr 7, 2020 at 11:34
• Follow-up question: mathematica.stackexchange.com/questions/218983/… Apr 7, 2020 at 11:50
• @Henrik Do you have any tips on the follow-up question I just linked? Apr 7, 2020 at 11:50
• Somehow I missed this: reference.wolfram.com/language/PDEModels/tutorial/… It almost certainly has comparable examples, but I will need some time to go through it. Apr 7, 2020 at 13:23
• This is explained, for example, in the section Partial Differential Equations with Variable Coefficients Apr 8, 2020 at 7:55

Yes, it is possible to solve the equation in a piecewise manner :

diffConst[x_]:=If[x<1,1,3];

fun = NDSolveValue[
{
D[z[x, t], t] == D[diffConst[x] D[z[x, t], x], x],
z[0, t] == 2,
z[2, t] == 1,
z[x, 0] == 2
},
z,
{x, 0, 2}, {t, 0, 20},
Method -> {"PDEDiscretization" -> {"MethodOfLines",
"SpatialDiscretization" -> {"FiniteElement"(*,
"MinPoints" -> 300*)}}}
]
Plot[{fun[x,0.4]},{x,0,2}] Note that the function diffConst[x] D[z[x, t], x] is continuous :

derivative = D[fun[x, .4], {x}]
Plot[{diffConst[x] derivative}, {x, 0, 2}] This is very likely the solution you are expecting, because it corresponds to most physical situations. (For example, in thermics, It means conservation of the flux of heat thru the point x=1)

• Specifying FiniteElement is critical though, otherwise we get a continuous derivative, not a continuous flux Apr 7, 2020 at 12:10
• related Apr 7, 2020 at 12:26
• @Szabolcs, please see the section What Triggers the Use of the Finite Element Method in the documentation. Apr 8, 2020 at 7:57