# Neumann boundary conditions are not satisfied

I am solving a system of 3 ODEs and I want to implement neumann BC (d/dx = 0 @x=0). However, the solution does not seem to obey neumann BC when parameter A1/A2/A3 is non zero, as I notice that the solution curves at x=0 is not flat.

List of constant parameters:

s = 1; (*ratio between Kc and Kb*)
p = 5; (*ratio between Db and Da*)
q = 1; (*ratio between Dc and Da*)
u = 0.5; (*ratio between surface concentration of A and initial \
concentration of active sites*)
TM = 1; (*Thiele modulus*)
Rba = 0.6; (*ratio between surface concentration of species B and A*)
\

Rca = 0.5; (*ratio between surface concentration of species C and A*)
\

Kb = 1.38*10^-23;  (*boltzman constant*)
T = 25 + 273; (*temperature*)
(*coefficients for function of non-ideal chemical potential*)
A1 = 1;
A2 = 1;
A3 = 1;
B1 = 1;
B2 = 1;
B3 = 1;
C1 = 1;
C2 = 1;
C3 = 1;


Some intermediate expressions:

(*Non-ideal chemical potential function of each species*)
chemicalpotential1[x_] =
A1/2*(1 - Tanh[B1 (x - C1)])*Kb*T; (*for species A*)
chemicalpotential2[x_] =
A2/2*(1 - Tanh[B2 (x - C2)])*Kb*T; (*for species B*)
chemicalpotential3[x_] =
A3/2*(1 - Tanh[B3 (x - C3)])*Kb*T; (*for species C*)

(*Expressions that involve non-ideal chemical potential function in \
the first order differential*)
diffintermediate1[x_] =
a[x]*x^2*chemicalpotential1'[x];(*for species A*)
diffintermediate2[x_] =
b[x]*x^2*chemicalpotential2'[x];(*for species B*)
diffintermediate3[x_] = (phi[x] - 1)*x^2*
chemicalpotential3'[x];(*for species C*)


The system of ODEs and BCs:

ODE1 = a''[x] + 2/x*a'[x] + 1/(Kb*T)/x^2*diffintermediate1'[x] -
TM^2*(1 + s)*phi[x]*a[x] == NeumannValue[0, x == 0]; ODE2  =
b''[x] + 2/x*b'[x] + 1/(Kb*T)/x^2*diffintermediate2'[x] +
TM^2/p*phi[x]*a[x] == NeumannValue[0, x == 0];
ODE3  = phi''[x] + 2/x*phi'[x] + 1/(Kb*T)/x^2*diffintermediate3'[x] -
s*TM^2*u/q*phi[x]*a[x] == NeumannValue[0, x == 0]; ODES = {ODE1,
ODE2 , ODE3}; (*Grouping the ODEs for input in NDSolveValue \
function*)

(*Dimensionless dirichlet BCs for species A, B, and C*)
bcs = {  a[1] == 1 , b[1] == Rba, phi[1] == 1 - u*Rca}


Solving and plotting the ODEs,

{asol, bsol, phisol} =
NDSolveValue[{ODES, bcs}, {a, b, phi}, {x, 0, 1},
Method -> "FiniteElement"];
csol[x_] = (1 - phisol[x])/u; (*transforming phi into dimensionless c*)

(*Plotting the concentration profiles*)
Plot[{asol[x], bsol[x], phisol[x], csol[x]}, {x, 0, 1},
PlotLegends -> "Expressions"]


Here is the code for testing the gradient at x=0, I realize that gradient@x=0 is more away from 0 when the A1/A2/A3 has larger magnitude, which oppose what I want for neumann BC.

(*Testing the gradient at x=0*)
differentiationdsd[x_] = phisol'[x];
differentiationdsd[0]

• Well, the derivative phisol is almost zero. I get -0.0302489. Since this is numerical method and not exact, may be trying to decrease the mesh spacing or other ways to improve the FEM accuracy could make it closer to zero? May 17 at 7:27
• But I agree that it does look a bit off zero, even for numerical method. But may be trying to improve the FEM accuracy could make this error smaller? Otherwise, it could be a bug. Just speculating. I am sure someone here with more knowledge of FEM in Mathematica will know for sure. May 17 at 7:47
• I'd like to remember everyone that the Neumann condition in FEM is imposed a bit... counterintuitively. Part of that is to interpret $\frac{\partial u}{\partial \nu}$ of a function $u \in H^1(\varOmega)$ on a smooth domain $\varOmega$ as an element in $H^{-1/2}(\partial \varOmega)$. Sloppily speakin, one derivative is lost due to taking the normal derivative, another one is lost due to restriction to the boundary $\partial \varOmega$ which has codimension one. May 17 at 11:53
• So it cannot be expected that $\frac{\partial u_h}{\partial \nu}$ of the discrete solution $u_h$ converges in $L^2$ or $L^\infty$ to $0$ under grid refinement $u_h \to 0$. Nonetheless, $u_h$ can still converge in $H^1$ to the true solution. May 17 at 11:56
• @HenrikSchumacher Oh this is a bit over my head… any recommended reference? May 17 at 12:15

Not sure how to explain the non-convergence of phisol'[0]. (Though phisol'[0] == -0.0302716 is close to zero, it won't get closer to zero with smaller MaxCellMeasure. This might be related to the interpolation method of FiniteElement. ) In this answer I just want to point out that

… the solution does not seem to obey neumann BC when parameter A1/A2/A3 is non zero, as I notice that the solution curves at x=0 is not flat

is more of an illusion. The solution found by FiniteElement method is reliable, I'll verify this by solving the problem with finite difference method (FDM).

I'll use pdetoae for generation of finite difference equations.

varlst = {a, b, phi};
newode = Map[# x &, ODES, {2}] /. NeumannValue[__] :> 0 // Expand
bcadd = Table[var'[0] == 0, {var, varlst}];

domain = {0, 1};
points = 25;
grid = Array[# &, points, domain];
difforder = 2;
(* Definition of pdetoae isn't included in this post,
ptoafunc = pdetoae[varlst[x], grid, difforder];
del = #[[2 ;; -2]] &;

ae = del /@ ptoafunc@newode;

MapThread[(guess[#][x_] = #2) &, {varlst, bcs[[All, -1]]}];

solrulelst =
FindRoot[Flatten@{ae, aebc},
Flatten[#, 1] &@Table[{var[x], guess[var][x]}, {var, varlst}, {x, grid}]];

solmatlst = ArrayReshape[solrulelst[[All, -1]], {3, points}];

solfunclst =
ListInterpolation[#, {grid}, InterpolationOrder -> difforder] & /@ solmatlst;

solfunclst[[3]]'[0]
(* -3.9968*10^-15 *)


As we can see, with 25 grid points and 2nd order difference formula, the obtained $$\phi'(0)$$ is very close to zero, while these solutions are visually the same as those obtained with FiniteElement:

ref = Plot[{asol[x], bsol[x], phisol[x]}, {x, 0, 1},
PlotLegends -> "Expressions"];
Show[ListPlot@solfunclst, ref]


Here I've made use of an undocumented syntax of ListPlot. See this post for more info.