I came across the paper Solidification dynamics of an impacted drop regarding a heat equations by Thiévenaz et.al and was interested in knowing how they obtained the graphs presented. From what I understand, they solved the PDE presented below, with a moving boundary using the finite differences method. I did solve them in matlab but was wondering how one would do this in mathematica given that the flexibility of matlab with matrices is much greater. The code in matlab may be found here. The relevant equations are: The two heat equations for the temperature field $T(z,t)$ $$\frac{\partial T}{\partial t}=D_s\frac{\partial^2 T}{\partial z^2}, \,\,z\leq0;\qquad \frac{\partial T}{\partial t}=D_i\frac{\partial^2 T}{\partial z^2}, \,\,0\leq z\leq h(t)$$ and the four boundary conditions:

$$T(0^-,t)=T(0^+,t);\quad \lambda_s\frac{\partial T(0^-,t)}{\partial t}=\lambda_i\frac{\partial T(0^+,t)}{\partial t};\quad T(h(t)^-,t)=T_m;$$ $$\lambda_i\frac{\partial T(h(t)^-,t)}{\partial t}=\rho_i L\frac{dh}{dt}$$ where $L$ denotes the latent heat of solidification, $\rho_k$ the thermal density, $\lambda_n$ the thermal conductivity, $D_k=\lambda_k/(\rho_kC_k)$ the diffusion coefficient. Here $k$ denotes $l$ in the liquid phase, $i$ in the solid phase and $s$ for the substrate (they are all constant parameters). It is also imposed a constant temperature $T_s$ deep in the substrate, i.e., $$\lim_{z\to-\infty}T(z,t)=T_s,\quad T(z,0)=T_s\quad z\leq0, h(0)=0$$ which means that at time $0$ the liquid is put into contact with the substrate.

A pictorial description of the model is showed below: enter image description here

Note: The constants mentioned above in my code are chosen arbitrarily (i.e. I chose the values assigned to them)

I did not post the code because it was somewhat long.

As requested in the comments, the scheme I used was:

$$\vec{T^{\,n+1}}=A^{-1}\left(\vec{T^{\,n}}+\vec{B}\right)\implies c\rho\frac{T_i^{n+1}-T_i^n}{\Delta t}=k\frac{T_{i-1}^{n+1}-2T_i^{n+1}+T_{i+1}^n}{\Delta x^2}$$ which is the implicit version. The explicit is similar, and looks like $$c\rho\frac{T_i^{n+1}-T_i^n}{\Delta t}=k\frac{T_{i-1}^{n}-2T_i^{n}+T_{i+1}^n}{\Delta x^2}$$ In my code you will see how I found $A$ and $B$. ($A$ is a matrix and $B$ a vector)

  • 1
    $\begingroup$ I'd say Mathematica and MATLAB are about the same when it comes to 95% of the major matrix ops, but maybe you mean element-wise manipulation? Also I think it'd be nice to at least post some code detailing the method used. It's a big ask for us to read, translate, and do all the mental processing without any sort of starting point. $\endgroup$
    – b3m2a1
    Commented Dec 10, 2019 at 6:03
  • $\begingroup$ @b3m2a1 Yes what I meant was element-wise manipulation. And it appears to me it feels more natural in matlab when doing so, or maybe I’m not experienced enough in mathematica. Regarding the method, as I mentioned above, I used the finite difference method which there are a lot of snippets around the web. I will try to post some code, maybe the general idea. Anyways, any help is appreciated. It’s just a curiosity $\endgroup$
    – DMH16
    Commented Dec 10, 2019 at 6:06
  • 2
    $\begingroup$ I'd say most element-wise tasks can be done better as vectorized operations and once you do that there's really no difference between Mathematica and MATLAB. That's really neither here nor there, though. Looks like you just need to use the FDM, from what I glanced at in that image. Lots of questions here pertain to that but maybe xczd will come by and help you. Also how about as a compromise for not having some preliminary code, you post the relevant equations here in LaTeX form? $\endgroup$
    – b3m2a1
    Commented Dec 10, 2019 at 6:09
  • 3
    $\begingroup$ Finally, we already have a few posts about moving b.c., have you read them?: 1. mathematica.stackexchange.com/q/58596/1871 2. mathematica.stackexchange.com/q/130962/1871 3. mathematica.stackexchange.com/q/95358/1871 4. mathematica.stackexchange.com/q/114967/1871 5. mathematica.stackexchange.com/q/184920/1871 $\endgroup$
    – xzczd
    Commented Dec 10, 2019 at 6:35
  • 1
    $\begingroup$ The question is much better after the edits. (Still, you need to show us the parameters and i.c.s. ) I'd say the MATLAB code isn't quite related to this question in my view. Is the finite difference scheme used by the paper necessary? If so, it'll be better to add the difference scheme to the question directly, because it's easier for us to understand the difference formula in traditional math notation. If not, then I believe the combination of DChange, pdetoode and NDSolve is a better choice for solving the problem, and again, we don't need to know how the MATLAB code looks like. $\endgroup$
    – xzczd
    Commented Dec 10, 2019 at 6:58

1 Answer 1


As pointed out in the comment, OP's MATLAB implementation doesn't seem to be correct, so I'll just give a solution to OP's problem.

