Solving a 2D heat equation on a square with Dirichlet boundary conditions

Question

I am trying to solve the following heat equation problem on the square [0,1]x[0,1]. \begin{equation*} \begin{gathered} u_t = u_{xx} + u_{yy} + f(x,y,t), \qquad u(x,y,0) = 0, \qquad u=0 \text{ on boundary} \\ f(x,y,t) = 10(\sin\pi x)^{10}(\sin\pi y)^{10}(\sin\pi t)^{10}. \end{gathered}\end{equation*}

I am trying to solve it with Mathematica using NDSolve but I am stuck and unsure of what I have done. I have

NDSolve[{Derivative[0, 0, 1][u][x, y, t] == 
   Derivative[0, 2, 0][u][x, y, t] + 
    Derivative[2, 0, 0][u][x, y, t] + 
    10*((Sin[Pi*x])^10)*((Sin[Pi*y])^10)*((Sin[Pi*t])^10), 
  u[x, y, 0] == 0}, u, {x, 0, 1}, {y, 0, 1}, {t, 0, 1}]

I am not sure that this is actually correct, and I do not know how to enforce the zero boundary condition on my square. Also, I think this equation has an exact solution. In the event that it does, can I also solve it analytically using Mathematica? I would appreciate any help on this as I am still a complete novice in Mathematica. I thank all helpers.

Answer 1

ClearAll[x, y, t, f];

f = 10*Sin[Pi*x]^10*Sin[Pi*y]^10*Sin[Pi*t]^10;
pde = D[u[x, y, t], t] == D[u[x, y, t], {x, 2}] + D[u[x, y, t], {y, 2}] + f;
ic = u[x, y, 0] == 0;
bc = {u[0, y, t] == 0, u[1, y, t] == 0, u[x, 0, t] == 0, u[x, 1, t] == 0};

sol = NDSolve[{pde, ic, bc}, u, {x, 0, 1}, {y, 0, 1}, {t, 0, 4}];

Manipulate[
 Plot3D[Evaluate[u[x, y, t0] /. sol], {x, 0, 1}, {y, 0, 1}, 
  PlotRange -> {Automatic, Automatic, {-0.1, 0.15}}, 
  PerformanceGoal -> "Quality"],
 {{t0, 0, "time"}, 0, 4, 0.01, Appearance -> "Labeled"},
 TrackedSymbols :> {t0}
 ]

NDSolve did give some warnings on console.

General::munfl: 6.17211*10^-153 7.58767*10^-160 is too small to represent 
as a normalized machine number; precision may be lost.

These look like from evaluating f at some points. May be with some additional options these can be eliminated.

For analytical solution, DSolve taking long time. So stopped it. It does not look like it can solve it analytically.

sol = DSolve[{pde, ic, bc}, u[x, y, t], {x, y, t}]

Mathematica 12.2 on windows 10

Answer 2

As an alternative, you can use HeatTransferPDEComponent and FEM to set up the model.

pde = HeatTransferPDEComponent[{u[x, y, t], t, {x, y}}, <||>] == 
   10*Sin[Pi*x]^10*Sin[Pi*y]^10*Sin[Pi*t]^10;
ic = u[x, y, 0] == 0;
bc = DirichletCondition[u[x, y, t] == 0, True];
sol = NDSolve[{pde, ic, bc}, 
   u, {x, y} \[Element] Rectangle[], {t, 0, 4}];

Also, note the Heat Transfer PDEs and Boundary Conditions guide page, the Heat Transfer tutorial, the Heat Transfer PDE models and the verification examples.