# Finite differences for the one-phase Stefan problem

I am trying to code the one-phase, one-dimensional Stefan problem using finite differences in Matlab, similarly to what has already been done in Mathematica (see https://mathematica.stackexchange.com/questions/58596/pde-with-stefan-conditions-a-k-a-variable-boundary).

The problem consists in finding $$(u,s)$$ satisfying the coupled PDE-ODE system

$$\begin{equation} \begin{cases} u_t - u_{xx} = 0 &\text{for} \ 0 < t < T, \ 0 < x < s(t),\\ u(0, t) = f(t) &\text{for} \ 0 < t < T, \\ u(s(t),t) = 0, &\text{for} \ 0 < t < T,\\ \dot{s}(t) = -u_x(s(t),t) &\text{for} \ 0 < t < T,\\ u(x,0)=u_0(x), \ s(0)=s_0, \end{cases} \end{equation}$$ where $$T, u_0, f, s_0$$ are given (positive) data (to assure the positivity of the state $$u$$ in the set where it solves the heat equation). One assumes that $$u=0$$ elsewhere in space.

Conducting the change of variable $$\xi = \frac{x}{s(t)}, \ z(\xi,t) = u(x,t),$$ the governing equations hold in the fixed space domain $$(0, 1)$$: $$\begin{equation} \begin{cases} z_t - \frac{1}{s^2}z_{\xi\xi} - \frac{\xi \dot{s}}{s}z_\xi = 0 & (\xi, t) \in (0,1)\times(0,T), \\ \dot{s}(t) = -\frac{1}{s(t)}u_x(1,t) &t \in (0, T), \end{cases} \end{equation}$$ now with a homogeneous Dirichlet boundary condition at $$x=1$$.

I used an explicit Euler scheme for discretizing $$\dot{s}$$, and an implicit scheme for the time derivative of the state $$z$$, and centered differences for its space derivatives. I chose the initial and boundary data $$f \equiv 1, u_0 \equiv 0.1$$ and $$s_0 = 0.05$$.

The problems I encountered are the following:

• The solution seems to equal $$0$$ after a very short time, no matter the choice of the initial data;

• The interface $$s$$ appears to have a non-smooth behavior when $$u_0 > s_0$$;

I would be very grateful if someone could help me out with the above issues. I am looking to obtain a similar output as the answer in Mathematica SE.

The code is below:

clc;

nx = 75; nT = 25;
dt = 1/nT; dx = 1/nx;
M1 = 1; M2 = 0;
s0 = 0.1;

U = zeros(nx+1,nT+1);
U(1,1) = s0; U(2:nx+1,1) = 0.05;

Id = eye(nx);
A = -2*Id + diag(ones(nx-1,1),-1) + diag(ones(nx-1,1),1); A = (1/dx^2)*A;
Ma = zeros(nx,1); Ma(1,1) = M1; Ma(nx,1) = M2; Ma = Ma/(dx^2);
G = diag(ones(1,nx-1),1) + diag(-ones(1,nx-1),-1); G = 1/(2*dx)*G;
Mg = zeros(nx,1); Ma(1,1) = -M1; Ma(nx,1) = M2; Mg = Mg/(2*dx);
x = linspace(0, 1, nx+2)';

for j=1:nT
U(1,j+1) = (dt/U(1,j))*U(nx+1,j)/dx + U(1,j);
ds = (U(1,j+1)-U(1,j))/dt;
U(2:nx+1,j+1) = (U(1,j)^2*Id-dt*A-U(1,j)*ds*dt*x(2:nx+1).*G)\(dt*Ma+U(1,j)*ds+dt*x(2:nx+1).*Mg+U(2:nx+1,j));
end

s = U(1,1:nT+1);
tint = linspace(0,1,nT+2);
plot(tint(2:nT+2),s, 'r--');
hold on;

[tp, esp] = meshgrid(tint(2:nT+2), x(2:nx+1));
% Going back to the original variable
esp2 = esp*diag(s);
y = U(2:nx+1,:);
%contourf(tp,esp2,y);
surf(tp,esp2,y);