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:


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));

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,:);

Your Answer

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

Browse other questions tagged or ask your own question.