Solve free boundary problem for heat equation

Question

How can I use Mathematica to compute/approximate and plot the solution of the following problem?

$\min\{u_t - u_{xx} -1, u \} = 0 \text{ in } (0,T)\times (-1,1)$

$u(0,\cdot) = 0 \text{ in } (-1,1)$

$u(\cdot, -1) = u(\cdot, 1) = 0 \text{ in } (0,T) $

Of course, for the standard heat equation (without free boundary), I could just write

    heq = D[u[x, t], t] == D[u[x, t], {x, 2}] + 1;
    ic = u[x, 0] == 0;
    bc1 = u[1,t] == 0;
    bc2 = u[-1,t] == 0;
    DSolve[{heq, ic, bc1, bc2},u[x,t], {x,t}]

Answer 1

Similar questions have been asked before but none of them are unanswered, partly because the goals in these posts are set too high:

PDE with `Min` in the definition

Convergence of approximate solutions to obstacle problem for the heat equation

Luckily OP's question is answerable. However, it's not hard to notice the solution for the specific problem given by OP is just the same as that of the standard heat equation, and the Dirichlet b.c. isn't a good approximation for b.c. at infinity either. To make the answer more interesting I'd like to modify the problem a little to:

\begin{align} &\min\{u_t - u_{xx} \color{red}{+}1, u \} = 0 \text{ in } (0,T)\times (-1,1)\\ &u\Big|_{t=0} = \color{red}{e^{-50 x^2}} \text{ in } (-1,1)\\ &\color{red}{\frac{\partial u}{\partial x}\Bigg|_{x=-1} = \frac{\partial u}{\partial x}\Bigg|_{x=1}} = 0 \text{ in } (0,T)\\ \end{align}

Then let's solve. The problem cannot be solved with NDSolve, so let's turn to FDM. I'll use pdetoae for the generation of difference equation:

lhsclassic = D[u[x, t], t] - (D[u[x, t], {x, 2}] - 1);
heq = 0 == (Min[lhsclassic, u[x, t]] // PiecewiseExpand // Simplify`PWToUnitStep);

domain@x = {xL, xR} = {-1, 1}; domain@t = {t0, tend} = {0, 1/5};
icfunc[x_] = Exp[-50 x^2];

ic = u[x, t0] == icfunc[x];
bc1 = Derivative[1, 0][u][xR, t] == 0;
bc2 = Derivative[1, 0][u][xL, t] == 0;

points@x = 50; points@t = 50;
(grid@# = Array[# &, points@#, domain@#]) & /@ {x, t};
difforder = 2;

(* Definition of pdetoae isn't included in this post,
   please find it in the link above. *)
ptoafunc = pdetoae[u[x, t], grid /@ {x, t}, difforder];

delL = Rest; delboth = #[[2 ;; -2]] &;

ae = delL /@ delboth@ptoafunc@heq;
aeic = ptoafunc@ic;
aebc = delL /@ ptoafunc@{bc1, bc2};

initial[x_, t_] = icfunc[x];
sollst = Partition[#, points@t] &@
   FindRoot[{ae, aeic, aebc} // Flatten, 
     Flatten[#, 1] &@Table[{u[x, t], initial[x, t]}, {x, grid@x}, {t, grid@t}], 
     MaxIterations -> 500][[All, -1]];

solfunc = ListInterpolation[sollst, grid /@ {x, t}];

Plot3D[solfunc[x, t], {x, xL, xR}, {t, t0, tend}, PlotRange -> All]

Let's compare it to the solution of standard heat equation:

solcompare = NDSolveValue[{lhsclassic == 0, ic, bc1, bc2}, u, {x, xL, xR}, {t, t0, tend}];

Manipulate[Plot[{solfunc[x, t], solcompare[x, t]}, {x, xL, xR}, PlotRange -> {-1/10, 1}, 
  PlotStyle -> {{Thick, Blue}, {Dashed, Red, Thick}}], {t, t0, tend}]

---

The solution above is fully applicable for the original problem, of course. We just need to modify the equation, i.c. and b.c. to the original:

lhsclassic = D[u[x, t], t] - (D[u[x, t], {x, 2}] + 1);

icfunc[x_] = 0;

bc1 = u[xR, t] == 0;
bc2 = u[xL, t] == 0;

The corresponding output is:

As we can see, the solution is just the same as the standard one, because $u$ is always non-negative in the domain of definition.