# Steady State Numerical Solution of Heat Diffusion Equation in 2D

I'v been trying to solve for the steady state of the heat diffusion equation numerically but I cant seem to get it to work. My code is

pde = D[u[x, y], {x, 2}] + D[u[x, y], {y, 2}] == 0
sol = NDSolve[{pde, u[x, 0] == 100, u[x, 10] == 400, u[0, y] == 0,
u[10, y] == 0}, u[x, y], {x, 0, 10}, {y, 0, 10}]


and I get this error

NDSolve::ivone: Boundary values may only be specified for one independent
variable. Initial values may only be specified at one value of the other
independent variable.


I think It has something to do with there not being an initial value condition, but I don't have one since I'm solving for the steady state. Could anyone please tell me where the above error comes from? Any help would be appreciated. Thanks.

• what is u[10,10]? 0 or 400? – chris Jun 24 '15 at 17:47
• You can cook up a boundary condition without the corner discontinuity issue and you still get that same error message. – george2079 Jun 24 '15 at 17:52
• It works by me ! Mathematica 10.0.2 on Windows 7 – andre314 Jun 24 '15 at 17:56
• works here as well: Mathematica 10.1 on Mac OS X 10.10.4 – sunt05 Jun 24 '15 at 17:59
• On Mathematica 8.0.4 (Windows 7), there is effectively the error – andre314 Jun 24 '15 at 18:48

If you have v.10 you can explicitly use the finite element method:

Needs["NDSolveFEM"]
mesh = ToElementMesh[Rectangle[{0, 0}, {10, 10}]]
sol = First@NDSolveValue[{Laplacian[w[x, y], {x, y}] == 0,
DirichletCondition[w[x, y] == 100, y == 0],
DirichletCondition[w[x, y] == 400, y == 10],
DirichletCondition[w[x, y] == 0, x == 0],
DirichletCondition[w[x, y] == 0, x == 10]}, {w}, {x, y} \[Element] mesh];

ContourPlot[sol[x, y], {x, 0, 10}, {y, 0, 10}] • You don't have to explicitly tell Mathematica 10 to use the finite element method. With the syntax the OP uses, it switches automatically to the finite element method. – andre314 Jun 24 '15 at 19:47

Since you are specifically asking about versions below 10, it may be useful to point out that this problem is equivalent to the electrostatics problem of finding the potential in a region bounded by conductors held at fixed voltages. This can be solved, e.g., with the simple relaxation method I implemented in this answer, where I actually allow for lots of other complications. The potential is solved on a rectangular grid, and the main point was to play around with a graphical method for specifying the boundary conditions. You can apply that method to the question here as follows (I'm first copying the relevant definitions from the lined answer):

defaultGridSize = 120;
step =
Compile[{{phi, _Real, 2}, {orig, _Real, 2}, {mask, _Real,
2}, {dchiX, _Real, 2}, {dchiY, _Real, 2}},
Module[{f = (RotateRight[phi] + RotateLeft[phi] +
RotateRight[phi, {0, 1}] + RotateLeft[phi, {0, 1}])*mask +
orig, ex, ey},
ex = (RotateRight[phi, {0, 1}] - RotateLeft[phi, {0, 1}]);
ey = (RotateRight[phi] - RotateLeft[phi]);

iterate =
Compile[{{gridArray, _Real, 2}, {originalArray, _Real,
2}, {maskArray, _Real, 2}, {dchiX, _Real, 2}, {dchiY, _Real,
2}, {tol, _Real}},
FixedPoint[step[#, gridArray, maskArray, dchiX, dchiY] &,
originalArray, 100000,
SameTest -> (Max@Abs@Flatten[#1 - #2] < tol &)]];

digitize[gr_, n_] :=
N@ImageData@
ColorConvert[
Image[Show[gr, Background -> Black,
BaseStyle -> {Antialiasing -> False}], ImageSize -> n],
"GrayScale"];

createLandscape[conductors_, chargePlus_, chargeMinus_,
suceptibility_, nGrid_] :=
gridConductors = digitize[conductors, nGrid];
gridRho = (digitize[chargePlus, nGrid] -

poissonSolver[conductors_, chargePlus_, chargeMinus_, suceptibility_,
nGrid_: defaultGridSize, tolerance_: 10^(-6)] :=
Block[{averagePotential, gridConductors, gridRho, gridChi, gridEps,
initialGrid}, {gridConductors, gridRho, gridChi, maskList} =
createLandscape[conductors, chargePlus, chargeMinus,
suceptibility, nGrid];
averagePotential = Mean[Select[Flatten@gridConductors, Positive]];
gridEps = 1. + gridChi;
gridList = gridConductors + gridRho/(4.*gridEps);
dChiYList = (RotateLeft[gridChi] -
RotateRight[gridChi])/(2. gridEps);
dChiXList = (RotateLeft[gridChi, {0, 1}] -
RotateRight[gridChi, {0, 1}])/(2. gridEps);
Reverse@
tolerance]];


Now for the specific boundary conditions. The region extends from 0 to 10 in both directions, but I need to add some padding on the outside in order to define areas with the desired boundary potential. So the simulation area is the square extending from -1 to 11 in both directions (this is where the figure shows rectangles in various gray shades):

conductors =
Graphics[{{GrayLevel[.5], Rectangle[{0, 10}, {10, 11}]},
{GrayLevel[.1], Rectangle[{10, 0}, {11, 11}]}, {GrayLevel[.2],
Rectangle[{-1, -1}, {11, 0}]},
{GrayLevel[.1], Rectangle[{-1, 0}, {0, 11}]}}, To explain how I chose the conductor potentials: in my graphical approach, the color black corresponds to the absence of any object. All other potentials are encoded in the grayscale value of the region. This means the value 0 can't be used as a potential, and therefore I shift all potentials by a value 0.1 to make them non-zero. Also, I rescaled all potentials so that they fit into the interval from .1 to 1 that can be directly translated to GrayLevel. This means I have to use the following correspondence between your boundary conditions and mine:

• 0 $\mapsto$ .1
• 100 $\mapsto$ .2
• 400 $\mapsto$ .5.

There is no need to specify any of the other electrostatic quantities, except to define their corresponding graphics objects to be empty. Then I invoke the solver:

chargePlus =
Graphics[{}, PlotRange -> (PlotRange /. FullOptions[conductors])];

chargeMinus =
Graphics[{}, PlotRange -> (PlotRange /. FullOptions[conductors])];

susceptibility =
Graphics[{}, PlotRange -> (PlotRange /. FullOptions[conductors])];

Timing[potential =
poissonSolver[conductors, chargePlus, chargeMinus, susceptibility];]

(* ==> {5.65001, Null} *)

ListPlot3D[potential, PlotRange -> All,
PlotStyle -> {Orange, Specularity[White, 10]},
DataRange->{{-1,11},{-1,11}}] For more details, you may want to consult the link I gave above.