4
$\begingroup$

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}]
$\endgroup$
3
$\begingroup$

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]

enter image description here

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

enter image description here


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:

enter image description here

enter image description here

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.

| improve this answer | |
$\endgroup$
  • $\begingroup$ (+1) nice, well done! $\endgroup$ – user21 Jun 3 at 5:26
  • $\begingroup$ Thank you very much! For completeness, could you also include the output for the problem proposed originally? $\endgroup$ – Zac Jun 3 at 10:20
  • $\begingroup$ @zac Check my update. $\endgroup$ – xzczd Jun 3 at 10:48

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.