8
$\begingroup$

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}]

enter image description here

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}
]

enter image description here

$\endgroup$
10
$\begingroup$

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}]  

enter image description here

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}]    

enter image description here

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)

| improve this answer | |
$\endgroup$

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.