NDSolve gives unexpected results when using the method of lines - Mathematica Stack Exchange most recent 30 from mathematica.stackexchange.com 2019-09-23T06:26:09Z https://mathematica.stackexchange.com/feeds/question/196921 https://creativecommons.org/licenses/by-sa/4.0/rdf https://mathematica.stackexchange.com/q/196921 0 NDSolve gives unexpected results when using the method of lines Carles BB https://mathematica.stackexchange.com/users/65195 2019-04-24T10:06:08Z 2019-05-29T11:01:41Z <p>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:</p> <p><span class="math-container">$u_t (t, r) = \alpha(r) \frac{1}{r} \partial_r (r u_r (t, r))$</span></p> <p>with boundary conditions <span class="math-container">$u_t(t, r_{in}) = p$</span> and <span class="math-container">$u(t, r_{max}) = 0$</span>, and initial condition <span class="math-container">$u(0,r)=0$</span> for <span class="math-container">$r\in(r_{in}, r_{max})$</span>. The function <span class="math-container">$\alpha$</span> represents the diffusion coefficient of the temperature in the media.</p> <p>The idea is that <span class="math-container">$r_{max}$</span> is large enough so that <span class="math-container">$u(t,r)$</span> is alsmost flat and equal to zero close to <span class="math-container">$r_{max}$</span> through all the time window integrated.</p> <p>I have tried different <span class="math-container">$\alpha(r)$</span> functions, which are</p> <pre><code>alpStep[r_] :=If[r &lt; rout, aA, aB]; alpLin[r_] :=If[r &lt; rout, aA + (aB - aA)(r-rin)/(rout - rin), aB]; alpA[r_] := aA; alpB[r_] := aB; </code></pre> <p>where <span class="math-container">$r_{out}\in(r_{ín},r_{max})$</span> is a critical radius beyond which the diffusion remains constant.</p> <p>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.</p> <p>The code is the following:</p> <pre><code>rin = 0.05; rout = 0.15; rmax = 5; p = 0.01; tend = 24*10; aA = 0.01; aB = 0.001; alpStep[r_] := If[r &lt; rout, aA, aB]; alpLin[r_] := If[r &lt; rout, aA + (aB - aA) (r - rin)/(rout - rin), aB]; alpA[r_] := aA; alpB[r_] := aB; opts = Method -&gt; {"MethodOfLines", "SpatialDiscretization" -&gt; {"FiniteElement", "MeshOptions" -&gt; {"MaxCellMeasure" -&gt; 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 -&gt; 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 -&gt; 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 -&gt; 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 -&gt; 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}] }}] </code></pre> <p>The picture with the temperature profiles for the different solutions alpStep (blue), alpLin (Orange), alpA (Green) and alpB (Red) at <span class="math-container">$r = r_{in}, r_{out}, 2 r_{out} - r_{in}$</span> from left to right is the following, where it can be seen that the profile associated to alpStep is much higher that the others:</p> <p><a href="https://i.stack.imgur.com/Uz67U.jpg" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/Uz67U.jpg" alt="enter image description here"></a></p> https://mathematica.stackexchange.com/questions/196921/-/197285#197285 1 Answer by Carles BB for NDSolve gives unexpected results when using the method of lines Carles BB https://mathematica.stackexchange.com/users/65195 2019-04-29T10:44:06Z 2019-04-29T10:44:06Z <p>Yes! The problem was that some information was missing in the problem I stated above. Specifically, due to the discontinuity of the function <span class="math-container">$\alpha$</span> at <span class="math-container">$r_{out}$</span>, 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 <span class="math-container">$k$</span>, which is not explicitly stated in the problem above. In general the thermal diffusivity <span class="math-container">$\alpha$</span> is equal to the conductivity <span class="math-container">$k$</span> over the heat capacity <span class="math-container">$c$</span> (i.e. <span class="math-container">$\alpha = k / c$</span>). In order to use properly the NDSolve function in this case (in the sense that the continuity of the flux is preserved at <span class="math-container">$r_{out}$</span>) one should avoid working with diffusion coefficient and use <span class="math-container">$k$</span> and <span class="math-container">$c$</span> instead, that is one should write</p> <pre><code>With[{u = u[t, r]}, eqn = k[r] ((1/r) D[r D[u, r], r]) - c[r] D[u, t]; robinbc = NeumannValue[p*k[rin], r == rin]; bc = DirichletCondition[u == 0, r == rmax]; ic = u == 0 /. {t -&gt; 0};] </code></pre>