I am looking for a way to find the largest subdomain $U$ of the unit disk $D(0,1)$ for which the solution to $$ \begin{align} \Delta u&=0.3-y,\quad (x,y)\in U\\ u&=0,\quad (x,y)\in\partial U \end{align} $$ is negative in $U$.

Some additional details

We start by looking at the problem in the disk. Solving the Dirichlet problem $$ \begin{align} \Delta u&=0.3-y,\quad (x,y)\in D(0,1)\\ u&=0,\quad (x,y)\in\partial D(0,1) \end{align} $$ using

usol1[x_, y_] =
NDSolveValue[{Laplacian[u[x, y], {x, y}] == 0.3 - y, 
 DirichletCondition[u[x, y] == 0, True]}, 
 u[x, y], 
 Element[{x, y},Disk[]]

gives a solution that is negative in some part (blue) and positive in some part (pink):

ContourPlot[Evaluate@usol1[x, y], Element[{x, y},Disk[]], 
 Contours -> {0, -0.05, -0.1}, 
 ContourStyle -> {Black, Black, {Thick, Black}}, 
 ContourShading -> {Darker[Blue], Blue, Lighter[Blue], Pink}, 
 PlotPoints -> 100, AspectRatio -> 1]

first contour plot, showing a solution changing signs

If we instead try a smaller domain $$ \Omega=D(0,1)\cap \{(x,y)\in D(0,1)~|~y<0.3\} $$ defined in Mathematica as

\[CapitalOmega] := 
  RegionDifference[Disk[], Rectangle[{-2, 0.3}, {2, 2}]]

We solve $$ \begin{align} \Delta u&=0.3-y,\quad (x,y)\in \Omega\\ u&=0,\quad (x,y)\in\partial \Omega \end{align} $$ with

usol2[x_, y_] = 
NDSolveValue[{Laplacian[u[x, y], {x, y}] == 0.3 - y, 
 DirichletCondition[u[x, y] == 0, True]}, 
 u[x, y], 
 Element[{x, y}, \[CapitalOmega]]

and show the corresponding contour plot with

ContourPlot[Evaluate@usol2[x, y], 
 Element[{x, y}, \[CapitalOmega]], 
 Contours -> {0, -0.05, -0.1}, 
 ContourStyle -> {Black, Black, {Thick, Black}}, 
 ContourShading -> {Darker[Blue], Blue, Lighter[Blue], Pink}, 
 PlotPoints -> 100, AspectRatio -> 1.3/2]

second contour plot, showing a positive solution

We find a strictly negative solution.

Now, clearly the domain $\Omega$ is not the largest possible, so I restate my question,

Is there some way, using Mathematica, probably with some numerical method, to find the largest subdomain $U\subset D(0,1)$ such that the Poisson problem $$ \begin{align} \Delta u&=0.3-y,\quad (x,y)\in U\\ u&=0,\quad (x,y)\in\partial U \end{align} $$ has a strictly negative solution in $U$?


Just to clarify, I do not want to find the area of the part where my usol1 is negative.

  • $\begingroup$ You might want to have a look at mathematica.stackexchange.com/questions/69473/… $\endgroup$
    – Lukas Lang
    Commented Aug 14, 2017 at 7:10
  • $\begingroup$ @Mathe172 Thank you for your suggestion. I don't see exactly how that applies. Do you mean that I should discretize the disc and then solve the Poisson equation on parts of the mesh in a trial and error way? $\endgroup$
    – mickep
    Commented Aug 14, 2017 at 7:16
  • $\begingroup$ Why are you sure that there exists a largest domain with this property? Does the union of two of such domains also have this property? Are all domains allowed or do you have only few-parameter families of domains in mind? $\endgroup$ Commented Aug 23, 2017 at 22:15
  • $\begingroup$ Do you restrict your attention to topological disks only (such that the boundary will always have a single component)? What do you mean by largest domain? Domain with largest area? Or union of all such domains? $\endgroup$ Commented Aug 23, 2017 at 22:20
  • $\begingroup$ @HenrikSchumacher Thank you for your comments, and many relevant questions. I am not really sure that there is a largest domain. Looking at the situation for different $\Omega_\alpha=\{(x,y)\in D(0,1)~|~y<\alpha\}$ somehow suggests it, but of course, if you have several components, that could be a different story. I would be happy to find a method that somehow starts with the domain in the second example in the question, and somehow pushes the boundary (I know that is not uniquely defined) as far as possible while keeping the negative solution. A bit vague, I know. By largest I mean in area. $\endgroup$
    – mickep
    Commented Aug 24, 2017 at 6:32

1 Answer 1


Define the unit disc centred on the origin. This is used only to highlight the unit disc in the contour plot below.

disk = Disk[{0, 0}, 1];

Define a grid of points on the boundary of the disk. The parameter n controls the fineness of the grid.

n = 20;
angles = Range[0, 2Pi-Pi/n, Pi/n];
boundarypts = AngleVector /@ angles;

From here on down one iteration of the algorithm is being described. Iterations 2, 3, ... should begin at this point. The solution converges quickly, so I have left the iterations as being manual, rather than wrap the code up in a loop (or whatever).

Define the corresponding boundary mesh region. Updated versions of this mesh region are used to define the region over which the partial differential equation is to be solved.

region = BoundaryMeshRegion[boundarypts, Line[Append[Range[2 n], 1]]];

Solve the PDE using the current boundary mesh region.

soln = NDSolve[{Laplacian[u[x, y], {x, y}] == 0.3 - y, DirichletCondition[u[x, y] == 0, True]}, u, {x, y} \[Element] region];

Substitute this solution into a function for ease of use.

u0[{x_, y_}] = u[x, y] /. soln[[1]];

Generate a contour plot of the solution. The boundary mesh is highlighted in red, and the original unit disc is highlighted in blue.

plot = ContourPlot[u0[{x, y}], {x, y} \[Element] disk, Contours -> {0, -0.05, -0.1}, BoundaryStyle -> Red, Epilog -> {Blue, Circle[]}, PlotRange -> {{-1, 1}, {-1, 1}}, PlotRangePadding -> 0.1]

Dirichlet Poisson solution - contour plot

Display the boundary mesh region - i.e. the solution to the posed problem.


Dirichlet Poisson solution - boundary mesh region

The above plot shows the solution after 20 update iterations, using the update procedure described below.

Compute the sign of the solution (the "signature") just inside each boundary point. Where the signature is negative/positive the boundary should be moved outwards/inwards (respectively).

signatures = u0 /@ (0.99 boundarypts);
maxsig = Max@Abs@signatures;

Update the boundary points:

  1. If a signature is negative then move the corresponding boundary point outwards a bit, and truncate it back to the unit circle if it jumps outside the circle.
  2. If the signature is positive then move the corresponding boundary point inwards a bit.

boundarypts = MapThread[ Switch[Sign[#1], -1, If[Norm[#] > 1, Normalize[#], #] &[(1 - #1) #2], 1, #2/(1 + #1), _, #2] &, {0.3/maxsig signatures, boundarypts}];

This algorithm for moving the boundary points is very crude, and with some experimentation it could no doubt be improved.

  • 1
    $\begingroup$ Thank you very much for your answer. I got your code working, and this is exactly in the line I was looking for! I will surely be able to tune it to fit my needs. Also, I'm happy you answered before the bounty expired, so that I can award it to you. $\endgroup$
    – mickep
    Commented Aug 24, 2017 at 19:50

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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