9
$\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.

Improve this question
$\endgroup$
4
  • 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, 2014 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, 2014 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, 2014 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, 2014 at 13:03

1 Answer 1

Reset to default
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}]
Improve this answer
$\endgroup$
5
  • $\begingroup$ Thanks for confirming the problem. What about the minus sign inside the Grad? Will be supported? $\endgroup$
    – unlikely
    Jan 16, 2015 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, 2015 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, 2015 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, 2015 at 20:49
  • $\begingroup$ @unlikely, agreed and that's something I am working on. $\endgroup$
    – user21
    Feb 14, 2015 at 14:07

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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