# Confusion with Neumann and Dirichlet

If I solve Laplace's equation with Neumann boundary conditions then everything is defined via derivatives. Consequently one needs to fix a point with a specific value to get a solution. However if I fix my value with a Dirichlet condition the solution is distorted. Where am I going wrong?

Edit to question I think I have been asking Mathematica to solve an impossible problem. See my answer below.

Continue with original question

Here is an example

 Needs["NDSolveFEM"];
x2 = 4; y2 = 1;

reg = ImplicitRegion[0 <= x <= x2 && 0 <= y <= y2, {x, y}];

mesh = ToElementMesh[reg,
"BoundaryMeshGenerator" -> {"Continuation"},
MaxCellMeasure -> .002,
"MaxBoundaryCellMeasure" -> 0.01];

Show[mesh["Wireframe"], Frame -> True, PlotRange -> All] Here is a Laplacian with Neumann boundary conditions so that everything is defined through derivatives.

sol = NDSolveValue[{
Laplacian[u[x, y], {x, y}] ==
NeumannValue[1, 0 <= y <= y2 && x == 0] +
NeumannValue[Cos[2 \[Pi] x/x2], 0 <= x <= x2 && y == y2]
}, u, {x, y} \[Element] mesh];


This gives a message along the expected lines

NDSolveValue::femibcnd: No DirichletCondition or Robin-type NeumannValue was specified for {u}; the result is not unique up to a constant. >>

Now we repeat with a Dirichlet condition that gives a value in the corner

 sol = NDSolveValue[{
Laplacian[u[x, y], {x, y}] ==
NeumannValue[1, 0 <= y <= y2 && x == 0] +
NeumannValue[Cos[2 \[Pi] x/x2], 0 <= x <= x2 && y == y2],
DirichletCondition[u[x, y] == 0, x == x2 && y == 0]
}, u, {x, y} \[Element] mesh];
Plot3D[sol[x, y], {x, y} \[Element] mesh, BoxRatios -> {x2, y2, 1}] The value in the corner is as expected but the whole solution is distorted to reach the value.

If I change the location of the Dirichlet point then the location of the distortion changes as perhaps might be expected.

sol = NDSolveValue[{
Laplacian[u[x, y], {x, y}] ==
NeumannValue[1, 0 <= y <= y2 && x == 0] +
NeumannValue[Cos[2 \[Pi] x/x2], 0 <= x <= x2 && y == y2],
DirichletCondition[u[x, y] == 0, x == x2 && y == y2]
}, u, {x, y} \[Element] mesh];
Plot3D[sol[x, y], {x, y} \[Element] mesh, BoxRatios -> {x2, y2, 1}] If I look at the Laplacian of the solution then the point continues to appear

 ClearAll[f2];
f2[x_, y_] := Evaluate[Laplacian[sol[x, y], {x, y}]]
Plot3D[f2[x, y], {x, y} \[Element] mesh, BoxRatios -> {x2, y2, 1},
PlotRange -> All]
Plot3D[f2[x, y], {x, y} \[Element] mesh, BoxRatios -> {x2, y2, 1}] The values in the above plots should be zero and they are small except where I have my Dirichlet point.

So what do I not understand? Is it something to do with Neumann values being associated with the normal to the boundary and this not being compatible with a defined point? Please enlighten me.

Following a comment from Young it may be that I have set up an impossible problem. As the equations are linear then the principle of superposition applies and we can look at the solutions from each Neumann condition separately. Starting with the second Neumann condition on its own, with the Dirichlet condition we have.

