# One term in non-linear equations prevents NDSolve to work

I am currently working on the stratification of the core of the planet Mercury, meaning the formation of a conductive layer at the top of the convective core, with a moving interface between both layers. After some variable changes to simplify my equations, here is the system of equations I want to solve:

$$\frac{\partial T}{\partial \tau}(\tau,x) = \frac{1}{(s(\tau)-1)^2}\frac{\partial^2 T}{\partial x^2}(\tau,x) +\left(\frac{2}{1+x(s(\tau)-1)}\frac{1}{s(\tau)-1}+\frac{x}{s(\tau)-1}\frac{\mathrm{d}s}{\mathrm{d}\tau}(\tau)\right) \frac{\partial T}{\partial x}(\tau,x)$$

$$\left(s(\tau) \left(\frac{T_{c0}}{T_{s0}-T_{c0}} + T_s(\tau)\right) s'(\tau) + \frac{1}{2 y}T_s'(\tau)\right)\left(2s(\tau) - \mathrm{e}^{y s^2(t)} \sqrt{\frac{\pi}{y}} \mathrm{erf}(\sqrt{y}s(\tau))\right) = 4 y s^3(\tau) \left(\frac{T_{c0}}{T_{s0}-T_{c0}}+T_s(\tau)\right)$$

$$\frac{1}{1-s(\tau)}\frac{\partial T}{\partial x}(\tau,0) = \frac{r_c}{k(T_{s0}-T_{c0})}q_c\left(\frac{\rho C_p r_c^2}{k}\tau+t_0\right)$$

$$T_s(\tau) = T(\tau,1)$$ $$\frac{1}{s(\tau)-1}\frac{\partial T}{\partial x}(\tau,1) = -2 y s(\tau)\left(\frac{T_{c0}}{T_{s0}-T_{c0}}+T_s(\tau)\right)$$

with $$T$$ the temperature profile in the conductive layer, $$s$$ the interface position, $$T_s$$ the interface temperature, $$T_{c0}$$, $$T_{s0}$$, $$r_c$$, $$k$$, $$\rho$$, $$C_p$$, $$t_0$$ and $$y = \frac{g_c \alpha r_c}{2 C_p}$$ constants.

I have discretised the spatial part of these equations in order to get a system of ODE's using the functions ptoo and ptoode (see here). Then I have used the function Solve in order to rewrite equations in the form '$$\frac{\mathrm{d}...}{\mathrm{d}\tau}=...$$' for all variables. And finally solve the equations using NDSolveValue. But I have got an error from NDSolve (see below).

If I delete the term $$\frac{\mathrm{d}s}{\mathrm{d}\tau}$$ in the right hand side of the first equation, everything goes fine and NDSolve solves my equations without complaining.

Is there something I can do in order to make NDSolve solving the system with the problematic term? I have tried to rearrange the equations in order to give simplified equations to NDSolve, changing the method (StiffnessSwitching, FixedStep, StartingStepSize or increasing the maximum number of steps) and I always have errors like 'max number of points reached' or 'stiff system'.

Here is my code:

(*constants*)
rc = 2050 10^3;
cp = 850;
rho = 7200;
alpha = 5 10^-5;
gc = 4;
k = 40;
y = (gc alpha rc)/(2 cp);

(*parameters*)
s0 = 2049 10^3;
Tc0 = 2100;
Ts0 = Exp[(-alpha gc)/(2 cp rc) (s0^2 - rc^2)] Tc0;
t0 = 0.099 10^9 365.25 24 3600;

qcmb[t_] =
With[{a = 0.004891658583550395, b = 0.34057028569554804,
c = 1.0021984665846737*^-15}, a + b E^(- c t)]

