# What causes the warning "The dependent variable in the boundary condition needs to be linear" when using NDSolve?

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 a warning

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 warning 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 warning. I do not have earlier versions of Mathematica to test on.

• I'm pretty sure this issue has raised up in a number of posts, but surprisingly I can't find a duplicate. Seems that it's never asked in a clear enough way as in your question. In short, When FiniteElement is used for spatial discretization, one cannot express Neumann and Robin b.c. with Derivative, and have to use NeumannValue, at least now. You might remember the discussion here: mathematica.stackexchange.com/questions/172972/… Commented Jun 27, 2020 at 11:28
• @xzczd thanks. I remember that question. But that was about using finite elements to solve time dependent ode. I am here not even asking to use Finite elements in the example where the error was generated. I see you said When "FiniteElement" is chosen. So what what is the bottom line here? One can not use Derivative at all as BC. with NDSolve Since when this started? May be if you know, you can make an answer giving the rules of thumb to follow, as I am now confused :) Commented Jun 27, 2020 at 11:39

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

I create a function allowfemdbc to facilitate setting Neumann/Robin b.c. when FiniteElement is chosen.

• Ok, thanks. this makes it more clear. I wonder why can't NDSolve, internally modify the Derivative to NeumannValue? This way, same code used for DSolve can be used with NDSolve without having to change the bc. This also assumes the user has to know that FEM will be used to solve the equations they are working on. Not every one knows these internal things by just looking at the equations and know what solver Mathematica choose to use in order to change the BC to make it happy. Commented Jun 27, 2020 at 12:43
• @Nasser You may take part in the discussion in 2nd link :) . Commented Jun 27, 2020 at 13:05

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.