# Boundary conditions being ignored or not obeyed by NDSolveValue

I have the following coupled ordinary differential equations (drift-diffusion with a loss term) that I am having trouble with solving. One of the boundary condition I supply is being ignored. When I provide a different one that should also hold, I flag an error.

mod1D[zMax_, {l0_, lb_, \[Gamma]s_}, {g_, \[Sigma]_}] := Module[{},

\[Tau]s = 1/\[Gamma]s;

(*dn/dt=0*)
pde1 = -D[D[u[z], {z, 1}] + lb/l0 v[z], {z, 1}];
(*ds/dt=0*)
pde2 = -D[D[v[z], {z, 1}] + lb/l0 u[z], {z, 1}] + v[z] -
( 2 g \[Tau]s)/(Sqrt[2 \[Pi]] l0^2 \[Sigma])
Exp[-(z^2/(2 \[Sigma]^2))];

(* Neumann bc for charge*)
Bu3 = NeumannValue[0, z == zMax];
Bu4 = NeumannValue[0, z == -zMax];

(* Dirichlet bc for spin*)
R1 = DirichletCondition[v[z] == 0, z == -zMax];
R2 = DirichletCondition[v[z] == 0, z == zMax];

(* Dirichlet bc for charge*)
R3 = DirichletCondition[u[z] == -u[-z], z == zMax];
R4 = DirichletCondition[u[z] == 0, z == 0];

(*Solve the PDE over the region*)
sol = NDSolveValue[{pde1 == Bu3 + Bu4, pde2 == 0, R1, R2,
R3}, {u, v}, {z, -zMax, zMax},
Method -> {"FiniteElement",
"MeshOptions" -> {"MaxCellMeasure" -> 0.0001}},
MaxStepSize -> 0.005];

]


mod1D[10, {1, 2, 1}, {1, 0.001}]
Plot[sol[[1]][z], {z, -10, 10}, PlotStyle -> {{Black, Thick}},
PlotRange -> {{-10, 10}, {-2, 2.1}},
FrameLabel -> {"y", "z", "u"},
PlotLabel -> "Steady-state solution for u"]


gives

which is not an odd function despite the boundary condition R3. If I use the boundary condition R4 instead, I get an error NDSolveValue::bcnop: No places were found on the boundary where z==0 was True, so DirichletCondition[u==0,z==0] will effectively be ignored. However, I know, from analytic solutions, solving same problem in Python using solve_bvp, and solving this same problem in Mathematica with an added dimension, that u[0] = 0.

What I expect to see if the same shape as shown above but shifted down to produce an odd function of z. I get same issue with NDSolve.

Using 13.2.

• In the current version NDSolve can't accept boundary condition in a form of R3, R4 using FEM. You can try solve this problem using FDM. Commented Jun 7 at 14:32
• @AlexTrounev I assume that is due to the conditions not being on the edge of the region? to clarify, to use FDM you mean by brute force and not as an option in NDSolve? I have done FDM in Python and been happy with it for this problem. My issue lies in that I want to add Poisson equation and Python has been fruitless and was hoping NDSolve might be more successful. Commented Jun 7 at 15:14
• It could be better if you show your problem with Poisson as well. Commented Jun 7 at 16:04

To have a boundary condition within the region you need to add an include point like so:

Needs["NDSolveFEM"]

mesh = ToElementMesh[Line[{{-zMax}, {zMax}}],
"IncludePoints" -> {{0}}, "MaxCellMeasure" -> 0.0001];
(*Solve the PDE over the region*)
sol = NDSolveValue[{pde1 == Bu3 + Bu4, pde2 == 0, R1, R2, R4}, {u,
v}, {z} \[Element] mesh];


To model a periodic boundary condition use PeriodicBoundaryCondition:

mod1D[zMax_, {l0_, lb_, \[Gamma]s_}, {g_, \[Sigma]_}] :=
Module[{}, \[Tau]s = 1/\[Gamma]s;
(*dn/dt=0*)pde1 = -D[D[u[z], {z, 1}] + lb/l0  v[z], {z, 1}];
(*ds/dt=0*)
pde2 = -D[D[v[z], {z, 1}] + lb/l0  u[z], {z, 1}] +
v[z] - (2  g  \[Tau]s)/(Sqrt[
2  \[Pi]]  l0^2  \[Sigma])  Exp[-(z^2/(2  \[Sigma]^2))];
(*Neumann bc for charge*)Bu3 = NeumannValue[0, z == zMax];
Bu4 = NeumannValue[0, z == -zMax];
(*Dirichlet bc for spin*)
R1 = DirichletCondition[v[z] == 0, z == -zMax];
R2 = DirichletCondition[v[z] == 0, z == zMax];
(*Dirichlet bc for charge*)
R3 = DirichletCondition[u[z] == -u[-z], z == zMax];
R4 = DirichletCondition[u[z] == 0, z == 0];
pbc1 =
PeriodicBoundaryCondition[-u[z], z == -zMax,
Function[z, z + 2 zMax]];
pbc2 =
PeriodicBoundaryCondition[-u[z], z == zMax,
Function[z, z - 2 zMax]];
mesh =
ToElementMesh[Line[{{-zMax}, {zMax}}], "IncludePoints" -> {{0}},
"MaxCellMeasure" -> 0.0001];
(*Solve the PDE over the region*)
sol = NDSolveValue[{pde1 == Bu3 + Bu4, pde2 == 0, R1, R2, R4, pbc1,
pbc2}, {u, v}, {z} \[Element] mesh];
]

Needs["NDSolveFEM"]

mod1D[10, {1, 2, 1}, {1, 0.001}]
Plot[sol[[1]][z], {z, -10, 10}, PlotStyle -> {{Black, Thick}},
PlotRange -> {{-10, 10}, {-2, 2.1}}, FrameLabel -> {"y", "z", "u"},
PlotLabel -> "Steady-state solution for u"]