PDE with Stefan Conditions, a.k.a variable boundary

Question

I want to solve the one-dimensional one-phase Stefan problem, but I don't know how to make Mathematica understand the conditions.

If you are not familiar with what I'm asking please refer to this wikipedia article: http://en.wikipedia.org/wiki/Stefan_problem#Mathematical_formulation

This is what I have so far. Clearly, it doesn't work.

NDSolve[
 {D[u[x, t], t] == D[u[x, t], {x, 2}],
  (-D[u[x, t], x] /. x -> 0) == 1,
  u[s[t], t] == 0,
  D[s[t]] == (-D[u[x, t], x] /. x -> s[t]),
  u[x, 0] == 0,
  s[0] == 0
  },
{u, s}, {x, 0, s[t]}, {t, 0, 10}]

I hope there is someone out there with a magical code

I'm using Mathematica 10.

Thanks!

Answer 1

One can do it semi-automatically. Let us introduce a normalized variable $$ \xi = \frac{x}{s(t)}, \quad \xi \in [0,1] $$ and make a simple finite difference method over $\xi$.

The differential equation in new variables is

ClearAll[u, s, x, t, ξ]
D[u[x/s[t], t], t] == D[u[x/s[t], t], x, x] /. x -> ξ s[t]

If we divide the interval $[0,1]$ by $n$ subintervals we come to the following finite difference scheme

n = 100;
δξ = 1./n;

