# Stefan problem with mixed bc

I am trying to solve through Mathematica the classical Stefan problem $$\left\{ \begin{array}{lll} \dot{v}(x,t)=v_{xx}(x,t) & x\in(0,s(t))\\ \dot{s}(t)=-v_x(s(t),t) & x = s(t)\\ v(0,t) = 0 & x = 0\\ v(x,0) = v_0(x) & \\ \end{array} \right.$$ with the novelty that the boundary condition on the moving front is not the usual $$v(s(t),t)=0$$, but rather $$\displaystyle\color{red}{v_x(s(t),t) - \frac{1}{s(t)} v(s(t),t) = F(t)}$$ with $$F(t)$$ some given function. I have tried to understand how the great answer of @ybeltukov could be modified, even (say) in the simple case of $$F(t)=t, \qquad v_0(x)=x$$ But I have been scratching my head uselessly for days ... I hope somebody could give me a grind on this!

Thanks so much!

• There're many other related posts in this site, e.g. mathematica.stackexchange.com/q/211080/1871 (Don't miss the links in comment), have you read them? Commented Mar 23, 2023 at 3:07
• Hi @xzczd, I think I've spent the whole week to read the various posts and suggested method of implementation of Stefan problem; nobody considers general types of bc on the moving front. Which is crucial for my research ... :/
– Josè
Commented Mar 23, 2023 at 9:32
• I don't think there's anything special in your b.c., it's merely a Robin b.c. involving $s(t)$ term. Please think carefully about change of variable, and play with DChange/DSolveChangeVariables. Commented Mar 23, 2023 at 9:43
• Another related post: mathematica.stackexchange.com/a/269866/1871 Commented Mar 23, 2023 at 9:52
• It is not clear what initial condition $v_0(x)=x$ means in this case? If $0\le x\le s(t)$ and $s(0)=0$ then $v_0(x)=0$. Should we suppose that $s(0)>0$? Commented Mar 23, 2023 at 13:01

This problem can be solved with using the Euler wavelets collocation method as follows. We introduce same normalized variable as here in the form $$x\rightarrow \frac {x}{s(t)}$$. First, we test code utilizing numerical solution from @ybeltukov answer