Since the heat flux continuity involves in, there's indeed something new in this problem compared to previous ones. We first handle the continuity condition. Notice that, strictly speaking, 1D heat conduction equation is not

$$\frac{\partial T}{\partial t}=\frac{\lambda}{\rho C}\frac{\partial^2 T}{\partial z^2}$$


$${\rho C}\frac{\partial T}{\partial t}=\frac{\partial}{\partial z}\left(\lambda\frac{\partial T}{\partial z}\right)$$

These 2 equations are equivalent only if $\lambda$ is constant, and piecewise constant is not constant. "So what? Why you mention this? " Because the latter equation owns a good property: when it's discretized directly, the heat continuity condition will be automatically satisfied.

Then we need to deal with the moving b.c.. As usual, we use the "change the variable -> discretize" method, with the help of pdetoode.

Given that OP hasn't shown the parameters, they'll be casually chosen as

ρL = 1; eps = 10^-2;
inf = 2; Tm = 1; Ts = 2; tend = 10;
ρc = Piecewise[{{2, z > 0}}, 1];
λ = Piecewise[{{3, z > 0}}, 1];
icexpr = Piecewise[{{Ts, z > 0}}, Tm];

Next we interpret the PDE, i.c. and b.c. to Mathematica code:

With[{T = T[z, t], dT = dT[z, t]}, 
 eq = ρc D[T, t] == D[dT, z];
 eqd = dT == λ D[T, z];
 bc = {{λ D[T, z] == ρL h'[t], T == Ts} /. z -> h[t], 
       T == Tm /. z -> -inf};
 icextra = h[0] == eps;
 ic = T == icexpr /. t -> 0;]

Then we change the variable with DChange:

unitStepExpand = Simplify`PWToUnitStep@PiecewiseExpand@# &;
(* Definition of pdetoode isn't included in this post,
   please find it in the link above. *)
{neweqd, neweq, newicmid, newbcmid} = 
  DChange[{eqd, eq, ic, bc}, ζ == scale z, z, ζ, {T[z, t], dT[z, t]}] /. 
    a_Piecewise :> unitStepExpand@a /. 
   scale -> unitStepExpand@Piecewise[{{1/h[t], ζ > 0}}, 1];

newic = newicmid /. t -> 0 /. Rule @@ icextra

newbc = {newbcmid[[1]] /. ζ -> 1, newbcmid[[2]] /. ζ -> -inf}

Discretize the system:

domain@1 = {-inf, 0};
domain@2 = {0, 1};
points = 400;
difforder = 2;
gridfunc = Array[# &, points, #] &;
grid = gridfunc@domain@1~Join~Rest@gridfunc@domain@2;

(*Definition of pdetoode isn't included in this post,
  please find it in the link above.*)
ptoofunc = pdetoode[{T, dT}[ζ, t], t, grid, difforder];
del = #[[2 ;; -2]] &;
With[{eq = neweq, eqd = neweqd, ic = newic, bc = newbc},
 ode = Block[{dT}, Set @@ ptoofunc /@ eqd;
 odeic = ptoofunc@ic;
 odebc = ptoofunc@bc;]    

The resulting system is a DAE system and NDSolve can't handle it well, so we create a toode function to transform the system to an ODE system (For more info check this and this post) and solve:

toode = With[{sf = 10}, diffbc[t, sf]];
{funch, funcTlist} = 
   MapAt[toode, {ode, icextra, odeic, odebc} /. T[i_] :> T@ToString@i // 
     Flatten, {{-1}, {-2}}], {h, T@ToString@# & /@ grid}, {t, 0, tend}, 
   MaxSteps -> Infinity];

funcT[1] = rebuild[funcTlist[[;; points]], grid[[;; points]], 2];
funcT[2] = rebuild[funcTlist[[-points ;;]], grid[[-points ;;]], 2];
func = Function[{z, t}, 
   Piecewise[{{funcT[1][z, t], z <= 0}, {Ts, z >= funch[t]}}, 
    funcT[2][z/funch[t], t]]];

Finally, just some illustration:

ListLinePlot[funch, PlotRange -> All]

enter image description here

Here funch is plotted using an undocumented syntax.

Manipulate[Plot[func[z, t], {z, -inf, funch[tend]}, PlotRange -> {0, 2}, 
  AxesLabel -> {"z", "T"}], {t, 0, tend}]

enter image description here

Animate[DensityPlot[func[z, t], {x, 0, 1}, {z, -inf, funch[tend]}, PlotRange -> {0, 2}, 
  PlotPoints -> 50, PerformanceGoal -> "Quality", ColorFunction -> "DeepSeaColors"], {t, 
  0., tend, 0.05}]

enter image description here

  • 2
    $\begingroup$ The process is very clear. I will have to look up the definitions to understand the code. Anyways, thank you for putting the time $\endgroup$
    – DMH16
    Commented Dec 13, 2019 at 23:23
  • $\begingroup$ @DMH16 I just noticed the moving grid method is wrong. See the corrected answer. (The difference is not obvious in this case, but it does exist. ) $\endgroup$
    – xzczd
    Commented Apr 24, 2023 at 13:08

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.