(*equations*)
s[τ] (Tc0/(Ts0 - Tc0) + Ts[τ]) s'[τ] +
1/(2 y) Ts'[τ]) (2 s[τ] -
E^(y s[τ]^2) Sqrt[π/y]
Erf[Sqrt[y] s[τ]]) == 4 y s[τ]^3 (Tc0/(Ts0 - Tc0) + Ts[τ]);
tempContinuityAdim = Ts[τ] == T[τ, 1];
heAdim = D[T[τ, x], τ] ==
1/(s[τ] - 1)^2 D[
T[τ, x], {x,
2}] + (2/(1 + x (s[τ] - 1)) 1/(s[τ] - 1) +
x/(s[τ] - 1) D[s[τ], τ]) D[T[τ, x], x];
1/(1 - s[τ]) D[T[τ, x], x] == rc/(k (Ts0 - Tc0)) qcmb[(rho cp rc^2)/k τ + t0] /. x -> 0;
1/(s[τ] - 1) D[T[τ, x], x] ==
-2 y s[τ] (Tc0/(Ts0 - Tc0) + Ts[τ]) /. x -> 1;

(*initial conditions*)
T[0, x] == (Exp[-y (2 x (s0/rc - 1) + x^2 (s0/rc - 1)^2)] - 1)/(
Ts0/Tc0 - 1);

(*parameters for transforming PDE's in ODE's*)
nbrPoints = 100;
scalingFactor = 1000;
xDiffOrder = 2;