ClearAll[dv, t];
dv[v_List] := With[{s = First@v, u = Rest@v},
   With[{ds = u[[-1]]/(s δξ),
     ξ = N@Range[n - 1]/n,
     d1 = ListCorrelate[{-0.5, 0, 0.5}/δξ, #] &,
     d2 = ListCorrelate[{1, -2, 1}/δξ^2, #] &},
    Prepend[d2[#]/s^2 + ξ ds d1[#]/s &@Join[{u[[1]] + s δξ}, u, {0.}], ds]
    ]];
s0 = 0.001;
v0 = Prepend[ConstantArray[0., n - 1], s0];
sol = NDSolve[{v'[t] == dv[v[t]], v[0] == v0}, v, {t, 0, 1}][[1, 1, 2]];

Here v contains s (the first element) and u (the rest list).

It remains only to decompose the interpolation function sol and return to the initial variable x

Needs@"DifferentialEquations`InterpolatingFunctionAnatomy`";
values = InterpolatingFunctionValuesOnGrid@sol;
valu = Transpose@Join[{#[[2]] + δξ #[[1]]}, Rest@#, {0 #[[1]]}] &@
   Transpose@InterpolatingFunctionValuesOnGrid@sol;
vals = InterpolatingFunctionValuesOnGrid[sol][[All, 1]];
t = Flatten@InterpolatingFunctionGrid@sol;
ξ = Range[0., n]/n;
s = ListInterpolation[vals, t];
uξ = ListInterpolation[valu, {t, ξ}];
u = If[#2 < s[#], uξ[#, #2/s[#]], 0.] &;

Visualization of the result

Show[{DensityPlot[u[t, x], {t, 0, 1}, {x, 0, 1}, FrameLabel -> {"t", "x"}], 
   Plot[s[t], {t, 0, 1}, PlotStyle -> {Red, Dashed}]}]

Answer 2

On the original question of how to solve classic Stefan-type problems using NDSolve with its high-level syntax (referred to in the posts as 'automatically'). This can be done. Essentially, first reformulate the problem as an optimal stopping problem: see van Moerbeke, Rocky Mountain J. Math. 4(3), 539-578, 1974. This leads to problems very similar to early exercise problems in Mathematical Finance. To solve these types of problems automatically with NDSolve, you can employ the WhenEvent methodology. Details are found in my recent book, "Option Valuation under Stochastic Volatility II", including a classic Stefan problem example. In Chapter 9, I treat a Stefan-type problem that is similar to, but not identical to the original question, namely: $$ \begin{align} w_{\tau} &= \frac{1}{2} w_{xx}, \quad x \in (s(\tau),\infty), \\ w(x,0) &= k \, 1_{\{x \ge 0 \}}, \quad (k > 0), \\ w(s(\tau),\tau) &= 0, \\ w_x(s(\tau),\tau) &= -2 \dot{s}(\tau). \end{align} $$ Here $w(x,\tau)$ is the water temperature for a water-ice phase transition and $s(\tau) \le 0$ is the moving boundary, which starts at $s(0)=0$.

The associated optimal stopping problem, using van Moerbeke, is $$ \begin{align} u_{\tau} &= \frac{1}{2} u_{xx}, \quad x \in (s(\tau),\infty), \\ u(x,0) &= (k+1) x^2 \, 1_{\{x > 0 \}}, \\ u(s(\tau),\tau) &= \tau, \\ u_x(s(\tau),\tau) &= 0. \end{align} $$

My main solver routine for the second problem is:

solverStefan[T_,k_,nexoppsperyear_,initsoln_,xmin_,xmax_,npts_]:=
Module[{xgrid,Soln,h,dt,hsoln,ecnt=0},
xgrid = xmin+(xmax-xmin)N[Table[i/npts,{i,0,npts}]]; 
dt = N[1/nexoppsperyear];
Soln =
 h/.NDSolve[{D[h[x,t],t] == 0.5 D[h[x,t],{x,2}],
 h[x,0] == initsoln[x],h[xmin,t] == t, h[xmax,t] == (k+1)(xmax^2+t),
 WhenEvent[Mod[t,dt] == 0, ecnt++;
 hsoln[x1_] := Max[h[x1,t], t];
 h[x,t]->Outer[hsoln[#1]&,xgrid]]},{h},{x,xmin,xmax},{t,0,T}, (* MUST USE {t,0,T} *)
  MaxSteps->1000000, MaxStepFraction->Min[dt/(10 T),1/npts],
  Method->{"MethodOfLines",
"SpatialDiscretization"->{"TensorProductGrid","Coordinates"->xgrid}}][[1]]; 
Return[Soln]]

What is going on is that there is a reward for stopping early, namely $\tau$. The WhenEvent method interrupts the solver repeatedly, and checks to see at what spatial points the current pde solution falls below the reward for stopping. If the stopping reward is higher, it replaces the solution with that reward.

In finance language, the reward $\tau$ is the payoff to the option buyer on early exercise of the option. If the option is not exercised early, the option buyer receives $u(x,0)$, which is a quadratic on $x>0$.

I call the solver routine with

StefanBoundary[T_,k_,xmin_,xmax_,npts_,noppsperyear_,size_]:=
  Module[{x1,initsoln,soln,g,MyTimeValue,eefunc,arg,alpha,s,
      eeexact,p1,p2a,p2},

Off[NDSolve::eerri];

(* initial soln *)
initsoln[x1_]=If[x1>0, (k+1) x1^2,0];

soln = solverStefan[T,k,noppsperyear,initsoln,xmin,xmax,npts];

(* reward function *)
g[x1_,t1_]:= If[t1>0,t1,initsoln[x1]];

 MyTimeValue[x1_,t1_] := soln[x1,t1]-g[x1,t1]; 
 eefunc[x1_,t1_] := If[Chop[MyTimeValue[x1,t1]]>0,1,0];

 p1 = ContourPlot[eefunc[x1,T-t1],{t1,0,T},{x1,xmin,2},
 Contours->2,ContourShading->True,PlotPoints->100,PlotLabel->"s(t) (PDE)",
 ImageSize->size]; 

(* Exact boundary *)
arg[a_] := a CN[a]-k E^(-a^2/2)/Sqrt[2 Pi];
alpha = a/.FindRoot[arg[a],{a,1}];
s[t_] := -alpha Sqrt[t];
eeexact[x1_,t1_] := If[x1>s[t1],1,0];

 p2a = ContourPlot[eeexact[x1,T-t1],{t1,0,T},{x1,xmin,2},
 Contours->2,ContourShading->True,PlotPoints->100,
 ImageSize->size,PlotLabel->"s(t) (Exact)"]; 

p2 = Show[p2a,Graphics[Inset[Framed[Style["Ice",12],
Background->LightYellow],{T/2,-3.0}]],
Graphics[ Inset[Framed[Style["Water",12],
Background->LightYellow],{T/2,0}]]];

Print["Stefan:finished: MMU=",MMU," GB"];
Return[GraphicsRow[{p1,p2},Spacings->0]]]

Some auxiliary functions used are:

MMU: = N[MaxMemoryUsed[]/10^9];
Cumnormal[xx_] := (1+Erf[xx/Sqrt[2]])/2;
CN = Cumnormal; 

Running

StefanBoundary[5,3,-5,5,500,500,175]

yields the NDSolve boundary on the left and the exact on the right: