what causes the error “The dependent variable in the boundary condition needs to be linear” when using NDSolve?

Question

V 12.1.1 on windows 10

Why the following works

ClearAll[u, x, y];
pde  = Laplacian[u[x, y], {x, y}];
bc   = {u[0, y] == 1, u[1, y] == 0};

sol  = NDSolve[{pde == NeumannValue[0, y == 0] + NeumannValue[0, y == 1], bc},
          u, {x, 0, 1}, {y, 0, 1}]

Plot3D[Evaluate[u[x, y] /. sol], {x, 0, 1}, {y, 0, 1},  AxesLabel -> {"x", "y", "u"}]

ps. I know NeumannValue[0, y == 0] above are not needed in this example, since by default they are zero. So the above could also be written as

sol = NDSolve[{pde == 0, bc}, u, {x, 0, 1}, {y, 0, 1}]

But here is the problem.

Now I just changed NeumannValue[0, y == 0] to use standard Derivative instead, and NDSolve gives now an error

ClearAll[u, x, y];
pde = Laplacian[u[x, y], {x, y}] == 0;
bc  = {u[0, y] == 1, u[1, y] == 0, Derivative[0, 1][u][x, 0] == 0, 
          Derivative[0, 1][u][x, 1] == 0};

sol = NDSolve[{pde, bc}, u, {x, 0, 1}, {y, 0, 1}]

But bc is not nonlinear. I have not even changed it. So why does the error says it is?

To show these boundary conditions are valid, changed NDSolve to DSolve and it gave same solution as the first one above using NeumannValue

ClearAll[u,x,y];
pde = Laplacian[u[x,y],{x,y}]==0;
bc  = {u[0,y]==1,u[1,y]==0,Derivative[0,1][u][x,0]==0,Derivative[0,1][u][x,1]==0};

sol = DSolve[{pde,bc},u[x,y],{x,y}];
sol = Simplify[Activate[sol]]

Plot3D[u[x,y]/.sol,{x,0,1},{y,0,1}]

Plot3D[u[x, y] /. sol, {x, 0, 1}, {y, 0, 1}]

Can't one use Derivative now with NDSolve and have to use NeumannValue? And why it says dependent variable in the boundary condition is not linear when it is?

ps. I tested this on 12 and 11.3 and they give same error. I do not have earlier versions of Mathematica to test on.

Answer 1

OK, let me extend my comment to an answer. You're solving Laplace equation, for which the new-in-v10 FiniteElement is the only available method in NDSolve at the moment. For more info about when FiniteElement will be automatically chosen, check the following post:

PDEs : automatic method choice : TensorProductGrid or FiniteElement?

When FiniteElement is used for spatial discretization, one cannot express b.c. (include but not limited to Neumann b.c. and Robin b.c. ) with Derivative, and has to use NeumannValue, at least now. [I tried to convince user21 to implement the automatic transform from Derivative to NeumannValue for regular domain (rectangle, disk, cube, ball, etc.) in comments under this post, but failed. ]

BTW, the following are examples solving Laplace equation in pre-v10 days:

Numerically solving an inhomogeneous partial differential equation

Optimizing a Numerical Laplace Equation Solver

Answer 2

You assume that NeumannValue and Derivative are interchangeable; they are generally not. This is explained in more detail in the section The Relation between NeumannValue and Boundary Derivatives in the Finite Element Method User Guide. Also the NeumannValue ref page has detailed information what it models.

Let's look at the message: The message states that Derivative[0,1][u][x,0] needs to be linear when used in DirichletConditon:

This means that NDSolve parsed your input Derivative[0,1][u][x,0] as a DirichletCondition and it tells you that it needs to be linear. You get the same message when you use

bc = {u[0, y] == 1, u[1, y] == 0, Sin[u[x, 0]] == 0,...}

Things like f[u[x,y]==... where f is not some sort of Indentity function are nonlinear and that includes Derivative[0,1][u].

Also, note that the message tag is fembdnl. So NDSolve choose to use the finite element method (all messages coming from the FEM subsystem start with fem).

What confuses me is that you show a screenshot of the message ref page but it seems you have not read the page completely. Because there is an example on that page that shows exactly this issue. Let me quote that (ref/message/NDSolve/fembdnl#506322725):

---

For the finite element, spatial discretization derivatives need to be modeled with NeumannValue:

NDSolve[{Laplacian[c[x, y, z], {x, y, z}] - c[x, y, z] == 0, 
  c[x, y, 0] == 1, Derivative[0, 0, 1][c][x, y, 1] == c[x, y, 1], 
  Derivative[0, 1, 0][c][x, 1, z] == c[x, 1, z], 
  Derivative[1, 0, 0][c][1, y, z] == c[1, y, z], 
  Derivative[0, 1, 0][c][x, 0, z] == c[x, 0, z], 
  Derivative[1, 0, 0][c][0, y, z] == c[0, y, z]}, c, {x, 0, 1}, {y, 0,
   1}, {z, 0, 1}]

The solution is to use NeumannValue:

NDSolveValue[{Laplacian[c[x, y, z], {x, y, z}] - c[x, y, z] == 
   NeumannValue[c[x, y, z], 
    z == 1 || y == 1 || x == 1 || y == 0 || x == 0], 
  DirichletCondition[c[x, y, z] == 1, z == 0]}, c, {x, 0, 1}, {y, 0, 
  1}, {z, 0, 1}]

---

I have moved this example further up - now it's the second example. I also added links to the relation between Derivative and NeumannValue section and added a link to the section: What Triggers the Use of the Finite Element Method

The fact that you can replace NDSolve with DSolve in your example is a coincidence and does not hold generally when making use of NeumannValue. Looking at it from the other side: You can not use your first example (with NeumannValue`) and use DSolve on it either.

So the answer to your question:

Can't one use Derivative now with NDSolve and have to use NeumannValue

You can not use Derivative as a boundary condition with the finite element method. This never worked so it is nothing that changed "now".

Concerning the comments: The automatic conversion of Derivative to NeumannValue is far from trivial. It's not impossible but seen that this has not been requested as a Future enhancements for the finite element method I am not sure it has a high priority.