# NDSolve gives unexpected results when using the method of lines

I've been trying to solve the following PDE using the NDSolve function but it seems something is not working properly. The PDE is a the heat equation on polar coordinates and assuming angular simmetry:

$$u_t (t, r) = \alpha(r) \frac{1}{r} \partial_r (r u_r (t, r))$$

with boundary conditions $$u_t(t, r_{in}) = p$$ and $$u(t, r_{max}) = 0$$, and initial condition $$u(0,r)=0$$ for $$r\in(r_{in}, r_{max})$$. The function $$\alpha$$ represents the diffusion coefficient of the temperature in the media.

The idea is that $$r_{max}$$ is large enough so that $$u(t,r)$$ is alsmost flat and equal to zero close to $$r_{max}$$ through all the time window integrated.

I have tried different $$\alpha(r)$$ functions, which are

alpStep[r_] :=If[r < rout, aA, aB];
alpLin[r_] :=If[r < rout, aA + (aB - aA)(r-rin)/(rout - rin), aB];
alpA[r_] := aA;
alpB[r_] := aB;



where $$r_{out}\in(r_{ín},r_{max})$$ is a critical radius beyond which the diffusion remains constant.

The unexpected result is the following. I set aA > aB, and I integrate the different PDEs (with the different alp functions) by means of NDSolve. It turns out that the solution associated to the alpStep function takes values much more higher than the other functions, and I would expect it to be somewhere in between the solutions associated to the functions alpA and alpB. May be the problem is due to the discontinuous diffusion, but I can't see why.

The code is the following:

rin = 0.05;
rout = 0.15;
rmax = 5;
p = 0.01;
tend = 24*10;

aA = 0.01;
aB = 0.001;
alpStep[r_] := If[r < rout, aA, aB];
alpLin[r_] := If[r < rout, aA + (aB - aA) (r - rin)/(rout - rin), aB];
alpA[r_] := aA;
alpB[r_] := aB;

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

With[{u = u[t, r]}, eqn = alpStep[r] ((1/r) D[r D[u, r], r]) - D[u, t];
robinbc = NeumannValue[p*alpStep[rin], r == rin];
bc = DirichletCondition[u == 0, r == rmax];
ic = u == 0 /. {t -> 0};]
solStepAB =
NDSolveValue[{eqn == robinbc, bc, ic},
u, {t, 0, tend}, {r, rin, rmax}, opts];

With[{u = u[t, r]}, eqn = alpLin[r] ((1/r) D[r D[u, r], r]) - D[u, t];
robinbc = NeumannValue[p*alpLin[rin], r == rin];
bc = DirichletCondition[u == 0, r == rmax];
ic = u == 0 /. {t -> 0};]
solLinAB =
NDSolveValue[{eqn == robinbc, bc, ic},
u, {t, 0, tend}, {r, rin, rmax}, opts];

With[{u = u[t, r]}, eqn = alpA[r] ((1/r) D[r D[u, r], r]) - D[u, t];
robinbc = NeumannValue[p*alpA[rin], r == rin];
bc = DirichletCondition[u == 0, r == rmax];
ic = u == 0 /. {t -> 0};]
solConA =
NDSolveValue[{eqn == robinbc, bc, ic},
u, {t, 0, tend}, {r, rin, rmax}, opts];

With[{u = u[t, r]}, eqn = alpB[r] ((1/r) D[r D[u, r], r]) - D[u, t];
robinbc = NeumannValue[p*alpB[rin], r == rin];
bc = DirichletCondition[u == 0, r == rmax];
ic = u == 0 /. {t -> 0};]
solConB =
NDSolveValue[{eqn == robinbc, bc, ic},
u, {t, 0, tend}, {r, rin, rmax}, opts];

Grid[{{
Plot[{solStepAB[t, rin], solLinAB[t, rin], solConA[t, rin],
solConB[t, rin]}, {t, tend/1000, tend}],
Plot[{solStepAB[t, rout], solLinAB[t, rout], solConA[t, rout],
solConB[t, rout]}, {t, tend/1000, tend}],
Plot[{solStepAB[t, rout + (rout - rin)],
solLinAB[t, rout + (rout - rin)],
solConA[t, rout + (rout - rin)],
solConB[t, rout + (rout - rin)]}, {t, tend/1000, tend}]
}}]



The picture with the temperature profiles for the different solutions alpStep (blue), alpLin (Orange), alpA (Green) and alpB (Red) at $$r = r_{in}, r_{out}, 2 r_{out} - r_{in}$$ from left to right is the following, where it can be seen that the profile associated to alpStep is much higher that the others:

• I am not a 100% I understand your question but have a look at FEMDocumentation/tutorial/FiniteElementBestPractice#588198981 in the help system or online here Apr 25, 2019 at 5:52
• Maybe similar to this Apr 25, 2019 at 6:23
• Thanks for the reply, I'll take a look. Apr 25, 2019 at 13:44
• Were you able to figure this out? Apr 29, 2019 at 5:16

Yes! The problem was that some information was missing in the problem I stated above. Specifically, due to the discontinuity of the function $$\alpha$$ at $$r_{out}$$, the solution fails to be differentiable at that point, and multiple solutions exist. NDSolve gives one of these following some criterium related to fluxes (the heat flux at one side of the discontinuity has to be equal to the heat flux at the other side). Notice that the flux depends on the thermal conductivity $$k$$, which is not explicitly stated in the problem above. In general the thermal diffusivity $$\alpha$$ is equal to the conductivity $$k$$ over the heat capacity $$c$$ (i.e. $$\alpha = k / c$$). In order to use properly the NDSolve function in this case (in the sense that the continuity of the flux is preserved at $$r_{out}$$) one should avoid working with diffusion coefficient and use $$k$$ and $$c$$ instead, that is one should write
With[{u = u[t, r]}, eqn = k[r] ((1/r) D[r D[u, r], r]) - c[r] D[u, t];