I am trying to find the solution to Laplace's equation:
$\dfrac{\partial^2 u}{\partial x^2}+\dfrac{\partial^2 u}{\partial y^2}=0$
over the following region of space:
the vertices of the triangle are (-1,1), (1,1) and (0,.2). I defined the region using RegionDifferrence
.
My boundary conditions are:
$u=0$ when $-1\leq x\leq 1, y=0$ and $u=0.1$ along the two edges of the triangle enclosing the domain.
I am interested in the electric field $\nabla u$. I setup the system in NDSolve using the following code:
uval = NDSolveValue[{D[u[x, y], x, x] + D[u[x, y], y, y] == 0,
DirichletCondition[u[x, y] == 0, -1 <= x <= 1 && y == 0],
DirichletCondition[u[x, y] == .1, 0 <= x <= 1 && y == .8 x + .2],
DirichletCondition[u[x, y] == .1, -1 <= x <= 0 && y == -.8 x + .2]},
u, {x, y} ∈ Ω,
InterpolationOrder -> All,
WorkingPrecision -> MachinePrecision,
AccuracyGoal -> $MachinePrecision,
PrecisionGoal -> $MachinePrecision];
I then use DensityPlot
to visualize the derivative $u'(x,y)$:
Clearly something doesn't look right in the derivative plot as the field between the tip and the plate seems non-monotonic along the line $x=0$. I would expect the field to decrease monotonically between 0.1 units at the tip to 0 units at the plate.
Question 1: Is this a problem with Plot
or with the solution given by NDSolve
itself?
I tried playing around with PlotDesntiy
and Maxrecursion
to no avail. We see the trouble even more clearly if we plot $u_y(0,y)$ i.e. the y component of the field along the line $x=0$:
ClearAll[dy]
dy[a_, b_] := D[uval[x, y], y]^2 /. {x -> a, y -> b};
Plot[dy[0, c], {c, 0, .2}, PlotStyle -> Red, PlotRange -> All]
Clearly the field in this physical problem cannot be discontinuous as we see above. I did some searching here and came across this post, in which Question 5 seems identical to my problem: Calculating a potential function using the finite element method
The poster above is using an FEM mesh which seemed to be the source of their jagged plot. Based on this, my hypothesis is that the interpolation/stepsize used by NDSolve here is not fine enough to give me a smooth derivative. So for my problem I tried changing a variety of NDSolve
parameters such as InterpolationOrder
, Method
, MaxStepFraction
and many others, but nothing seemed to change the output.
Question 2: Should I solve my problem with a 2D FEM scheme and adopt the solution in the link above?
Finally, when I use the solution determined by NDSolve
to test it with the boundary condition for example uval[1,1]
I get the following error message:
InterpolatingFunction::dmval: "Input value {1.,1.} lies outside the range of data in the interpolating function. Extrapolation will be used."
This seems to suggest that the boundary of my defined region is not included in the domain of {x,y}? How do I include the boundary too? The way I define my region is as below, although it does output the correct value of 0.1:
Ω =
RegionDifference[
Rectangle[{-1, 0}, {1, 1}],
Triangle[{{-1, 1}, {1, 1}, {0, .2}}]];