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:
Is this behavior of
NDSolve
intended, or kind of a mistake?Is this behavior controlable? I mean, can we predict what's continuous in the solution, just from the form of the equation?
x[x,0]=Exp[x]
is not compatible with flux continuity. Same probleme here and here !! $\endgroup$ – andre314 Dec 18 '16 at 14:12