# How to solve a system of PDEs with zip condition?

I have a system of PDEs in the following form,

$$\frac{\partial T_1}{\partial t}=\frac{\partial^2 T_1}{\partial x^2},\,\,0<x<S(t)$$ with $T_1(x,0)=-10,\,\,T_1(0,t)=-10,\,\, T_1(S(t),t)=10$, and $$\frac{\partial T_2}{\partial t}+\frac{dS(t)}{dt}\frac{\partial T_2}{\partial x}=\frac{\partial^2 T_2}{\partial x^2},\,\,x>S(t)$$ with $T_2(x,0)=10,\,\,T_2(S(t),t)=10,\,\, T_2(20,t)=20$.

The above two equations are zip together by $$\frac{dS(t)}{dt}=\frac{\partial T_1}{\partial x}\Biggl|_{x=S(t)}-\frac{\partial T_2}{\partial x}\Biggl|_{x=S(t)}$$ with $S(0)=0$.

Eq1 = D[T1[x, t], t] == D[T1[x, t], x, x]
Eq2 = D[T2[x, t], t]+D[S[t], t] D[T2[x, t], x]== D[T2[x, t], x, x]
T1[x, 0] == -10,T1[0,t] == -10, T1[S(t), t] == 10,T2[x, 0] == 10,T2[S(t), t] == 10,T2[20, t] == 20

Is there a way to solve the above system in Mathematica?

This problem is quite similar to PDE with Stefan Conditions, a.k.a variable boundary, and can be solved more easily because now we have DChange and pdetoode.

We first use DChange to transform $$T_1(x_1,t)$$ and $$T_2(x_2,t)$$ to $$T_1(\xi,t)$$ and $$T_2(\xi,t)$$, where

$$\frac{x_1}{S(t)}=\xi$$ $$\frac{x_2-S(t)}{20-S(t)}=\xi$$

xR = 20;
eqL = D[T1[x, t], t] == D[T1[x, t], x, x];
eqR = D[T2[x, t], t] + D[S[t], t] D[T2[x, t], x] == D[T2[x, t], x, x];
{icL, icR} = {T1[x, 0] == -10, T2[x, 0] == 10};
{bcL, bcR} = {{T1[0, t] == -10, T1[S@(t), t] == 10},
{T2[S@(t), t] == 10, T2[xR, t] == 20}};
slopeL = (D[T1[x, t], x] /. x -> S[t]);
slopeR = (D[T2[x, t], x] /. x -> S[t]);
bcmidfunc = S'[t] == # - #2 &;
With[{eps = 10^-10}, icmid = S[0] == eps];
(* Definition of DChange isn't included in this code piece,
changeL = DChange[#, x/S@t == ξ, x, ξ, T1[x, t]] &;
changeR = DChange[#, (x - S[t])/(xR - S[t]) == ξ, x, ξ, T2[x, t]] &;

(* Definition of diffbc isn't included in this code piece,
toode = With[{sf = 100}, diffbc[t, sf]];

{neweqL, newicL, newbcL, newslopeL} =
changeL@{eqL, icL, toode@bcL, slopeL} /. S[0] -> S[t]
{neweqR, newicR, newbcR, newslopeR} =
changeR@{eqR, icR, toode@bcR, slopeR} /. S[0] -> S[t] // Simplify

Remark

1. The /. S[0] -> S[t] is added because currently DChange can't handle coefficients involving the independent variable very well.

2. Check this post if you don't understand why we need toode.

After the change of variable, $$T_1$$ and $$T_2$$ are both defined in $$[0,1]\times[0,t_{end}]$$, while $$S$$ is defined in $$[0,t_{end}]$$, so NDSolve still can't solve the equation set directly. We need to discretize the PDEs to a set of ODEs. I'll use pdetoode for the task:

points = 25;
xdifforder = 4;
{ξL, ξR} = domain = {0, 1};
grid = Array[# &, points, domain];

(* Definition of pdetoode isn't included in this code piece,
ptoo = pdetoode[{T1, T2}[ξ, t], t, grid, xdifforder];
del = #[[2 ;; -2]] &;

{odeL, odeicL, odebcL, odeslopeL, odeR, odeicR, odebcR, odeslopeR} =
MapAt[del, ptoo /@ {neweqL, newicL, newbcL, newslopeL, neweqR, newicR, newbcR,
newslopeR}, {{1}, {5}}];

odebcmid = bcmidfunc[odeslopeL, odeslopeR]
odeicmid = icmid;
tend = 40;
sollst = NDSolveValue[{odeL, odeR, odeicL, odeicR, odebcL, odebcR,
odebcmid, odeicmid}, {T1 /@ grid, T2 /@ grid, S}, {t, 0, tend}];

solL = rebuild[#, grid, -1]&@sollst[[1]];
solR = rebuild[#, grid, -1]&@sollst[[2]];
solS = sollst[[-1]];
sol = {x, t} \[Function]
Piecewise[{{solL[x/solS@t, t], ξL <= x/solS@t <= ξR}},
solR[(x - solS[t])/(xR - solS[t]), t]];

Manipulate[Plot[sol[x, t], {x, 0, xR}, PlotRange -> {-10, 20}], {t, 0, tend}]
(* Plot3D[sol[x, t], {x, 0, xR}, {t, 0, tend}, PlotRange -> {-10, 20},
PlotPoints -> 50] *)

• Thanks for your efforts. I am unable to run your code without error.
– zhk
Commented Nov 12, 2016 at 15:54
• @mmm Have you already execute the definition of DChange and pdetoode? Commented Nov 12, 2016 at 16:01
• Now it is working. Thanks
– zhk
Commented Nov 12, 2016 at 16:23
• Is it possible to plot S vs t?
– zhk
Commented Nov 12, 2016 at 16:25
• @mmm Sure, just run Plot[solS[t], {t, 0, tend}] Commented Nov 13, 2016 at 2:24