2

For my project I am studying a paper, namely "Perturbation solutions for the nonlinear Poisson–Boltzmann equation with a higher order-accuracy Debye–Huckel approximation" by Cunlu Zhao, Qiuwang Wang and Min Zeng, Zeitschrift für angewandte Mathematik und Physik volume 71, 140 (2020), DOI: 10.1007/s00033-020-01367-9.
I am stuck at a problem: I tried to solve a nonlinear P-B equation in an unbounded domain, in order to determine the EDL potential distributions around a spherical particle, by using the finite difference scheme given in that paper. But I found I am unable to figure out how specify boundary conditions on an infinite domain. Precisely the equation is the following $$ \frac{1}{r^2}\frac{d}{dr}\bigg(r^2 \frac{d\psi}{dr}\bigg)=\sinh(\psi), $$ while the boundary conditions are $$ \begin{cases} \psi|_{r=k}=\zeta\\ \psi|_{r\rightarrow\infty}=0 \end{cases}. $$

I am not figuring out how to construct the grid points in an infinite domain and solve it. Please help me dealing with this equation. Thanks in advance!

Daniele Tampieri
  • 10,752
Deepak Gupta
  • 21
  • When you say you want to construct grid points does that mean you want to solve this numerically ? If so : restrict yourself to a finite and large enough domain. – Kurt G. Apr 01 '22 at 09:28
  • What do you mean by finite large enough domain since ultimately I have to solve by some fix large number so how can we use boundary condition which is at infinity. – Deepak Gupta Apr 01 '22 at 09:36
  • Please answer first my question (yes or no) if your task is to solve this equation numerically . – Kurt G. Apr 01 '22 at 09:50
  • Aparently it is (given the title). So : in that case there is no way to construct grid points for the infinite domain . Choose a large upper bound , choose finite grid points and experiment with those parameters until your solution doesn't change anymore . Done. – Kurt G. Apr 01 '22 at 09:57
  • Yes...! I have to solve it numerically. Thank you I understood what you want to say but I have one more doubt is there that how would we take care of our 2nd boundary condition ? Cn you please tell me the steps in details. – Deepak Gupta Apr 02 '22 at 05:32
  • Is the "2nd boundary condition" $\psi|_{r=k}=\zeta$ ? If so, what exactly is the problem there to "take care" of it? – Kurt G. Apr 02 '22 at 13:37
  • No.. this is the first 2nd is at infinity if we do consider that condition then we are unable to solve system of Linear since in that case there will be n+1 variable and n equation so we can not find the solution – Deepak Gupta Apr 05 '22 at 13:44

2 Answers2

1

For the case of a single flat plate, the potential varies as $\psi\propto \exp(-\kappa y)$, both for small applied surface potential ("Debye-Hückel theory") and for the full nonlinear Poisson Boltzmann equation (Gouy-Chapman). In the paper you mentioned, the authors have nondimensionalized all lengths by the Debye length $\kappa^{-1}$. Hence, $\psi(\bar{y})\propto \exp(-\bar{y})$. This means that already at $\bar{y}=5$, say, the potential is $\exp{(5)}\approx 148$ times smaller than at the electrode at $\bar{y}=0$ where the potential is applied. So presumably you will see that the difference between applying your second boundary condition either at $\psi(\bar{y}=5)=0$ and $\psi(\bar{y}=6)=0$ are very small. As Kurt G. says:

Choose a large upper bound , choose finite grid points and experiment with those parameters until your solution doesn't change anymore . Done.

The case of a spherical particle will be very similar. Chose a certain dimensionless sphere radius $K$ and apply the second boundary condition at $r=K+10$. You will see that the potential at $K+5$ will already be close to zero. Now set the second boundary condition at $r=K+6$ and observe that the differences between these two solution will be small.

adriaanJ
  • 649
0

For large $r$ the equation gets close to $$ ψ''=\sinhψ\implies ψ'^2=2\coshψ+C $$ and with the asymptotic convergence at infinity $C=-2$ so that $$ ψ'^2-(2\sinh(ψ/2))^2=0. $$ The desired stable behavior for the far-field solution corresponds to the factor $$ψ'+2\sinh(ψ/2)=0.$$ This can be taken as upper boundary condition at some largish but finite upper interval end.

Let's try an example with $ψ(0.5)=1$.

def sys(r,u):
    ψ,θ = u
    return θ, sinh(ψ)-2*θ/r

def bc(ua,ub):
    ψa,θa = ua
    ψb,θb = ub
    return ψa-1, θb+2*sinh(ψb/2)


a,b = 0.5, 5 
r = a+logspace(-3,0,10)(b-a)
u = [1-r/b, -1/b+0r]


res = solve_bvp(sys,bc,r,u)

This works and gives a simple falling and asymptotic curve

enter image description here

Lutz Lehmann
  • 131,652