10
$\begingroup$

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["NDSolve`FEM`"]; 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]

Mathematica graphics

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}]

Mathematica graphics

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}]

Mathematica graphics

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}]

Mathematica graphics

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.

Improve this question
$\endgroup$

3 Answers 3

Reset to default
3
$\begingroup$

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["NDSolve`FEM`"]; 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}]

Mathematica graphics 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]

Mathematica graphics 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}]

Mathematica graphics

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]

Mathematica graphics 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.

Improve this answer
$\endgroup$
1
  • $\begingroup$ 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. $\endgroup$ Commented Jul 16, 2016 at 16:53
3
$\begingroup$

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["NDSolve`FEM`"]; 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]

enter image description here

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]

enter image description here

Improve this answer
$\endgroup$
6
  • 2
    $\begingroup$ 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. $\endgroup$ Commented Jul 14, 2016 at 18:11
  • $\begingroup$ @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? $\endgroup$ Commented Jul 14, 2016 at 19:50
  • $\begingroup$ @Hugh Updated again, I placed the DirichletCondition on the wrong edge. Now x==x2 $\endgroup$ Commented Jul 14, 2016 at 20:15
  • $\begingroup$ 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. $\endgroup$ Commented Jul 14, 2016 at 20:28
  • 2
    $\begingroup$ 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. $\endgroup$ Commented Jul 15, 2016 at 8:15
0
$\begingroup$

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["NDSolve`FEM`"]; 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}]
Improve this answer
$\endgroup$
1
  • $\begingroup$ 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. $\endgroup$ Commented Jul 14, 2016 at 20:38

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.