{xL, xR} = domain = {0, 1};
grid = Array[# &, nbrPoints, domain];

ptoo = pdetoode[T[τ, x], τ, grid, xDiffOrder];
toode[expr_Equal] :=
With[{sf = scalingFactor}, sf # + D[#, τ] & /@ expr];
toode[expr_List] := toode /@ expr;

heODE = ptoo[heAdim]; del = #[[2 ;; -2]] &; heODE = del[heODE];

(*rewritting the equations like 'd.../dτ = ...'*)
dTdtau[τ] =
Flatten@{ptoo@Derivative[1, 0][T][τ, x], Ts'[τ],
s'[τ]};
solveDerivative =
Solve[Flatten@{Collect[heODE, ptoo@T[τ, x]],
Collect[bc1ODE,
Flatten[{ptoo@T[τ, x],
ptoo@Derivative[1, 0][T][τ, x]}]],
Collect[bc2ODE // Simplify,
Flatten[{ptoo@T[τ, x],
ptoo@Derivative[1, 0][T][τ, x]}]], tempContinuityODE,


Solving the equations using the default method gives:

result =
dTdtau[τ] == (dTdtau[τ] /. solveDerivative[])],
iniTODE // Simplify, iniTsODE, inisODE}, {T /@ grid, Ts, s} //
Flatten, {τ, 0,
k/(rho cp rc^2) (10^9 365.25 24 3600 4.5 - t0)}];
(*NDSolveValue::mxst : Maximum number of 10000 steps reached at the point τ == 3.640731908397398*^-12.*)


Using the StiffnessSwitching method, I am going a bit further and the error is different but we are still far from the end value $$\tau_{end} = 0.22086$$:

result =
dTdtau[τ] == (dTdtau[τ] /. solveDerivative[])],
iniTODE // Simplify, iniTsODE, inisODE}, {T /@ grid, Ts, s} //
Flatten, {τ, 0,
k/(rho cp rc^2) (10^9 365.25 24 3600 4.5 - t0)},
Method -> "StiffnessSwitching"];
(*NDSolveValue::ndsz: At τ == 6.36963610146291*^-11, step size is effectively zero; singularity or stiff system suspected.*)


Changing the number of points in the space grid (nbrPoints) or the scaling factor (scalingFactor) does not help: results are not converging:

And for a number of points larger than 250, I got the message error:

Cannot solve to find an explicit formula for the derivatives. Consider using the option Method->{\"EquationSimplification\"->\"Residual\"}.

Edit: the model

• heat equation in the conductive layer: $$\rho C_p \frac{\partial T}{\partial t}(t,r) = \frac{1}{r^2} \frac{\partial}{\partial r}\left(r^2 k \frac{\partial T}{\partial r}(t,r)\right)$$

• energy budget at the interface $$s(t)$$: cooling of the convective layer equals the heat conducted through the interface: $$-\rho C_p \int_0^{s(t)} \frac{\partial T_a}{\partial t}(t,r) \mathrm{d}V = -4\pi s^2(t) k \frac{\partial T_a}{\partial r}(t,s(t))$$

• given heat flux at the top of the core $$r_c$$: $$-k \frac{\partial T}{\partial r}(t,r_c) = q_c(t)$$

• temperature and heat flux continuity at the interface: $$T_s(t) = T(t,s(t))$$ and $$-k \frac{\partial T}{\partial r}(t,s(t)) = -k \frac{\partial T_a}{\partial r}(t,s(t))$$

with $$\rho$$, $$C_p$$, $$\alpha_c$$, $$g_c$$ and $$k$$ constants (density, specific heat, thermal expansivity, gravity and thermal conductivity respectively), $$T(t,r)$$ the temperature profile in the conductive layer, $$T_a(t,r)$$ the temperature profile in the convective layer, $$T_s(t) = T_a(t,s(t))$$ the temperature at the interface $$s(t)$$ and $$\frac{\partial T_a}{\partial r}(t,r) = -\frac{\alpha_c g_c}{r_c C_p}r T_a(t,r)$$ the adiabatic gradient.

In order to simplify the equations, I have considered the following variable changes: $$\tau = \frac{k}{\rho C_p r_c^2}(t-t_0)$$

$$x = \frac{r}{s(t)} \mathrm{ if} r \leq s(t) \mathrm{and} \frac{r-r_c}{s(t)-r_c} \mathrm{if} r\geq s(t)$$

$$T(\tau,x) = \frac{T(t,r) - T_{c0}}{T_{s0}-T_{c0}}$$

$$s(\tau) = \frac{s(t)}{r_c}$$ with $$T_{c0}$$, $$T_{s0}$$ and $$t_0$$ constants.

Edit 2: version 8

Using version 8 of mathematica, I can solve the equations but I do not have spatial convergence (example with scaling factor = 1, nx = nbrPoints): • Are you sure the underlying model itself is correct? Since you've mentioned "moving interface", I guess this problem is somewhat related to this one, right? Then I guess $s(\tau)$ is probably between $0$ and $1$? Nevertheless, with nbrPoints = 25; scalingFactor = 1; and MaxSteps -> 10 10^5 option I find $s(\tau)$ hits zero at about t = 0.005882695744315565. BTW the Simplify in NDSolveValue can be taken away, it only slows down the code. – xzczd Feb 27 '19 at 11:53
• @xzczd I'm quite sure of the model. And yes the problem is related to the one you mentioned. That's where I have found the useful tools ptoo and ptoode. The interface moves between 0 and 1. But I just tried with your modifications and I still have an error: "NDSolveValue::ndsz: At [Tau] == 1.1510886780232115*^-9, step size is effectively zero; singularity or stiff system suspected." Have you changed something else? – Mariel Feb 27 '19 at 13:31
• No, I don't make other modification. This seems to be another backslide of v11. I can reproduce the issue in v11.2, but my previous test is done in v9.0.1. So far I haven't found a way to adjust v11.2 to produce the result of v9.0.1. – xzczd Feb 27 '19 at 14:28
• The heat equation, and the other equations, are in spherical coordinates, explaining why it is $r^2$. I've done the variable change manually. I agree it's not the safest way but I did it so many times that I should have eliminate all the errors now. – Mariel Mar 1 '19 at 11:43
• How do you transform the integro-differential equation to pure PDE? – xzczd Mar 1 '19 at 13:35

This system of equations can be solved by explicit Euler in time.

(*constants*)rc = 2050 10^3;
cp = 850;
rho = 7200;
alpha = 5 10^-5;
gc = 4;
k = 40;
y = (gc alpha rc)/(2 cp);

(*parameters*)
s0 = 2049 10^3;
Tc0 = 2100;
Ts0 = Exp[(-alpha gc)/(2 cp rc) (s0^2 - rc^2)] Tc0;
t0 = 0.099 10^9 365.25 24 3600;
tm = k/(rho cp rc^2) (10^9 365.25 24 3600 4.5 - t0);

qcmb[t_] :=
With[{a = 0.004891658583550395, b = 0.34057028569554804,
c = 1.0021984665846737*^-15}, a + b E^(-c t)];

T[x_] := (Exp[-y (2 x (s0/rc - 1) + x^2 (s0/rc - 1)^2)] -
1)/(Ts0/Tc0 - 1);
T[-1][x_] := (Exp[-y (2 x (s0/rc - 1) + x^2 (s0/rc - 1)^2)] -
1)/(Ts0/Tc0 - 1);
s = s0/rc;
n = 200; tn = tm/4000;
Do[s[i] =
s[i - 1] +
tn*(4 y s[
i - 1]^3 (Tc0/(Ts0 - Tc0) + T[i - 1])/(2 s[i - 1] -
E^(y s[i - 1]^2) Sqrt[\[Pi]/y] Erf[Sqrt[y] s[i - 1]]) -
1/(2 y) (T[i - 1] - T[i - 2])/tn)/(s[
i - 1]*(Tc0/(Ts0 - Tc0) + T[i - 1])); np = i;
If[s[i] <= 0, Break[]];
T[i] = NDSolveValue[{(Ti[x] - T[i - 1][x])/tn ==
1/(s[i - 1] - 1)^2 *
Ti''[x] + (2/(1 + x (s[i - 1] - 1)) 1/(s[i - 1] - 1) +
x/(s[i - 1] -
1) ((4 y s[
i - 1]^3 (Tc0/(Ts0 - Tc0) +
T[i - 1])/(2 s[i - 1] -
E^(y s[i - 1]^2) Sqrt[\[Pi]/y] Erf[
Sqrt[y] s[i - 1]]) -
1/(2 y) (T[i - 1] - T[i - 2])/tn)/(s[
i - 1]*(Tc0/(Ts0 - Tc0) + T[i - 1])))) Ti'[x],
1/(1 - s[i - 1]) Ti' ==
rc/(k (Ts0 - Tc0)) qcmb[(rho cp rc^2)/k tn*i + t0],
1/(s[i - 1] - 1) Ti' == -2 y s[
i - 1] (Tc0/(Ts0 - Tc0) + T[i - 1])}, Ti, {x, 0, 1}];, {i,
1, n}] // Quiet

T3 = Table[{tn*i, x, T[i][x]}, {i, 0, np - 1}, {x, 0, 1, .02}];

T2 = Interpolation[Flatten[T3, 1]];

{Plot3D[T2[t, x], {t, 0, tn*(np - 1)}, {x, 0, 1}, Mesh -> None,
ColorFunction -> Hue, AxesLabel -> {"t", "x", ""}, PlotLabel -> "T",
PlotRange -> All],
ListLinePlot[Table[{tn*i, s[i]}, {i, 0, np - 1}], PlotRange -> All,
AxesLabel -> {"t", "s"}]} • I'm sorry, but the solution actually hasn't converged. If we choose tn = tm/8000, we'll see the calculation stops at about τ=0.0012, and with tn = tm/16000, the calculation stops at about τ=0.0008, with tn = tm/32000, stops at about τ=0.0005 – xzczd Feb 28 '19 at 5:45
• This is only the first step of the algorithm. In the scheme we need to add a variable step and increase the accuracy of the calculation Ti[x]. – Alex Trounev Feb 28 '19 at 11:59
• Well, the ODE solver of NDSolve is already adaptive (and long-tested), but as shown in OP's question, it just fails to go to τ=0.22. I really doubt if a self-made adaptive ODE solver will behave better. Though OP claims the underlying model is correct, I still suspect something is wrong with the equation system itself. – xzczd Feb 28 '19 at 12:18
• I agree that in the Stefan problem, the equation relates $ds/dt$ and heat flow from both sides proportionally$\lambda \nabla T$. In this case, we see the connection $ds/dt$ and $\partial T/\partial t$ at the mobile border. – Alex Trounev Feb 28 '19 at 13:27
• Thanks for your answer. I think there's an error in the discretisation of ds/dt: it should be 1/(2 y (1 - 2 s[i - 1])) and not 1/(2y)`. With this correction the result is less beautiful and there is still no convergence. – Mariel Feb 28 '19 at 15:21