3
$\begingroup$

I am trying to solve Laplace's equation in 2 dimensions with potencial boundary conditions at the edges of an external square and an internal circle. Everything seems to work fine until the point where I try to obtain a DensityPlot to visualize the the obtained potential at all points. When I execute the density plot line, I am returned with a diagram that displays the solution of the equation for only a fraction of the region where I solved the equation. How can I fix this ? The code is as follows:

α = 
 RegionDifference[Rectangle[{-50, -50}, {50, 50}], Disk[{0, 0}, 5]]
sol = NDSolveValue[{D[u[x, y], x, x] + D[u[x, y], y, y] == 0, 
   DirichletCondition[u[x, y] == 100, x^2 + y^2 == 25], 
   u[x, -50] == u[x, 50] == u[-50, y] == u[50, y] == 0}, 
  u, {x, y} ∈ α]

DensityPlot[sol[x, y], {x, y} ∈ α, Mesh -> None, 
 ColorFunction -> "Rainbow", PlotRange -> All, 
 PlotLegends -> Automatic]

This is what I mean by "only a fraction of the region where I solved the equation is shown.

$\endgroup$
  • $\begingroup$ I filed this as a bug. $\endgroup$ – user21 Nov 9 '18 at 8:26
3
$\begingroup$

It's a bad day for Regions, I guess...

Until NDSolve will dicretize it correctly, you can discretize the region yourself:

α = DiscretizeRegion[
   RegionDifference[Rectangle[{-50, -50}, {50, 50}], 
    Disk[{0, 0}, 5]],
   MaxCellMeasure -> {1 -> 1}
   ];
sol = NDSolveValue[{D[u[x, y], x, x] + D[u[x, y], y, y] == 0, 
    DirichletCondition[u[x, y] == 100, x^2 + y^2 == 25], 
    u[x, -50] == u[x, 50] == u[-50, y] == u[50, y] == 0}, 
   u, {x, y} ∈ α];
DensityPlot[sol[x, y], {x, y} ∈ α, Mesh -> None, 
 ColorFunction -> "Rainbow", PlotRange -> All, 
 PlotLegends -> Automatic]

enter image description here

$\endgroup$
  • 1
    $\begingroup$ Using DiscretizeRegion is not ideal for FEM based computations. See my answer for details. $\endgroup$ – user21 Nov 9 '18 at 8:23
  • $\begingroup$ Thanks a lot for your continuous edits of my lazily written posts! $\endgroup$ – user21 Nov 9 '18 at 8:29
  • 1
    $\begingroup$ Dito. That's my way of showing that I read them! ;) $\endgroup$ – Henrik Schumacher Nov 9 '18 at 8:30
4
$\begingroup$

Using DiscretizeRegion is not the best solution. This is because DiscretizeRegion will only provide a linear approximation for a region. A better approximation is to use ToElementMesh which can generate a second order accurate approximation to the region:

α = RegionDifference[Rectangle[{-50, -50}, {50, 50}], Disk[{0, 0}, 5]]
sol = NDSolveValue[{D[u[x, y], x, x] + D[u[x, y], y, y] == 0, 
   DirichletCondition[u[x, y] == 100, x^2 + y^2 == 25], 
   u[x, -50] == u[x, 50] == u[-50, y] == u[50, y] == 0}, 
  u, {x, y} ∈ α]
sol["ElementMesh"]["Wireframe"]

enter image description here

DensityPlot understands ElementMesh:

DensityPlot[sol[x, y], {x, y} ∈ sol["ElementMesh"], 
 Mesh -> None, ColorFunction -> "Rainbow", PlotRange -> All, 
 PlotLegends -> Automatic]

enter image description here

Now to look at the approximation quality difference:

N[exact = Integrate[1, {x, y} ∈ α]]
9921.460183660256`

exact - Area[DiscretizeRegion[α]]
-2.610161994172813`

exact - Total[sol["ElementMesh"]["MeshElementMeasure"], 2]
-0.004173680044914363`

Clearly the use of an ElementMesh is superior here as it much better approximates the curved boundaries.

This is explained in more detail in the section Comparing ElementMesh and MeshRegion in the Element Mesh Generation tutorial.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.