# How to solve 2D Heat equation on rectangle with NDSolve?

I am trying to solve the 2D heat equation $$\Delta u =u_t$$ with Dirichlet boundary conditions on the rectangle $$(0,a)\times (0,b)$$ and with initial condition $$u(x,y,0)=g(x,y)$$. This is my code

     g[x_,y_]:=x^2+y^2
a=1;
b=2;
s = NDSolve[{Laplacian[un[x, y, t], {x, y}] == D[un[x, y, t], {t, 1}],
DirichletCondition[un[x, y, t] == 0, True],
un[x, y, 0] == g[x, y]},
un, Element[{x, y},Rectangle[{0, a}, {0, b}]]] // Flatten

ParametricPlot3D[{x, y, un[x, y, 1] /. s}, {x, 0, a}, {y, 0, b}, PlotRange -> All, BoundaryStyle -> Directive[Blue, Thick]]


The code is not working, what is my mistake?

• You are using numerical solver but have undefined g[x, y]. How to expect it to be solved numerically? Mathematica is very powerful software but even it can not guess what $g(x,y)$ could be without telling it. Commented Jun 15, 2022 at 19:52
• That's not my mistake, I had defined a function $g[x,y]$ but I forgot to put it here. The code still didn't work.
– user86304
Commented Jun 16, 2022 at 8:18
• but I forgot to put it here OK. I voted to open it now. Once it is open I'll post the answer. I just solved it for you. could not do it without know what your $g(x,y)$ is. Or you can post new question with the correct input instead. Commented Jun 16, 2022 at 8:55

The code is not working, what is my mistake?

Your Rectangle was wrong Rectangle[{0, a}, {0, b}]. Should be Rectangle[{0, 0},{a, b}]]. You did not have time range specification in the call to NDSolve. Only the spatial. That is not enough for numerical solver.

ClearAll[x, y, t, a, b];
g[x_, y_] := x^2 + y^2
a = 1;
b = 2;
pde = Laplacian[u[x, y, t], {x, y}] == D[u[x, y, t], t]
bc = DirichletCondition[u[x, y, t] == 0, True]
ic = u[x, y, 0] == g[x, y]
sol=NDSolveValue[{pde, bc, ic}, u,Element[{x, y},Rectangle[{0, 0},{a, b}]] ,{t, 0, 2}]


Animate[Plot3D[sol[x, y, t],
Element[{x, y}, Rectangle[{0, 0}, {a, b}]],
PlotRange -> {Automatic, Automatic, {0, 5}} ,
PerformanceGoal -> "Quality", ColorFunction -> "SolarColors",
Mesh -> True], {t, 0, 2}, AnimationRate -> .3]


This can also be solved exactly by DSolve

ClearAll[x, y, t, a, b];
g[x_, y_] := x^2 + y^2
a = 1;
b = 2;
pde = Laplacian[u[x, y, t], {x, y}] == D[u[x, y, t], t]
bc = {u[0, y, t] == 0, u[a, y, t] == 0, u[x, 0, t] == 0, u[x, b, t] == 0}
ic = u[x, y, 0] == g[x, y]
sol = DSolveValue[{pde, bc, ic}, u[x, y, t], {x, y, t}]