OEm[m_, x_] :=
Sqrt[2 m +
1] Sum[(-1)^(m - k) x^k Binomial[m, k] Binomial[m + k, k], {k, 0,
m}]; UE[m_, t_] := OEm[m, t];
psi[k_, n_, m_, t_] :=
Piecewise[{{2^((k - 1)/2) UE[m, 2^(k - 1) t - n + 1], (n - 1)/
2^(k - 1) <= t < n/2^(k - 1)}, {0, True}}];
PsiE[k_, M_, t_] :=
Flatten[Table[psi[k, n, m, t], {n, 1, 2^(k - 1)}, {m, 0, M - 1}]]
k0 = 3; M0 = 4;
With[{k = k0, M = M0},
nn = Length[Flatten[Table[1, {n, 1, 2^(k - 1)}, {m, 0, M - 1}]]]];
dx = 1/(nn); xl = Table[l*dx, {l, 0, nn}]; xcol =
Table[(xl[[l - 1]] + xl[[l]])/2, {l, 2, nn + 1}];
Psijk = With[{k = k0, M = M0}, PsiE[k, M, t1]]; Int1 =
With[{k = k0, M = M0}, Integrate[PsiE[k, M, t1], t1]];
Int2 = Integrate[Int1, t1];
Psi[y_] := Psijk /. t1 -> y;
int1[y_] := Int1 /. t1 -> y;
int2[y_] := Int2 /. t1 -> y;
wA = Table[wa[i][t], {i, nn}]; wB = Table[wb[i][t], {i, 2}];
v2[x_] := wA . Psi[x]; v1[x_] := wA . int1[x] + wB[[1]];
v0[x_] := wA . int2[x] + wB[[1]] x + wB[[2]];
equ = {D[v0[x], t] + (x v1[1] v1[x])/s[t]^2 == v2[x]/
s[t]^2}; eqnS = {s'[t] == -v1[1]/s[t]};
eqx = Table[equ, {x, xcol}];
eqs = Join[Flatten[eqx], eqnS];

bc = Join[{v1[0] + s[t] == 0}, {v0[1] == 0}]; icx = {v0[x] == 0 /.
t -> 0}; ic = Table[icx, {x, xcol}] // Flatten;
varAll = Join[wA, wB, {s[t]}];
icn = Join[ic, bc /. t -> 0, {s[0] == 10^-3}]; eqnN =
Join[eqs, D[bc, t]]; var1 = D[varAll, t];
{vec, mat} = CoefficientArrays[eqnN, var1];

f = Inverse[mat // N] . (-vec);

vr0 = varAll /. t -> 0; {w0, mat0} = CoefficientArrays[icn, vr0];

s0 = Inverse[mat0 // N] . (-w0); rul0 =
Table[vr0[[i]] -> s0[[i]], {i, Length[vr0]}];
f0 = f /. t -> 0 /. rul0; m = Length[f];

sol = NDSolve[{Table[var1[[i]] == f[[i]], {i, Length[var1]}],
Table[vr0[[i]] == s0[[i]], {i, Length[vr0]}]}, varAll, {t, 0, 1}];


Visualization and comparison with @ybeltukov solution (dashed lines)

lst = Table[{{t, x}, v0[x] /. sol[[1]]}, {t, 0, 1, .1}, {x, 0,
1, .1}]; u = Interpolation[Flatten[lst, 1]]; S =
Interpolation[Table[{t, s[t] /. sol[[1]]}, {t, 0, 1, .01}]];
{Show[{DensityPlot[
u[t, x/S[t]] UnitStep[S[t] - x], {t, 0, 1}, {x, 0, 1},
FrameLabel -> {"t", "x"}, ColorFunction -> "Rainbow"],
Plot[s[t] /. sol[[1]], {t, 0, 1}, PlotStyle -> {Red, Dashed}]}],Plot[{u[t, 0], u[1, t/S[1]] UnitStep[S[1] - t], S[t]}, {t, 0, 1},
PlotLegends -> Automatic, Frame -> True,
PlotStyle -> {Red, Blue, Black}]}


Note, that two numerical solution are in a good agreement, therefore we can try to solve problem proposed by Josè. In this case we consider solution with initial data $$v(x,0)=vini(x)=x, s(0)=3$$, we have

s0 = 3; F[t_] := t; vini[x_] := x; tmax = 1.;

OEm[m_, x_] :=
Sqrt[2 m +
1] Sum[(-1)^(m - k) x^k Binomial[m, k] Binomial[m + k, k], {k, 0,
m}]; UE[m_, t_] := OEm[m, t];
psi[k_, n_, m_, t_] :=
Piecewise[{{2^((k - 1)/2) UE[m, 2^(k - 1) t - n + 1], (n - 1)/
2^(k - 1) <= t < n/2^(k - 1)}, {0, True}}];
PsiE[k_, M_, t_] :=
Flatten[Table[psi[k, n, m, t], {n, 1, 2^(k - 1)}, {m, 0, M - 1}]]
k0 = 3; M0 = 4;
With[{k = k0, M = M0},
nn = Length[Flatten[Table[1, {n, 1, 2^(k - 1)}, {m, 0, M - 1}]]]];
dx = 1/(nn); xl = Table[l*dx, {l, 0, nn}]; xcol =
Table[(xl[[l - 1]] + xl[[l]])/2, {l, 2, nn + 1}];
Psijk = With[{k = k0, M = M0}, PsiE[k, M, t1]]; Int1 =
With[{k = k0, M = M0}, Integrate[PsiE[k, M, t1], t1]];
Int2 = Integrate[Int1, t1];
Psi[y_] := Psijk /. t1 -> y;
int1[y_] := Int1 /. t1 -> y;
int2[y_] := Int2 /. t1 -> y;
wA = Table[wa[i][t], {i, nn}]; wB = Table[wb[i][t], {i, 2}];
v2[x_] := wA . Psi[x]; v1[x_] := wA . int1[x] + wB[[1]];
v0[x_] := wA . int2[x] + wB[[1]] x + wB[[2]];
equ = {D[v0[x], t] + (x v1[1] v1[x])/s[t]^2 ==
v2[x]/s[t]^2}; eqnS = {s'[t] == -v1[1]/s[t]};
eqx = Table[equ, {x, xcol}];
eqs = Join[Flatten[eqx], eqnS];

bc = Join[{v0[0] == 0}, {v1[1] - v0[1] == s[t] F[t]}]; icx = {v0[x] ==
vini[x] s0 /. t -> 0}; ic = Table[icx, {x, xcol}] // Flatten;
varAll = Join[wA, wB, {s[t]}];
icn = Join[ic, bc /. t -> 0, {s[0] == s0}]; eqnN =
Join[eqs, D[bc, t]]; var1 = D[varAll, t];

f = var1 /. Solve[eqnN, var1][[1]];

vr0 = varAll /. t -> 0; s00 = vr0 /. Solve[icn, vr0][[1]];

sol = NDSolve[{Table[var1[[i]] == f[[i]], {i, Length[var1]}],
Table[vr0[[i]] == s00[[i]], {i, Length[vr0]}]},
varAll, {t, 0, tmax}];


Visualization

S = Interpolation[
Table[{t, s[t] /. sol[[1]]}, {t, 0, tmax, .01}]]; lst =
Table[{{t, x}, v0[x] /. sol[[1]]}, {t, 0, tmax, .05}, {x, 0,
1, .05}];
Show[{DensityPlot[
u[t, x/S[t]] UnitStep[S[t] - x], {t, 0, tmax}, {x, 0, S[tmax]},
FrameLabel -> {"t", "x"}, ColorFunction -> "Rainbow",
PlotRange -> All, PlotLegends -> Automatic],
Plot[S[t], {t, 0, tmax}, PlotStyle -> {Red, Dashed},
PlotRange -> All]}]


Update 1. Second method to solve this problem is FDM, same as @ybeltukov, but much clear

ns = 101; h = 1/(ns - 1);
x[t_] = Table[Symbol["x" <> ToString[i]][t], {i, 1, ns}];
xgrid = Table[k h, {k, 0, ns - 1}]; M2 =
NDSolveFiniteDifferenceDerivative[Derivative[2], xgrid,
DifferenceOrder -> 2]@"DifferentiationMatrix"; M1 =
NDSolveFiniteDifferenceDerivative[Derivative[1], xgrid,
DifferenceOrder -> 2]@"DifferentiationMatrix";
eqnS = {s'[t] == -(M1 . x[t])[[ns]]/s[t]}; eq =
Table[(D[x[t], t])[[i]] +
xgrid[[i]] (M1 . x[t])[[ns]] (M1 . x[t])[[i]]/
s[t]^2 == (M2 . x[t])[[i]]/s[t]^2, {i, 2, ns - 1}];
bc = {(M1 . x[t])[[1]] == -s[t], x[t][[ns]] == 0};
ini = Join[{s[0] == .001}, Table[x[0][[i]] == 0, {i, 2, ns - 1}]];

sol = NDSolve[{Join[eqnS, eq, D[bc, t]], Join[ini, bc /. t -> 0]},
Join[{s[t]}, x[t]], {t, 0, 1}];


Visualization and comparison with previous solution (dashed lines)

u = Interpolation[
Flatten[Table[{{t, xgrid[[i]]}, x[t][[i]] /. sol[[1]]}, {t, 0,
1, .1}, {i, ns}], 1]];

S = Interpolation[Table[{t, s[t] /. sol[[1]]}, {t, 0, 1, .001}]];

{Show[{DensityPlot[
u[t, x/S[t]] UnitStep[S[t] - x], {t, 0, 1}, {x, 0, 1},
FrameLabel -> {"t", "x"}, ColorFunction -> "Rainbow"],
Plot[S[t], {t, 0, 1}, PlotStyle -> {Red, Dashed}]}],Plot[{u[t, 0], u[1, t/S[1]] UnitStep[S[1] - t], S[t]}, {t, 0, 1},

Frame -> True, PlotStyle -> {Red, Green, Blue},
PlotLegends -> Automatic]}


Finally we can solve problem proposed by Josè with initial data $$v(x,0)=v_0(x)=x,s(0)=3$$, and for $$F(t)=t$$, we have

F[t_] := t; v0[x_] := x;

ns = 101; h = 1/(ns - 1);
x[t_] = Table[Symbol["x" <> ToString[i]][t], {i, 1, ns}];
xgrid = Table[k h, {k, 0, ns - 1}]; M2 =
NDSolveFiniteDifferenceDerivative[Derivative[2], xgrid,
DifferenceOrder -> 2]@"DifferentiationMatrix"; M1 =
NDSolveFiniteDifferenceDerivative[Derivative[1], xgrid,
DifferenceOrder -> 2]@"DifferentiationMatrix";
eqnS = {s'[t] == -(M1 . x[t])[[ns]]/s[t]}; eq =
Table[(D[x[t], t])[[i]] +
xgrid[[i]] (M1 . x[t])[[ns]] (M1 . x[t])[[i]]/
s[t]^2 == (M2 . x[t])[[i]]/s[t]^2, {i, 2, ns - 1}];
bc = {(M1 . x[t])[[ns]] - x[t][[ns]] == s[t] F[t], x[t][[1]] == 0};
ini = Join[{s[0] == 3},
Table[x[0][[i]] == 3 v0[xgrid[[i]]], {i, 2, ns - 1}]];

sol = NDSolve[{Join[eqnS, eq, D[bc, t]], Join[ini, bc /. t -> 0]},
Join[{s[t]}, x[t]], {t, 0, 1}];


Visualization and comparison with wavelets solution (dashed lines)

Show[{DensityPlot[
u[t, x/S[t]] UnitStep[S[t] - x], {t, 0, 1}, {x, 0, 3},
FrameLabel -> {"t", "x"}, ColorFunction -> "Rainbow"],
Plot[S[t], {t, 0, 1}, PlotStyle -> {Red, Dashed}]}];

Plot[{u[1, t/S[1]] UnitStep[S[t] - t], S[t]}, {t, 0, 1},
PlotStyle -> Dashed, PlotLegends -> Automatic];

u = Interpolation[
Flatten[Table[{{t, xgrid[[i]]}, x[t][[i]] /. sol[[1]]}, {t, 0,
1, .1}, {i, ns}], 1]];

S = Interpolation[Table[{t, s[t] /. sol[[1]]}, {t, 0, 1, .001}]];
{Show[{DensityPlot[
u[t, x/S[t]] UnitStep[S[t] - x], {t, 0, tmax}, {x, 0, s0},
FrameLabel -> {"t", "x"}, ColorFunction -> "Rainbow",
PlotRange -> All, PlotLegends -> Automatic],
Plot[S[t], {t, 0, tmax}, PlotStyle -> {Red, Dashed},
PlotRange -> All]}],Plot[{u[1, t/S[1]] UnitStep[S[t] - t], S[t]}, {t, 0, 1},
PlotStyle -> {Red, Blue}, PlotLegends -> Automatic]}


• But this is just amazing! you can not even imagine how grateful am I …. I would have NEVER arrived here! I will play with it and will give feedback on how it goes :) thanks heathfully!
– Josè
Commented Mar 24, 2023 at 15:51
• @Josè You are welcome! Commented Mar 24, 2023 at 15:56
• Ok so. The 1st part where you redo the solution of the Stefan problem is ok. 2nd part where you solve my problem doesn't run on my machine (Mathematica 12). Specifically the (first) problematic line is: f = Inverse[mat // N].(-vec); Delivering several: General::munfl: 1.493*10^-237 (-7.40643*10^-238) is too small to represent as a normalized machine number; precision may be lost. altogether with: Inverse::sing: Matrix {<<1>>} is singular. Any clue? Thanks so much ...
– Josè
Commented Mar 25, 2023 at 12:26
• Ah, I see that last part has no definition for S, u (done). Also, please run second part on the fresh kernel. Commented Mar 25, 2023 at 13:28
• @Josè Please, see also FDM solution as a control test. Commented Mar 26, 2023 at 11:26

Let me add my FDM solution. The spirit of this solution is the same as that of Alex's FDM solution, except that the process is automated with pdetoode and DChange. Since there's nothing new in the following code compared to the previous moving boundary problems in this site, I won't explain much. (Anyway, if you still feel it difficult to follow, feel free to ask in the comment, but please be specific. )

I've chosen the same initial data as in Alex's answer so we can easily verify the result.

s0 = 3;
F = t; v0 = x;
With[{v = v[x, t],
s = s[t]}, {eq, ic,
bc} = {{D[v, t] == D[v, x, x],
D[s, t] == -D[v, x] /. x -> s}, {v == v0, s == s0} /. t -> 0, {v == 0 /. x -> 0,
D[v, x] - 1/s v == F /. x -> s}}];
(* Definition of DChange isn't included in this post,
{neweq, newicmid, newbc} = DChange[{eq, ic, bc}, x/s[t] == ξ, x, ξ, v[x, t]]
newic = {newicmid[[1]] /. t -> 0 /. s[0] -> s0, newicmid[[2]]};

With[{ic = newic, eq = neweq, bc = newbc},
domain = {0, 1}; points = 25; grid = Array[# &, points, domain]; difforder = 2;
(* Definition of pdetoode isn't included in this post,
ptoofunc = pdetoode[v[ξ, t], t, grid, difforder];
del = #[[2 ;; -2]] &;

ode = {del@ptoofunc@eq[[1]], ptoofunc@eq[[2]]};
odebc = ptoofunc@bc;
odeic = {ic[[1]] // ptoofunc // del, ic[[2]]} /. t -> 0;
tend = 1;
{ssol, vsollst} =
NDSolveValue[{ode, odeic, odebc}, {s, v /@ grid}, {t, 0, tend}]; // AbsoluteTiming]

vscaledsol = rebuild[vsollst, grid, 2];
vsol = {x, t} |-> vscaledsol[x/ssol[t], t];

reg = DiscretizeRegion@ImplicitRegion[x <= ssol[t], {{t, 0, 1}, {x, 0, 3}}];

DensityPlot[vsol[x, t], {t, x} ∈ reg]


Manipulate[
Plot[vsol[x, t], {x, 0, ssol[t]}, PlotRange -> {{0, s0}, {-1, 5}}], {t, 0, 1}]


• Thank you. It is nice application of your perfect code pdetoode (+1). Commented Apr 24, 2023 at 14:22