Solution of Dimensionless Partial Differential Equation Over a Semi-Circular Domain

Question

My problem is to solve the following BVP.

and I want to recreate a contour plot like the following graph using Mathematica.

I am trying with the following code.

ClearAll["Global`*"];
k = 5;
α = 3.0*π/4.0;
η = 0;
pde = Laplacian[u[r, θ], {r, θ}, "Polar"] == 
   k^2*u[r, θ];
bc = u[1, θ] == 
   1 - (1 - η)*HeavisideTheta[θ - α];
sol = NDSolveValue[{pde, bc}, u, {r, 0, 1}, {θ, 0, Pi}]
ContourPlot[sol[Sqrt[x^2 + y^2], ArcTan[x, y]], {x, -1, 1}, {y, 0, 1},
  Contours -> 50, ColorFunction -> "Pastel", AspectRatio -> Automatic,
  PlotLegends -> Automatic]

But the code is continuously running without showing any result.

Answer 1

Another approach:

When dealing with numerical solvers, discontinuities from functions like HeavisideTheta can introduce major hurdles. Such discontinuities might result in convergence issues or even failures in securing a solution. A typical method to mitigate this is to "smooth out" the discontinuity, allowing the numerical solver to work more effectively.

The logistic function (or sigmoid function) is one such tool to achieve this smoothing:

$$f(\theta)=\frac{1}{1+e^{-a (\theta -\alpha )}}$$

Here:

• $\alpha$ is the point of discontinuity.
• $a$ controls the "sharpness" of the transition.

For a large value of $a$, the logistic function closely approximates a step function, making it a suitable choice for this purpose. The advantage of the logistic function is that it transitions smoothly from $0$ to $1$, thereby avoiding the abrupt change that can confound numerical methods.

By using the logistic function, we create a gentle transition around the point of discontinuity, allowing the numerical solver to handle the problem more robustly and often leading to better solutions. The sharpness of this transition can be adjusted by varying $a$, but care should be taken: making the transition too sharp (by choosing a very high value for $a$) might reintroduce numerical issues.

Based on this approach, let's apply the logistic function to smooth out the discontinuity in your differential equation. Here's the solution implemented using a logistic function instead of a piecewise function:

k = 5;
α = 3*π/4;
η = 0;
a = 50;  (*Controls the sharpness of the transition*)

(*Logistic function to smooth the discontinuity*)
smoothStep[θ_, α_, a_] := 1/(1 + Exp[-a (θ - α)]);

pde = Laplacian[u[r, θ], {r, θ}, "Polar"] == k^2*u[r, θ];

bc = u[1, θ] == 1 - (1 - η)*smoothStep[θ, α, a];

sol = NDSolveValue[{pde, bc}, u, {r, 0, 1}, {θ, 0, Pi}];

ContourPlot[sol[Sqrt[x^2 + y^2], Mod[ArcTan[x, y], Pi]], {x, -1, 1}, {y, 0, 1}, 
Contours -> 20, ColorFunction -> "Pastel", AspectRatio -> Automatic, 
PlotLegends -> Automatic, PerformanceGoal -> "Quality", PlotPoints -> 75]

I hope you find it useful!

Answer 2

Replacing HeavisideTheta with Piecewise makes it work. it could be a limitation of NDSolve. I do not know.

ClearAll["Global`*"];
k=5;
α=3*π/4
η=0;
pde=Laplacian[u[r,θ],{r,θ},"Polar"]==k^2*u[r,θ];
bc=u[1,θ]==1-(1-η)*Piecewise[{{0,θ<α},{1,True}}]
sol=NDSolveValue[{pde,bc},u,{r,0,1},{θ,0,Pi}]

ContourPlot[sol[Sqrt[x^2 + y^2], ArcTan[x, y]], {x, -1, 1}, {y, 0, 1},
  Contours -> 20, ColorFunction -> "Pastel", AspectRatio -> Automatic,
  PlotLegends -> Automatic]

Note that HeavisideTheta can be written using Piecewise

Plot[1 - (1 - η)*   HeavisideTheta[θ - α], {θ, -4, 4}]

Plot[1 - (1 - η)*   Piecewise[{{0, θ < α}, {1, True}}], {θ, -4,   4}]

How to remove the black dots that appears on the boundary of the semicircle?

Just add more PlotPoints:

ContourPlot[sol[Sqrt[x^2 + y^2], ArcTan[x, y]], {x, -1, 1}, {y, 0, 1},
  Contours -> 20, ColorFunction -> "Pastel", AspectRatio -> Automatic,
  PlotLegends -> Automatic, PerformanceGoal -> "Quality", 
 PlotPoints -> 75]

This below shows the problem more clearly

ClearAll["Global`*"];
k = 5;
α = 3*π/4;
η = 0;
pde = Laplacian[u[r, θ], {r, θ}, "Polar"] == 0;
bc = u[1, θ] == 1 - (1 - η)*HeavisideTheta[θ - α];
sol = DSolveValue[{pde, bc}, u[r, θ], {r, θ}]

Now we do the same thing, but using NDSolve instead of DSovle and get an error

ClearAll["Global`*"];
k = 5;
α = 3*π/4;
η = 0;
pde = Laplacian[u[r, θ], {r, θ}, "Polar"] == 0;
bc = u[1, θ] == 1 - (1 - η)*HeavisideTheta[θ - α];
sol = NDSolveValue[{pde, bc}, u[r, θ], {r, 0, 1}, {θ, 0, Pi}]

Changing to Piecewise fixes the above.

Maybe this is a bug or could be a feature.