# System of PDEs with coupled Neumann BC--how to use NeumannValue?

I have the following system of PDEs for functions $$u_1,u_2$$ in a region $$\Omega\subset\mathbb R^2$$: $$\begin{cases} 2\partial_{xx}u_1+\partial_{yy}u_1+\partial_{xy}u_2-u_1=f_1 & (x,y)\in\Omega\\ \partial_{xx}u_2+2\partial_{yy}u_2+\partial_{xy}u_1-u_2=f_2 & (x,y)\in\Omega\\ \nu\cdot\nabla u_1=\nu_x-\nu\cdot\partial_x(u_1,u_2) & (x,y)\in\partial\Omega\\ \nu\cdot\nabla u_2=\nu_y-\nu\cdot\partial_y(u_1,u_2) & (x,y)\in\partial\Omega, \end{cases}$$ where $$f_1,f_2$$ are known functions on $$\Omega$$ and $$\nu$$ is the outward unit normal to the boundary $$\partial\Omega$$.

I am interested in solving this problem numerically using Mathematica for generic $$f_1,$$ $$f_2,$$ and $$\Omega$$, but for ease of presentation, I here is my attempt at a solution when $$f_1(x,y)=x$$, $$f_2(x,y)=y$$, and $$\Omega$$ is the unit disk:

lhs1 = 2 D[u1[x, y], {x, 2}] + D[u1[x, y], {y, 2}] + D[u2[x, y], x, y] - u1[x, y];
lhs2 = 2 D[u2[x, y], {y, 2}] + D[u2[x, y], {x, 2}] + D[u1[x, y], x, y] - u2[x, y];
f1 = x;
f2 = y;
bc1 = NeumannValue[x - {x, y} . D[{u1[x, y], u2[x, y]}, x], True];
bc2 = NeumannValue[y - {x, y} . D[{u1[x, y], u2[x, y]}, y], True];
NDSolveValue[{lhs1 == f1 + bc1, lhs2 == f2 + bc2}, {u1,u2}, {x,y} \[Element] Disk[]];


Mathematica doesn't seem to like derivatives appearing in Nuemann boundary values; it gives the error:

NDSolve::fembcdepderiv: Derivatives of dependent variables in boundary conditions are not supported with the Finite Element Method in this version of NDSolve.


I do not see an analytical way to decouple the Neumann boundary conditions or transform the boundary conditions to no longer include derivatives on the right hand side. Are there any numerical workarounds?

• You misunderstand the usage of NeumannValue. It's not because the Neumann value is coupled, it's because Derivative cannot appear inside NeumannValue, see discussions in e.g. mathematica.stackexchange.com/a/224815/1871 mathematica.stackexchange.com/a/267744/1871 mathematica.stackexchange.com/a/245309/1871 But sadly knowing the correct usage of NeumannValue isn't enough to resolve the problem, because there doesn't seem to be an obvious way to express your b.c. as legal NeumannValue. Commented Aug 18, 2022 at 15:49
• @xzczd Right, that is the root of the problem. My understanding $u_1$ and $u_2$ can appear in NeumannValue but not their derivatives. Commented Aug 18, 2022 at 16:51
• Just to confirm, $\nu_x=x,\ \nu_y=y$, right? Commented Aug 19, 2022 at 5:03
• @xzczd In the case that $\Omega$ is the unit disk, yes. Commented Aug 19, 2022 at 14:01