# Laplace equation with robin boundary conditions

I want to solve the following steady state heat transfer problem with robin boundary condition at the bottom: The following is the code for the transient solution, but how should I change the code for the steady state solution?

Clear["Global*"]
<< NDSolveFEM
L = 100;
a = 0.1;
b = 0.6;
region = Rectangle[{-L/2, -L/2}, {L/2, L/2}];
bmesh = ToElementMesh[region];

sol = NDSolveValue[{D[u[x, y, t], x, x] + D[u[x, y, t], y, y] - D[u[x, y, t], t] ==
NeumannValue[a*u[x, y, t], y == -L/2] + NeumannValue[0., x == -L/2],
u[x, L/2, t] == b,u[L/2, y, t] == 0,u[x, y, 0] == 0.
}, u, {x, y} \[Element] bmesh, {t, 0, 10},
Method -> {"MethodOfLines", "SpatialDiscretization" -> {"FiniteElement"}}, AccuracyGoal -> \[Infinity]]

DensityPlot[sol[x, y, 10], {x, y} \[Element] bmesh, PlotPoints -> 10,
MaxRecursion -> 0, Mesh -> All, ColorFunction -> "Rainbow",
PlotRange -> All, PlotLegends -> Automatic]


## Using DSolve

V 12.1 can solve this exactly.

ClearAll[u, x, y];
pde = Laplacian[u[x, y], {x, y}] == 0;
L0 = 100;
a = 1/10;
b = 6/10;
leftSide = Derivative[1, 0][u][0, y] == 0;
rightSide = u[L0, y] == 0;
bottomSide = Derivative[0, 1][u][x, 0] == a*u[x, 0];
topSide = u[x, L0] == b;
bc = {leftSide, rightSide, bottomSide, topSide};
sol = DSolve[{pde, bc}, u[x, y], {x, y}];
sol0 = sol /. K -> n (*more terms makes it more accurate*)
sol1 = u[x, y] /. First@(sol0 /. Infinity -> 200);

Plot3D[Activate[sol1], {x, 0, L0}, {y, 0, L0}] ## Using FEM

Clear["Global*"]
<< NDSolveFEM
L = 100;
a = 0.1;
b = 0.6;
region = Rectangle[{0, 0}, {L, L}]
bmesh = ToElementMesh[region]
bmesh["Wireframe"] But triangle elements can be more accurate

bmesh = ToElementMesh[Rectangle[{0, 0}, {L, L}], "MeshElementType" -> TriangleElement]
bmesh["Wireframe"] pde = Laplacian[u[x, y], {x, y}];
sol = NDSolveValue[{pde == NeumannValue[a*u[x, y], y == 0] ,
u[x, L] == b, u[L, y] == 0}, u, {x, y} \[Element] bmesh,
Method -> "FiniteElement"];
Plot3D[sol[x, y], {x, 0, L}, {y, 0, L}] DensityPlot[sol[x, y], {x, y} \[Element] bmesh, PlotPoints -> 10,
MaxRecursion -> 0, Mesh -> All, ColorFunction -> "Rainbow",
PlotRange -> All, PlotLegends -> Automatic] • A great answer. – Please Correct GrammarMistakes Mar 28 at 23:45
• Thanks, your answer helps a lot. – Jiangming Mar 30 at 15:28
• I just noticed that the FEM solution at the bottom is smaller than zero which is not reasonable. I tried to refine the mesh, but it does not help. Any ideas? – Jiangming Mar 31 at 4:38
• @Jiangming I moved NeumannValue to the right side of == in the PDE specification. Now it matches OK the analytical. Also changed the mesh to triangles, I think it is more accurate, at least that what I would expect. – Nasser Mar 31 at 6:42
• @Nasser, you could compare the solutions for the quad and triangle mesh with the analytical solution. That would tell you which one is more accurate the this problem, given the same number of elements or better the same number of degrees of freedom (Length[mesh["Coordinates"]]) – user21 Mar 31 at 7:00