# Incorrect solution of 3D PDE with Neumann conditions?

I would like to solve the following PDE: $$\Delta \mu (x,y,z)=\mu$$ Over the parallelepiped: $x \in [0,100]$, $y \in [0,100]$ and $z \in [-10,0]$

With the following Neumann boundary conditions: $$\mathbf n \cdot \nabla \mu |_{z=0}= x y, \ \ \ \ \ \ \ \mathbf n \cdot \nabla \mu |_{ rest \ of \ boundary}= 0$$

I use Mathematica 11.0. First, I specify the region:

reg = Parallelepiped[{0, 0, -10}, {{100, 0, 0}, {0, 100, 0}, {0, 0, 10}}]


Then I solve the diffusion equation:

sol = NDSolveValue[Laplacian[mu[x, y, z], {x, y, z}] - mu[x, y, z] == NeumannValue[- x y, z == 0]},
mu, {x, y, z} \[Element] reg]


I get a solution (however, I get the mistake NDSolveValue::femibcnd: No DirichletCondition or Robin-type NeumannValue was specified for {mu}; the result may not be unique.). Afterwards, I want to check whether the second boundary condition is satisfied (in documentation it is said that if boundary condition is not specified, it is automatically assumed to be Neumann 0). Therefore, I try to calculate the flow of the function at a random point of the boundary at the lower boundary:

D[sol[1, 1, z], z] /. z -> -10


And obtain -0.00125585, and not 0 as I expected. What is wrong with my solution?

Also, I am confused with the result not being smooth. See the following:

Plot[Evaluate[sol[1, 1, z], {z, -10, 0}]]


Another check shows that probably accuracy of the solution is unsufficient even at the boundary, however, I am not sure how to improve it without significant cost to the speed of computation:

Plot[{Evaluate[sol[1, y, 0]], y}, {y, 0, 100}, PlotRange -> All]


I tried using right-hand side in the following form: NeumannValue[- x y, z == 0] + NeumannValue[0, z != 0], which doesn't affect the answer.

• Is it normal that you don't have time in the "diffusion" equation? Also, did you notice the error NDSolveValue::femibcnd: No DirichletCondition or Robin-type NeumannValue was specified for {mu}; the result may not be unique.? – anderstood Mar 5 '18 at 14:49
• You could also try with NDEigensystem, since you are trying to find the eigenfunctions of the Laplacian for the eigenvalue 1. But that does not explain the result you observe... – anderstood Mar 5 '18 at 14:56
• Did you try NeumannValue[x y, z == 0] + NeumannValue[0, z != 0] as right hand side of your pde? – Ulrich Neumann Mar 5 '18 at 15:09
• @HenrikSchumacher could you please incorporate on your answer? What is the linear system being built in the background and why it is like this? What can possibly be done? – Gretchen Mar 5 '18 at 15:15
• NDSolve employs the finite element element method: The differential equation gets discretized into a finite dimensional linear system and this gets solved afterwards. If the infinite dimensional equations are not well formulated then the finite dimensional system will be badly behaved. But I have correct myself: Since the differential operator is $\Delta -1$, the system does not have a null space at all. Long story short: I do not know why this happens. I would have expected that Ulrich's answer would settle this problem... – Henrik Schumacher Mar 5 '18 at 15:31

I think this is somewhat an artefact of the discretization. Note that Neumann boundary conditions are only applied implicitely in the weak formulation of the PDE. In the weak formulation, the Neumann condition may be indeed fulfilled, but not in the way you compute it by differentiating with (interpolated) solution of the discrete system. Note, that we can decrease the mesh size for the discretization, e.g. as follows.

reg0 = Parallelepiped[{0, 0, -10}, {{100, 0, 0}, {0, 100, 0}, {0, 0, 10}}];
reg = DiscretizeRegion[reg0, MaxCellMeasure -> 5]; sol =
NDSolveValue[
Laplacian[mu[x, y, z], {x, y, z}] - mu[x, y, z] ==
NeumannValue[-x y, z == 0], mu, {x, y, z} \[Element] reg];


When I compute the Neumann differential as you did, then its absolute is smaller than before:

