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.
Inactive
interacts withNDSolve
, but I would have first triedNDSolveValue[{-Inactivate[Div[Grad[u[x], {x}], {x}], Div | Grad] == <..>]
, which seems to work. (Maybe it's a bug?) $\endgroup$Inactive
equation should be equivalent to the equation obtainedActivate
-ing it (except Neumann boundary conditions...). $\endgroup$NDSolve
sees the things that are inactive (likeDiv
,Grad
) there are a few lines that then parse out the coefficients which are then put intoInitializePDECoefficients
. There simply was no rule for Times[fact_, Grad[...]] which will work in a future version, if fact isNumberQ
. For all other cases one should useDot
. Hope this helps. $\endgroup$