# Issue with solving coupled partial differential equations with NDSolve for certain boundary conditions

I am trying to solve the coupled PDEs for functions u, v and w. The main problem I am facing is that I am unable to understand how to input boundary condition of this type in NDSolve: D[u[x, y], y] + D[v[x, y], x] == (-u[x, y]/constant) /. y -> ymax

I read about NeumannValue in the Mathematica Documentation, but still could not realize about implementing it in my context.

This is my code:

(*Equations: *)
eqn1 := D[u[x, y], x] + D[v[x, y], y] == 0
eqn2 := \[Alpha]*D[w[x, y], x] + \[Beta] Laplacian[u[x, y], {x, y}] ==
u[x, y]
eqn3 := \[Alpha]*D[w[x, y], y] + \[Beta] Laplacian[v[x, y], {x, y}] ==
v[x, y]

(* CONSTANTS: *)
xmin := -50
xmax := 50
ymin := -5
ymax := 5

(* \[CapitalOmega] is the area we are working on*)
\[CapitalOmega] = Rectangle[{xmin, ymin}, {xmax, ymax}];

(*Boundary Conditions: *)
bcnd1 := v[x, ymax] == 0
bcnd2 := v[x,
ymin] == \[Gamma] DiracDelta[x + x0] - \[Gamma] DiracDelta[
x - x0]
bcnd3 := D[u[x, y], y] + D[v[x, y], x] == -u[x, y]/\[Delta] /.
y -> ymax
bcnd4 := D[u[x, y], y] + D[v[x, y], x] == u[x, y]/\[Delta] /.
y -> ymin

(*other boundary condition to determine w*)
bcnd5 := w[xmax, y] == w0
bcnd6 := w[-xmax, y] == w0

solN = NDSolve[{eqn1, eqn2, eqn3, bcnd1, bcnd2, bcnd3, bcnd4, bcnd5,
bcnd6} /. {\[Alpha] -> 1, \[Beta] -> 1, \[Gamma] -> 1, \[Delta] ->
1, x0 -> 4, w0 -> 0.1}, {w, u,
v}, {x, y} \[Element] \[CapitalOmega]];

(*Output: *)
(*
NDSolve::fembdnl: The dependent variable in (u^(0,1))[x,5]+(v^(1,0))[x,5]==-u in the boundary condition DirichletCondition[(u^(0,1))[x,5]+(v^(1,0))[x,5]==-u,y==5.] needs to be linear.

NDSolve::fembdnl: The dependent variable in (u^(0,1))[x,5]+(v^(1,0))[x,5]==-u in the boundary condition DirichletCondition[(u^(0,1))[x,5]+(v^(1,0))[x,5]==-u,y==5.] needs to be linear.
*)



In the end, I have to obtain the vector {u[x,y],v[x,y]} and the function w[x,y] and plot it in the region [CapitalOmega] I have spent much time to try and do this, but could not. I will be very grateful, if somebody helps me to do so.

• just like message says DirichletCondition can not be nonlinear. How does [x,ymin] == \[Gamma] DiracDelta[x + x0] - \[Gamma] DiracDelta[x - x0] supposed to be physically? These are the BC NDSolve seems to complain about. see ref/message/NDSolve/fembdnl Commented Jun 1 at 8:26
• @Nasser The differential equations describe electron flow in a rectangular device. That particular boundary condition means that there are two current probes, one for sending current, and one for receiving. Current is injected at -x0 of strength gamma, and the other probe is at x0. I think the problem lies with the boundary conditions 3 and 4 (slip-boundary conditions). These are necessary to describe the fluid flow as these are proportional to the off-diagonal terms of the stress tensor in context of a fluid. The constant Delta is sort of a 'slip-length' from the boundaries. Commented Jun 1 at 12:09
• maybe doesn't like the deltas? Could try making those narrow gaussian functions. Commented Jun 1 at 12:35
• @neesh With the current data solution seems to be unstably. Have you any estimation for u? Commented Jun 1 at 15:13
• @BeauGeste I tried changing it, but I still had the same output, which says that it does not accept the boundary conditions 3 and 4, which are very necessary to model the correct behaviour I have shared the important part of the code here(as it did not allow me to write more than the character limit), rest is the same as mentioned in the question. bcnd2 := v[x, ymin] == \[Gamma] 1/(Sqrt[2 Pi]*s) Exp[-(x + x0)^2/(2 s^2)] - \[Gamma] 1/(Sqrt[2 Pi]*s) Exp[-(x - x0)^2/(2 s^2)] /. s -> 0.01  Commented Jun 2 at 5:38

