5
$\begingroup$

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]
Improve this question
$\endgroup$
5
  • 2
    $\begingroup$ 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? $\endgroup$ Commented May 17, 2022 at 7:27
  • 1
    $\begingroup$ 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. $\endgroup$ Commented May 17, 2022 at 7:47
  • $\begingroup$ 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. $\endgroup$ Commented May 17, 2022 at 11:53
  • $\begingroup$ 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. $\endgroup$ Commented May 17, 2022 at 11:56
  • $\begingroup$ @HenrikSchumacher Oh this is a bit over my head… any recommended reference? $\endgroup$ Commented May 17, 2022 at 12:15

1 Answer 1

Reset to default
7
$\begingroup$

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, please find it in the link above. *) ptoafunc = pdetoae[varlst[x], grid, difforder]; del = #[[2 ;; -2]] &; ae = del /@ ptoafunc@newode; aebc = ptoafunc@{bcadd, bcs}; 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]

enter image description here

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

Improve this answer
$\endgroup$

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.