# Trouble understanding NeumannValue and Inactive/Formal PDEs

I read all documentation about the Finite Element Method in Mathematica 10, and I read some questions here, but I'm still unable to properly understand how to use the NeumanValue expression to assign a Neumann boundary condition.

Starting with the documentation tutorial "Solving Partial Differential Equations with Finite Elements", in the section "Partial Differential Equations and Boundary Conditions" I read:

The value $g - q u$ prescribes a flux over the outward normal on some part of the boundary: $\vec{n}\cdot(c\nabla u + \alpha u + \gamma) = g - q u$

Apart some missing vector sign, I think this shoud be a $-\gamma$. With this fix I agree with the following discussion and with the consequent weak form of the PDE.

From this discussion I understand the first Part of the NeumannValue expression is related to the prescribed outward flux of the quantity inside the Div operator.

This leads to the so-called "Formal Partial Differential Equations" section where is pointed that at some time the use of Inactive is mandatory to prevent evaluation of differential operators so that NDSolve can properly interpret boundary conditions based on NeumannValue expression.

It is not so evident (for me) how things are interpreted when Div and NeumannValue are on the same side of the equation and/or prefixed with a minus sign.

To be more precise, my expectation is that, when NeumannValue expression is present on a PDE with also a Div operator, after rearranging the equation so that Div and NeumannValue are on opposite side of the equation and no minus sign is prepended to eiter expression, the first Part of the NeumannValue expression assign the outward flux of the quantity inside the Div. My expectation is that an operator like Laplacian should be probably treated as "Div @* Grad".

I'm probably wrong, because I'm unable to translate simple 1-D Poisson equations like these:

NDSolveValue[{-u''[x] == 1, u[0] == 0, u'[1] == 1}, u, {x, 0, 1}]
Plot[%[x], {x, 0, 1}]


and

NDSolveValue[{-u''[x] == 1, u[0] == 0, u'[1] == -1}, u, {x, 0, 1}]
Plot[%[x], {x, 0, 1}]


to use "Formal" notation and inactive Div, Grad, Laplacian, and to understand the following results. For example:

NDSolveValue[{Inactivate[-Div[Grad[u[x], {x}], {x}], Div | Grad] ==
1 + NeumannValue[1, x == 1], u[0] == 0}, u, {x, 0, 1},
Method -> "FiniteElement"]
Plot[%[x], {x, 0, 1}]


Ok, like the first plot: but why? NeumannValue first Part assign the flux of the first argument of Div whatever sign Div has? Maybe...

Moving the minus sign inside the Div gives me a result I cannot understand:

NDSolveValue[{Inactivate[Div[-Grad[u[x], {x}], {x}], Div | Grad] ==
1 + NeumannValue[1, x == 1], u[0] == 0}, u, {x, 0, 1},
Method -> "FiniteElement"]
Plot[%[x], {x, 0, 1}]


In my expectation, the differential equation is unchanged, it's still $-u''=1$. Assuming I missed the meaning of NeumannValue the Neumann boundary condition can be changed. But the solution plotted appear a solution for $u''=1$...

Moving minus sign inside the Grad gives errors (the equation is still $-u_{xx}=1$):

NDSolveValue[{Inactivate[Div[Grad[-u[x], {x}], {x}], Div | Grad] == 1,
u[0] == 0}, u, {x, 0, 1}, Method -> "FiniteElement"]
Plot[%[x], {x, 0, 1}]


Any help undestanding this matter is appreciated.

• I'm not sure how Inactive interacts with NDSolve, but I would have first tried NDSolveValue[{-Inactivate[Div[Grad[u[x], {x}], {x}], Div | Grad] == <..>], which seems to work. (Maybe it's a bug?) – Michael E2 Dec 21 '14 at 22:52
• @MichaelE2 I'm not sure too. But I think an Inactive equation should be equivalent to the equation obtained Activate-ing it (except Neumann boundary conditions...). – unlikely Dec 22 '14 at 8:43
• @unlikely, I fixed normal flux the sign issue in the documentation. I'll need a little more time to look at the other issue you mention. – user21 Dec 29 '14 at 16:24
• @MichaelE2, basically NDSolve sees the things that are inactive (like Div, Grad) there are a few lines that then parse out the coefficients which are then put into InitializePDECoefficients. There simply was no rule for Times[fact_, Grad[...]] which will work in a future version, if fact is NumberQ. For all other cases one should use Dot. Hope this helps. – user21 Dec 31 '14 at 13:03

I think you understood the NeumannValue documentation/concept quite well. The fact the -Grad does not give the expected result is due to a parsing error in NDSolve (which is already fixed in the development version).

NDSolveValue[{Inactivate[Div[-Grad[u[x], {x}], {x}], Div | Grad] ==
1 + NeumannValue[1, x == 1], u[0] == 0}, u, {x, 0, 1},
Method -> "FiniteElement"]


You can use:

NDSolveValue[{Inactivate[Div[{{-1}}.Grad[u[x], {x}], {x}],

• @unlikely, I had another look at this. I do not think that I will support the - inside the Inactive[Grad][- u[x]]. The mathematical proper way to specify the diffusion coefficient is via Inactivate[Div[{{-1}}.Grad[u[x],{x}],{x}],Div|Grad] The fact that -1*Inactive[Grad] is supported is a convenience but I feel that Inactive[Grad][-1 u[x],{x}] would be too much. You can try to convince me other wise, if you want. – user21 Feb 13 '15 at 16:48