8
$\begingroup$

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

Mathematica graphics

and

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

Mathematica graphics

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

Mathematica graphics

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

Mathematica graphics

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

Mathematica graphics

Any help undestanding this matter is appreciated.

$\endgroup$
  • 1
    $\begingroup$ 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?) $\endgroup$ – Michael E2 Dec 21 '14 at 22:52
  • $\begingroup$ @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...). $\endgroup$ – unlikely Dec 22 '14 at 8:43
  • $\begingroup$ @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. $\endgroup$ – user21 Dec 29 '14 at 16:24
  • 1
    $\begingroup$ @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. $\endgroup$ – user21 Dec 31 '14 at 13:03
3
$\begingroup$

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}], 
    Div | Grad] == 
      1 + NeumannValue[1, x == 1], u[0] == 0}, u, {x, 0, 1}, 
  Method -> "FiniteElement"]
Plot[%[x], {x, 0, 1}]
$\endgroup$
  • $\begingroup$ Thanks for confirming the problem. What about the minus sign inside the Grad? Will be supported? $\endgroup$ – unlikely Jan 16 '15 at 10:54
  • $\begingroup$ Can you eventually add a skeleton of the current best/safest way to input the equation and Neumann boundary condition? $\endgroup$ – unlikely Jan 16 '15 at 10:55
  • $\begingroup$ @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. $\endgroup$ – user21 Feb 13 '15 at 16:48
  • $\begingroup$ I don't think it's necessary to support every possible placement of a scalar multiplier of any PDE coefficent. But, I think it would be useful to expand a bit the documentation to explain what are the "supported" way(s) to input a PDE with Inactive differential operators: after many readings I wasn't able to figure the proper way and what was wrong with my approach. $\endgroup$ – unlikely Feb 13 '15 at 20:49
  • $\begingroup$ @unlikely, agreed and that's something I am working on. $\endgroup$ – user21 Feb 14 '15 at 14:07

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.