Needs["NDSolveFEM"];
x2 = 4; y2 = 1;
reg = ImplicitRegion[0 <= x <= x2 && 0 <= y <= y2, {x, y}];
mesh = ToElementMesh[reg,
"BoundaryMeshGenerator" -> {"Continuation"},
MaxCellMeasure -> .002,
"MaxBoundaryCellMeasure" -> 0.01];
sol = NDSolveValue[{
Laplacian[u[x, y], {x, y}] ==
NeumannValue[Cos[2 \[Pi] x/x2], 0 <= x <= x2 && y == y2],
DirichletCondition[u[x, y] == 0, x == x2 && y == 0]
}, u, {x, y} \[Element] mesh];
Plot3D[sol[x, y], {x, y} \[Element] mesh, BoxRatios -> {x2, y2, 1}] There is no distortion at {4,0}. Also we can check the solution by looking at the Laplacian

ClearAll[f2];
f2[x_, y_] := Evaluate[Laplacian[sol[x, y], {x, y}]]
Plot3D[f2[x, y], {x, y} \[Element] mesh, BoxRatios -> {x2, y2, 1},
PlotRange -> All] The Laplacian is good with small values everywhere.

Now we look at the solution with the first Neumann condition and the Dirichlet condition.

sol = NDSolveValue[{
Laplacian[u[x, y], {x, y}] ==
NeumannValue[1, 0 <=  y <= y2 && x == 0],
DirichletCondition[u[x, y] == 0, x == x2 && y == 0]
}, u, {x, y} \[Element] mesh];
Plot3D[sol[x, y], {x, y} \[Element] mesh, BoxRatios -> {x2, y2, 1}] The distortion is back. However do the boundary conditions prescribe a possible configuration for Laplace's equation? We have a gradient on one boundary and zero gradient on the other boundaries. From a physical fluid mechanics viewpoint this is like having an inlet on one edge but no outlet. Such a condition is not physically possible for an incompressible fluid. We may illustrate the solution using a StreamPlot.

ClearAll[f];
f[x_, y_] := Evaluate[Grad[sol[x, y], {x, y}, "Cartesian"]];
StreamPlot[f[x, y], {x, 0, x2}, {y, 0, y2}, AspectRatio -> Automatic] There is an inflow on the left and an outflow through the Dirichlet point. This looks like a physical solution but the distortion around the Dirichlet point has become essential to let the fluid out.

In conclusion it looks like I was trying to solve an impossible problem and Mathematica gave me a possible answer to a slightly different problem. I guess it would be difficult to flag up a warning message that the impossible was being asked for.

• I think this may be more of a question for scicomp. I am not sure this is a Mathematica issue per se but more of a FEM issue in general. If you find out what makes this NeumannCondition so special I'd appreciate if you could let me know. Jul 16 '16 at 16:53

Update

As a follow-up to my comment, my hypothesis is that the first Neumann boundary is problematic and a variation of the problem could be solved.

Another issue seems to be associated with the meshed region. I lowered the "MeshOrder" to 1 and the Laplacian of the solution becomes perfectly 0.

Needs["NDSolveFEM"];
x2 = 4; y2 = 1;
reg = ImplicitRegion[0 <= x <= x2 && 0 <= y <= y2, {x, y}];
mesh = ToElementMesh[reg, "BoundaryMeshGenerator" -> {"Continuation"},
MaxCellMeasure -> .002, "MaxBoundaryCellMeasure" -> 0.01, "MeshOrder" -> 1];
Show[mesh["Wireframe"], Frame -> True, PlotRange -> All]

sol = NDSolveValue[{Laplacian[u[x, y], {x, y}] ==
NeumannValue[Cos[2 π x/x2], 0 <= x <= x2 && y == y2],
DirichletCondition[u[x, y] == 0, x == x2 && y == 0]},
u, {x, y} ∈ mesh];

