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] + 
  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.


2 Answers 2


enter image description here

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

 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


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

pde = HeatTransferPDEComponent[{u[x, y, t], t, {x, y}}, <||>] == 
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.

  • $\begingroup$ @kroner, you are welcome. There is lots of stuff to explore. $\endgroup$
    – user21
    Feb 16, 2021 at 5:02

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.