To reproduce Figure 3 from the paper Non-local transport and the hydrodynamic shear viscosity in graphene we use FDM code from here.

 (*CONSTANTS:*)rule = {\[Alpha] -> 1, \[Beta] -> 1, \[Gamma] ->
1, \[Delta] -> 1, x0 -> 1, w0 -> .0, sigma -> 2,
p0 -> 0.}; xmin = -5;
xmax = 5;
ymin = -.5;
ymax = .5;

XYgrid[dom_List, pts_List] :=
N@Range[Sequence @@ #1, Abs[Subtract @@ #1]/#2] &, {dom,
pts - 1}];
BoundaryIndex[xgridlen_, ygridlen_] :=
Module[{tmp, left, right, bot, top},
tmp = Table[(n - 1) ygridlen + Range[1, ygridlen], {n, 1,
xgridlen}]; {left, right} = tmp[[{1, -1}]]; {bot, top} =
Transpose[{First[#], Last[#]} & /@ tmp]; {top, right[[2 ;; -2]],
bot, left[[2 ;; -2]]}];
FDMat[deriv_, xygrid_, difforder_] :=
Map[NDSolveFiniteDifferenceDerivative[#, xygrid,
"DifferenceOrder" -> difforder]["DifferentiationMatrix"] &, deriv]
{domain, pts, difforder} = {{{-5, 5}, {-.5, .5}}, {100, 10}, 4};
xygrid = XYgrid[domain, pts]; {nx, ny} =
Map[Length, xygrid]; {top, right, bot, left} =
BoundaryIndex[nx, ny]; {dx, dy, dx2, dy2, dxy} =
FDMat[{{1, 0}, {0, 1}, {2, 0}, {0, 2}, {1, 1}}, xygrid,
difforder]; boundaries = Join[top, right, bot, left]; sgrid =
Flatten[Outer[List, Sequence @@ xygrid], 1];
Attributes[MakeVariables] = {Listable};
MakeVariables[var_, n_] := Table[Unique[var], {n}];
{u, v, w, p} = MakeVariables[{uu, vv, ww, pp}, nx ny];

sx1 = Table[
If[Abs[sgrid[[bot]][[i, 1]] + x0] <= 1/20, 1, 0] /. rule, {i,
Length[bot]}]; sx2 =
Table[If[Abs[sgrid[[bot]][[i, 1]] - x0] <= 1/20, 1, 0] /. rule, {i,
Length[bot]}];

eqn1 = (dx2 + dy2) . p - dx . u - dy . v; eqn4 = (dx2 + dy2) . w;
eqn2 = \[Alpha]*dx . w + \[Beta] (dx2 + dy2) . u - u;
eqn3 = \[Alpha]*dy . w + \[Beta] (dx2 + dy2) . v - v;
(*bcnd1:*)eqn3[[top]] = v[[top]];
(*bcnd2:=*)eqn3[[bot]] = v[[bot]] - \[Gamma] sx1 + \[Gamma] sx2;
(*bcnd3:*)eqn2[[top]] = (dy . u + dx . v)[[top]] + u[[top]]/\[Delta];
(*bcnd4:*)eqn2[[bot]] = (dy . u + dx . v)[[bot]] - u[[bot]]/\[Delta];
(*bcnd5:=*)eqn1[[left]] = p[[left]] - p0;
(*bcnd6:=*)eqn1[[right]] = p[[right]] - p0;
eqn1[[top]] = (dy . p)[[top]]; eqn1[[bot]] = (dy . p)[[bot]];
eqn2[[left]] = (dx . u)[[left]]; eqn2[[right]] = (dx . u)[[right]];
eqn3[[left]] = (dx . v)[[left]]; eqn3[[right]] = (dx . v)[[right]];
eqn4[[left]] = (dx . w)[[left]] - w0;
(*bcnd6:=*)eqn4[[right]] = (dx . w)[[right]] - w0;
eqn4[[top]] = (dy . w)[[top]];
eqn4[[bot]] = (dy . w)[[bot]] -
sigma (\[Gamma]  sx1 - \[Gamma]  sx2); eqs =
Join[eqn1, eqn2, eqn3, eqn4]; var = Join[u, v, w, p]; {vec, mat} =
CoefficientArrays[eqs /. rule, var];
sol = LinearSolve[mat, -vec]; ns = Length[sol]; phi =
Take[sol, -ns/4]; sol1 =
sol - Join[dx . phi, dy . phi, 0  phi, 0  phi];

{U, V, W, P} =
Map[Interpolation@Join[sgrid, Transpose@List@#, 2] &,
Partition[sol1, Length[sgrid]]];


Visualization

Show[DensityPlot[W[x, y], {x, -2, 2}, {y, -.5, .5},
PlotLegends -> Automatic, ColorFunction -> "Rainbow",
PlotRange -> Automatic, AspectRatio -> Automatic,
ColorFunctionScaling -> True],
StreamPlot[
Evaluate[-Grad[W[x, y], {x, y}]], {x, -2, 2}, {y, -.5, .5},
AspectRatio -> Automatic, StreamColorFunction -> None,
StreamStyle -> LightGray, PlotLegends -> Automatic,
PlotRange -> Automatic]]

Show[DensityPlot[W[x, y], {x, -2, 2}, {y, -.5, .5},
AspectRatio -> Automatic, ColorFunction -> "Rainbow",
PlotLegends -> Automatic, PlotRange -> Automatic],
StreamPlot[-Evaluate[{U[x, y], V[x, y]}], {x, -2, 2}, {y, -.5, .5},
AspectRatio -> Automatic, StreamColorFunction -> None,
StreamStyle -> LightGray, PlotLegends -> Automatic,
PlotRange -> Automatic]]


• Thank-you soo much for this code @Alex Trounev You just solved my entire problem, thanks. I will try to understand the code. Thanks man! You're a saviour. Commented Jun 2 at 16:24
• @neesh Please note that code has been edited. Figure 3b added. Commented Jun 2 at 16:44
• Dear user @Alex Trounev I tried to change the geometry of the problem. I wanted the current probes to be on the opposite edges facing each other. Hence, the boundary conditions will change as: bcnd1 = v[x, ymax] == \[Gamma] DiracDelta[x - x0] ; bcnd2 = v[x, ymin] == - \[Gamma] DiracDelta[ x - x0] ; Commented Jun 14 at 18:07
• Follow-up of previous comment Hence, I changed the following in your code to try to simulate this: sx1 = Table[ If [ Abs [sgrid[[top]] [[i,1]] - x0] <= 1/20, 1, 0] /. rule, {i, Length[top]}]; sx2 = Table[ If [ Abs [sgrid[[bot]] [[i,1]] - x0] <= 1/20, 1, 0] /. rule, {i, Length[bot]}]; eqn3[[top]] = v[[top]] - \[Gamma] sx1; eqn3[[bot]] = v[[bot]] + \[Gamma] sx2; eqn4[[top]] = (dy . w)[[top]] - sigma \[Gamma] sx1 ; eqn4[[bot]] = (dy . w)[[bot]] + sigma \[Gamma] sx2 ;` However, the result I am getting is weird. I'm sorry to ask again. I have tried many times, but no avail. Commented Jun 14 at 18:08
• Follow up of previous comment: In particular, I am trying to see the variation of ΔΦ with respect to changing β. Here, ΔΦ is Φ[xmax/3 , ymax] - Φ[xmax/3 , ymin]; Ideally, it should be negative in some range, and positive in some range. But I am getting a constant result. Commented Jun 14 at 18:32