Plot3D[sol[x, y], {x, y} ∈ mesh, BoxRatios -> {x2, y2, 1}, PlotRange -> All,
ColorFunction -> "Rainbow", Mesh -> None] ClearAll[f2];
f2[x_, y_] := Evaluate[Laplacian[sol[x, y], {x, y}]]
Plot3D[f2[x, y], {x, y} \[Element] mesh, BoxRatios -> {x2, y2, 1}, PlotRange -> All,
ColorFunction -> "Rainbow", Mesh -> None] • According to Help if no boundary condition is given then automatically a Neumann condition of zero gradient is used. This is what I want. Your Dirichlet conditions gives a value to the boundary and not a gradient. This avoids my problem but is solving a different problem. But thanks for your thoughts.
– Hugh
Jul 14 '16 at 18:11
• @Hugh I think there's connection to the Mesh Order. The default mesh order for ToElementMesh is 2 and a lower order seems appropriate. Thoughts? Jul 14 '16 at 19:50
• @Hugh Updated again, I placed the DirichletCondition on the wrong edge. Now x==x2 Jul 14 '16 at 20:15
• Unfortunately you have put the Dirichlet condition along the whole edge not at a point. If you put it at a point the distortion reappears. Thanks for trying.
– Hugh
Jul 14 '16 at 20:28
• As this is a linear problem solutions superimpose. The second Neumann condition on its own with the Dirichlet condition is fine. The first Neumann condition on its own with the Dirichlet condition is not good. You may be onto something here.
– Hugh
Jul 15 '16 at 8:15

If you use Dirichlet boundary condition(BC) you are imposing that your solution must be a known value at that particular point. In your case you don't know the value to specify. Your solution will have a discontinuity unless you specify a correct value. You can minimize the problem refining the mesh at the borders. Take a look at the code below:

   Needs["NDSolveFEM"];
x2 = 4; y2 = 1;
(*reg=BoundaryDiscretizeRegion[Rectangle[{0,0},{x2,y2}]];*)
reg = ImplicitRegion[0 <= x <= x2 && 0 <= y <= y2, {x, y}];
mesh = ToElementMesh[reg, "BoundaryMeshGenerator" -> {"Continuation"},
MaxCellMeasure -> 0.5, "MaxBoundaryCellMeasure" -> 0.001,
"MeshOrder" -> 1];
Show[mesh["Wireframe"], Frame -> True, PlotRange -> All]
op = Laplacian[u[x, y], {x, y}]
(*Subscript[\[CapitalGamma], D]=DirichletCondition[u[x,y]\[Equal]0, x\
\[Equal]0&&0\[LessEqual]y\[LessEqual]y2]*)
Subscript[\[CapitalGamma], D] =
DirichletCondition[u[x, y] == 0, x == x2  && y == 0]
(*Subscript[\[CapitalGamma], N]=NeumannValue[-1,y\[Equal]1 &&0\
\[LessEqual]x\[LessEqual]4]*)
Subscript[\[CapitalGamma], N] =
NeumannValue[1, x == 0 && 0 <= y <= y2] +
NeumannValue[Cos[2 \[Pi] x/x2], 0 <= x <= x2 && y == y2]
uif = NDSolveValue[{op == Subscript[\[CapitalGamma], N],
Subscript[\[CapitalGamma], D]}, u, {x, y} \[Element] mesh]
ContourPlot[uif[x, y], {x, y} \[Element] mesh,
ColorFunction -> "Temperature", AspectRatio -> Automatic]
Plot3D[uif[x, y], {x, y} \[Element] mesh, BoxRatios -> {x2, y2, 1}]
ClearAll[f2];
f2[x_, y_] := Evaluate[Laplacian[uif[x, y], {x, y}]]
Plot3D[f2[x, y], {x, y} \[Element] mesh, BoxRatios -> {x2, y2, 1},
PlotRange -> All]
Plot3D[f2[x, y], {x, y} \[Element] mesh, BoxRatios -> {x2, y2, 1}]

• I don't see why my solution has to have a discontinuity. If I add a constant to a potential function it makes no difference to the Laplacian or the Neumann values both of which are derivatives. Hence I conclude that I can add any value to the solution and it should still be a solution. Examination of your solution reveals a distortion although the Laplacian is very good.
– Hugh
Jul 14 '16 at 20:38