3
$\begingroup$

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!

Improve this question
$\endgroup$
7
  • $\begingroup$ 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? $\endgroup$
    – xzczd ♦
    Commented Mar 23, 2023 at 3:07
  • $\begingroup$ 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 ... :/ $\endgroup$
    – Josè
    Commented Mar 23, 2023 at 9:32
  • $\begingroup$ 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. $\endgroup$
    – xzczd ♦
    Commented Mar 23, 2023 at 9:43
  • $\begingroup$ Another related post: mathematica.stackexchange.com/a/269866/1871 $\endgroup$
    – xzczd ♦
    Commented Mar 23, 2023 at 9:52
  • $\begingroup$ 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$? $\endgroup$
    – Alex Trounev
    Commented Mar 23, 2023 at 13:01

2 Answers 2

Reset to default
4
$\begingroup$

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}]}

Figure 1

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]}]

Figure 2

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 = 
 NDSolve`FiniteDifferenceDerivative[Derivative[2], xgrid, 
   DifferenceOrder -> 2]@"DifferentiationMatrix"; M1 = 
 NDSolve`FiniteDifferenceDerivative[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]}

Figure 3

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 = 
 NDSolve`FiniteDifferenceDerivative[Derivative[2], xgrid, 
   DifferenceOrder -> 2]@"DifferentiationMatrix"; M1 = 
 NDSolve`FiniteDifferenceDerivative[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]}

Figure 4

Improve this answer
$\endgroup$
11
  • 1
    $\begingroup$ 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! $\endgroup$
    – Josè
    Commented Mar 24, 2023 at 15:51
  • $\begingroup$ @Josè You are welcome! $\endgroup$
    – Alex Trounev
    Commented Mar 24, 2023 at 15:56
  • $\begingroup$ 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 ... $\endgroup$
    – Josè
    Commented Mar 25, 2023 at 12:26
  • $\begingroup$ Ah, I see that last part has no definition for S, u (done). Also, please run second part on the fresh kernel. $\endgroup$
    – Alex Trounev
    Commented Mar 25, 2023 at 13:28
  • 1
    $\begingroup$ @Josè Please, see also FDM solution as a control test. $\endgroup$
    – Alex Trounev
    Commented Mar 26, 2023 at 11:26
2
$\begingroup$

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,
   please find it in the link above. *)       
{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,
     please find it in the link above. *)
  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]

enter image description here

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

enter image description here

Improve this answer
$\endgroup$
1
  • $\begingroup$ Thank you. It is nice application of your perfect code pdetoode (+1). $\endgroup$
    – Alex Trounev
    Commented Apr 24, 2023 at 14:22

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.