D[sol[1, 1, z], z] /. z -> -10

-0.000587538


This is somewhat an indicator for the fact that the two notions of Neumann boundary conditions (weakly and strongly formulated) coincide only in the limit $\text{mesh size} \to 0$.

Up front, I can not explain exactly what is going on but I do not think this is a bug. Let's look at a simplified 2D version of this:

Needs["NDSolveFEM"]
w = 100;
l = 10;
r = Rectangle[{0, -l}, {w, 0}];
m = ToElementMesh[r,(*"MeshElementType"\[Rule]TriangleElement,*)
"MeshOrder" -> 2, MaxCellMeasure -> 10];
(*m["Wireframe"]*)

nval = y;
sol = NDSolveValue[{Inactive[
Div][-1. Inactive[Grad][mu[y, z], {y, z}], {y, z}] +
mu[y, z] == NeumannValue[nval, z == 0],
DirichletCondition[mu[y, z] == 0, z == -l]},
mu, {y, z} \[Element] m];


I have made a few changes to the PDE. I put this in inactive form (does not make a difference) and I also added a DirichletCondition at z==0 - does not matter either but rules out that what we see is an issue due to an neumann only BC. It does not seem so.

D[sol[1, z], z] /. z -> -l
0.00008033698364150563


Looks reasonable and not too far off compared to a Neumann zero value at that end.

Plot[sol[1, z], {z, -l, 0}, PlotRange -> All]


Note, that we also see this slight kink; but much less than with TriangleElement. This is also true for the 3D version. Using Cuboid and HexahedronElement makes this behave better.

Plot[Evaluate[{sol[y, z] - nval} /. {mu -> sol, z -> 0}], {y, 0, w},
PlotRange -> All]


Note the drift in the difference between the neumann value and the computed value and also the jump at both ends. The jumps at both ends, I think, are there because the neumann zero BC also needs to be satisfied. The drift can be reduced is the mesh is refined but it remains.

Also, the

ContourPlot[sol[y, z], {y, 0, w}, {z, -l, 0}, PlotRange -> All]
Plot3D[sol[y, z], {y, 0, w}, {z, -l, 0}, PlotRange -> All]


look much smoother when using QuadElement (or HexahedronElement in 3D)

I have checked with other FEM software and I get a similar enough result to believe that this is not an implementation bug.

• thank you a lot for your contribution! Answering to your question: I need to solve a similar problem in my research, and $x y$ is just a sample function for understanding how to solve such equations in Mathematica. Could you please specify what form of PDE exactly you used to obtain these plots? How did you incorporate the BC in the Dirichlet form? What is the variable nval which you use in for the last plot? I struggle to understand your answer without this information :) – Gretchen Mar 7 '18 at 15:33
• @Gretchen, sorry forgot the NDSolve call. Hope this is complete now. – user21 Mar 7 '18 at 15:44

I tried to solve the problem inside a smaller parallelepiped(based on the approach from Henrik Schumacher) because I assumed a scaling problem.

reg0 = Parallelepiped[{0,0, -1 }, {{1 , 0, 0}, {0, 1 , 0}, {0, 0, 1 }}];
reg = DiscretizeRegion[reg0 , MaxCellMeasure -> {"Volume" -> .001}] ;
sol = NDSolveValue[Laplacian[mu[x, y, z], {x, y, z}] - mu[x, y, z] ==
NeumannValue[x y, z == 0] + NeumannValue[0, z != 0],mu, {x, y, z} \[Element] reg]


The neumann conditions can be visualized

Plot3D[Evaluate[Derivative[0, 0, 1][sol][x, y, 0]], {x, 0, 1}, {y, 0,1}, PlotRange -> {-1, 1}] (* z=0: mu= x y *)


Plot3D[Evaluate[Derivative[0, 0, 1][sol][x, y, -1]], {x, 0, 1}, {y, 0,1}, PlotRange -> {-1, 1}]


and both plots agree well with the conditions given in NDSolve. Finally the examplary verification

Plot [ sol[1, 1, z] , {z, -1, 0 } , PlotRange -> {-1, 0}]
`

becomes smooth.

I'm still unsure what happens, perhaps scaling or meshing issue?