# Solving Laplace's equation in 2D using region primitives

Fixed in 12.1

I was looking at examples of solving problems in electromagnetism and ran across this question, which was a great example of just that. However, looking at the Dirichlet condition, both answers define the boundary of the rectangle explicitly:

DirichletCondition[u[x, y] == 100.,
x == 40 && 40 <= y <= 60 ||
x == 60 && 40 <= y <= 60 ||
40 <= x <= 60 && y == 40 ||
40 <= x <= 60 && y == 60]


I thought to use Mathematica region capabilities to achieve the same result, and did this:

RegionBoundary @ Rectangle[{40, 40}, {60, 60}]


Much simpler and should be equivalent. But when I tried to solve the equation (same setup as in the linked question) it resulted in an error:

sol =
NDSolveValue[
{Laplacian[V[x, y], {x, y}] == 0,
DirichletCondition[
V[x, y] == 100, {x, y} ∈ RegionBoundary @ Rectangle[{40, 40}, {60, 60}]],
V[x, 0] == V[x, 100] == V[0, y] == V[100, y] == 0},
V,
{x, y} ∈
RegionDifference[Rectangle[{0, 0}, {100, 100}],
Rectangle[{40, 40}, {60, 60}]]]


No places were found on the boundary where {x,y}∈Line[{{40,40},{60,40},{60,60},{40,60},{40,40}}] was True, so DirichletCondition[V==100,{x,y}∈Line[{{40,40},{60,40},{60,60},{40,60},{40,40}}]] will effectively be ignored.

This doesn't make any sense, so I dug a little deeper. Calling RegionMember on my constructed region (FullSimplify to clean it up) gives

(x | y) ∈ Reals &&
((40 <= y <= 60 && (x == 40 || x == 60)) || (40 <= x <= 60 && (y == 40 || y == 60)))


But, if I remove the first condition

(x | y) ∈ Reals


and just copy-paste the rest in, it works (as expected, since it's the same as the explicit one given in the beginning). Just to be clear, this works:

sol =
NDSolveValue[
{Laplacian[V[x, y], {x, y}] == 0,
DirichletCondition[
V[x, y] == 100,
((40 <= y <= 60 && (x == 40 || x == 60)) || (40 <= x <=60 && (y == 40 || y == 60)))],
V[x, 0] == V[x, 100] == V[0, y] == V[100, y] == 0},
V,
{x, y} ∈
RegionDifference[
Rectangle[{0, 0}, {100, 100}],
Rectangle[{40, 40}, {60, 60}]]]


Thus the culprit seems to be this little bit:

(x | y) ∈ Reals


My question is: why? Both $$x$$ and $$y$$ are real numbers, so to me that little bit should only be extraneous, saying what's implied by the other conditions. Clearly this is not the case.

• @xzczd, thanks for checking / updating this. – user21 Mar 20 '20 at 10:26
• @user21 Oh I just added the tag, it's OP that checks it :) . – xzczd Mar 20 '20 at 10:58

Use a decimal point for your coordinate values so that regions are converted to real numbers.

ri = Rectangle[{40., 40.}, {60., 60.}];
rbi = Element[{x, y}, RegionBoundary@ri];
ro = Rectangle[{0., 0.}, {100., 100.}];
rdiff = RegionDifference[ro, ri];
rbo = Element[{x, y}, RegionBoundary@ro];
eqn = Laplacian[V[x, y], {x, y}] == 0;
dci = DirichletCondition[V[x, y] == 100, rbi];
dco = DirichletCondition[V[x, y] == 0, rbo];
eqns = {eqn, dci, dco};
sol = NDSolveValue[eqns, V, Element[{x, y}, rdiff]]
Plot3D[sol[x, y], Element[{x, y}, rdiff]] • This works, thanks. Is there some underlying reason why using integers doesn't work or is that just a bug(-ish) of Mathematica? – Sami Jan 13 '20 at 17:00
• You're welcome. I would suspect that the reason is that NDSolve is calling ToElementMesh under the hood and that is compiled C or Fortran code and that an "integer" is converted to a double precision number. That double precision number may not be the same as a Mathematica Real number. – Tim Laska Jan 13 '20 at 18:59
• @imas145, this issue seems fixed in the upcoming release. Sorry for the trouble. I think Tim's suggestion is a good workaround. – user21 Jan 14 '20 at 7:15