# A New Approach for the Solution of the Electrostatic Potential Differential Equations

Authors: Fadi Awawdeh; M. Khandaqji; Z. Mustafa
Journal: Mathematical Problems in Engineering
 (2009)

---

## Abstract

This paper is focused on deriving an explicit analytical solution for the prediction of the electrostatic potential, commonly used on electrokinetic research and its related applications. Different from all other analytic techniques, this approach provides a simple way to ensure the convergence of series of solution so that one can always get accurate enough approximations. This new approach can be a useful tool in electrical field applications such as the separation of a mixture of macromolecules and the removal of contaminants in soil cleaning processes.


---

## Body

## 1. Introduction

In every electrokinetic related application, the electrostatic potential is a key function to be identified in order to understand its effect on flow regimes. The methods to calculate this electrostatic potential function depend largely on approximations based on the range of values of the applied electrical field; therefore, it would be extremely useful to identify a procedure that yields a more general solution valid for a wide range of magnitude of the electrostatic potential. Thus, the same solution will be explicit for the electrostatic potential and valid for the small and large values without restrictions. To our best knowledge, such type of solution has not yet been reported in literature [1–5].It becomes very much desired that an explicit mathematical expression for the electrostatic potential should be available. This is not necessarily easy to accomplish since the conservation of charge yields a nonlinear differential equation. This is based on the fact that the equation that is the common starting point for the description of the electrostatic potential is the Poisson-Boltzmann equation [6]:
(1.1)y′′=λ2sinhy,
where y is the dimensionless electrical potential and λ is the dimensionless inverse Debye length. The nonlinear term on the right-hand side of (1.1) is related to the free charge density.The main focus of the previous studies was concentrated on how to approximate the hyperbolic function, on the right-hand side of (1.1), in order to obtain an effective analytical solution of the Poisson-Boltzmann differential equation. A very common simplification invokes the Debye-Hückel approximation is usually written as
(1.2)sinhy≈y.As a result, the nonlinear Poisson-Boltzmann equation reduces to the following linear equation:(1.3)y′′=λ2y.
Such approximation is valid for the case when -1≤y≤+1 [7]. Other contributions that extend the range of valid solution to -∞≤y≤+∞ propose splitting electrical potential values in three regions [8] as a way of simplifying the hyperbolic sine function:
(1.4)sinhy≈{-0.5exp(-y),y<-1,y,-1≤y≤1,0.5exp(y),y>+1.
Based on the Debye-Hückel approximation, a new solution strategy for the differential equations of the electrostatic potential was proposed by Oyanader and Arce [4, 9]. The authors introduced the correction function fAO to the inverse Debye length, λ. The fAO function improves the Debye-Hückel approximation. It is a recursive function of the electrical potential and has a polynomial form whose coefficients are adjusted to predict the correct values of the hyperbolic sine function.In this paper, we will suggest an alternative approach to the search for an explicit analytical solution for (1.1) by homotopy analysis method (HAM) [10–21]. In particular, a series solution is obtained and the constant coefficients are given by recursive formulas.The homotopy analysis method is based on homotopy, a fundamental concept in topology and differential geometry. Briefly speaking, by means of the HAM, one constructs a continuous mapping of an initial guess approximation to the exact solution of considered equations. An auxiliary linear operator is chosen to construct such kind of continuous mapping, and an auxiliary parameter is used to ensure the convergence of solution series. The method enjoys great freedom in choosing initial approximations and auxiliary linear operators. By means of this kind of freedom, a complicated nonlinear problem can be transferred into an infinite number of simpler, linear subproblems, as shown by Liao and Tan [19].


## 2. Approach Based on the HAM

Consider the nonlinear Poisson-Boltzmann equation:(2.1a)y′′=λ2sinhy,(2.1b)y(-1)=A,y(+1)=A,y′(0)=0.According to (2.1a) and (2.1b), it is obvious that the solution y(x) can be expressed by the base functions
(2.2){xn:n∈ℕ}
as
(2.3)y(x)=∑n=1∞bnxn,
where bn is a coefficient. It provides the Solution Expression of y(x).Based on (2.1a), we define a nonlinear operator
(2.4)N[ϕ(x;q)]=∂2ϕ(x;q)∂x2-λ2sinh[ϕ(x;q)],
where q∈[0,1] denotes the embedding parameter. Let y0(x) denote an initial guess of the exact solution y(x) which satisfies the boundary conditions (2.1b), h≠0 an auxiliary parameter, and L an auxiliary linear operator. All of y0(x), L, and h will be chosen later with great freedom. Then, we construct a one-parameter family of differential equations
(2.5)(1-q)L[ϕ(x;q)-y0(x)]=qhN[ϕ(x;q)]
subject to the boundary conditions
(2.6)ϕ(-1;q)=A,ϕ(+1;q)=A,∂ϕ(x;q)∂x|x=0=0.
Obviously, when q=0, because of the property L(0)=0 of any linear operator L, (2.5) and (2.6) have the solution
(2.7)ϕ(x;0)=y0(x),
and when q=1, since h≠0, (2.5) and (2.6) are equivalent to the original ones, (2.1a) and (2.1b), provided that
(2.8)ϕ(x;1)=y(x).
Thus, based on (2.7) and (2.8), as the embedding parameter q increases from 0 to 1, ϕ(x;q) varies continuously from the initial approximation y0(x) to the exact solution y(x). This kind of deformation ϕ(x;q) is totally determined by the so-called zeroth-order deformation equations (2.5) and (2.6).By Taylor's theorem,ϕ(x;q) can be expanded in a power series of q as follows:
(2.9)ϕ(x;q)=y0(x)+∑m=1∞ym(x)qm,
where
(2.10)ym(x)=Dm[ϕ(x;q)],Dm[f(q)]=1m!∂mf(q)∂qm|q=0.Dm is called the mth-order homotopy-derivative of f.Fortunately, the homotopy-series (2.9) contains an auxiliary parameter h. Besides, we have great freedom to choose the auxiliary linear operator L, as illustrated by Liao [16]. If the auxiliary linear parameter L and the nonzero auxiliary parameter h are properly chosen so that the power series (2.9) of ϕ(x;q) converges at q=1, then, under these assumptions, we have the so-called homotopy-series solution:
(2.11)y(x)=∑m=0∞ym(x).According to the fundamental theorems in calculus, each coefficient of the Taylor series of a function is unique. Thus,ym(x) is unique and is determined by ϕ(x;q). Therefore, the governing equations and boundary conditions of ym(x) can be deduced from the zeroth-order deformation equations (2.5) and (2.6). For brevity, define the vectors
(2.12)y⃗n(x)={y0(x),y1(x),y2(x),…,yn(x)}.Differentiating the zero-order deformation equation (2.5) m times with respect to q and then dividing by m! and finally setting q=0, we have the so-called high-order deformation equation:
(2.13)L[ym(x)-χmym-1(x)]=hℜm(y⃗m-1(x)),ym(-1)=0,ym(+1)=0,ym′(0)=0,
where
(2.14)ℜm(y⃗m-1(x))=Dm-1(N[ϕ])=1(m-1)!∂m-1N[ϕ(x;q)]∂qm-1|q=0,(2.15)χm={0,m≤1,1,m>1.
Using the definitions (2.4) and (2.14), we have the explicit expression
(2.16)Dm(N[ϕ])=Dm(ϕ′′)-λ2Dm(sinhϕ).
Now we have to find Dm(sinhϕ). It is clear from (2.10) that Dm(ϕ′′)=ym′′, and by using the definition of the operator Dm, it holds
(2.17)D0(eϕ)=ey0.
Besides, one has
(2.18)∂eϕ∂q=eϕ∂ϕ∂q.
Thus, Leibnitz's rule for derivatives of product implies that
(2.19)1m!∂meϕ∂qm=1m!∂m-1∂qm-1(eϕ∂ϕ∂q)=1m∑k=0m-11(m-k-1)!∂keϕ∂qk∂m-kϕ∂qm-k=∑k=0m-1(1-km)(1k!∂keϕ∂qk)(1(m-k)!∂m-kϕ∂qm-k).
Setting q=0 in the above expression, one has
(2.20)Dm(eϕ)=∑k=0m-1(1-km)Dk(eϕ)Dm-k(ϕ).
It is clear that D0(sinhϕ)=sinhy0, and it holds for any integer m≥1 that
(2.21)Dm(sinhϕ)=12Dm(eϕ-e-ϕ)=∑k=0m-1(1-km)Dm-k(ϕ)Dk(eϕ+e-ϕ2)=∑k=0m-1(1-km)Dm-k(ϕ)Dk(coshϕ).
Similarly,
(2.22)D0(coshϕ)=coshy0,Dm(coshϕ)=∑k=0m-1(1-km)Dm-k(ϕ)Dk(sinhϕ),m≥1.
In this line we have that
(2.23)ℜm(y⃗m-1(x))=Dm-1(N[ϕ])=ym-1′′-λ2Sm(x)
with the following recurrence formulas:
(2.24)Sn=An-Bn2,n≥0,A0=1,B0=1,An=∑m=0n-1(1-mn)Amyn-m-1,n≥1,Bn=-∑m=0n-1(1-mn)Bmyn-m-1,n≥1.
From the above recurrence formulas, we have
(2.25)ℜ1(y⃗0)=y0′′(x)-λ2y0,ℜ2(y⃗1)=y1′′(x)-λ2y1,ℜ3(y⃗2)=y2′′(x)-λ2(y2+16y03),
and so on. So, by means of symbolic computation software such as Mathematica and Matlab, it is not difficult to get ℜm(y⃗m-1(x)) for large values of m.Note that the high-order deformation equations (2.13) are linear ODEs. So, according to (2.11), the original nonlinear problem is transferred into an infinite number of linear ODEs.Both the auxiliary linear operatorL and the initial guess y0(x) are chosen under the so-called Rule of Solution Expression: the auxiliary linear operator L and the initial guess y0(x) must be chosen so that the solutions of the high-order deformation equations (2.13) exist and besides they obey the Solution Expression (2.3). So, for the solutions to obey the Solution Expression (2.3) and the boundary conditions (2.1b), we choose the initial guess of the solution:
(2.26)y0(x)=A2(1+x2).
Because the original equation (2.1a) is of second order, we simply choose such a second-order auxiliary linear operator:
(2.27)L[y]=y′′,
which has the property
(2.28)L[C1+C2x]=0
for any real coefficients C1 and C2. Let ŷm(x) denote a particular solution of (2.13), then its general solution is expressed by
(2.29)ym(x)=ŷm(x)+χmym-1(x)+C1+C2x,
where C1 and C2 are determined by the conditions in (2.13), that is,
(2.30)C2=-ŷm′(0)-χmym-1′(0),C1=0.
In this way, we get ym(x) one by one in the order m=1,2,3,….At theMth-order approximation, the solution can be expressed as follows:
(2.31)y(x)=∑m=0M∑k=0m+1αm,kxk,
where αm,k is a coefficient.


## 3. Illustrative Example

Consider the system consisting of two parallel walls, each subjected to a dimensionless electrostatic potential equal toψw. In order to take advantage of the symmetry, the coordinates have been placed in the vertical center of the capillary domain. In dimensionless coordinates, the wall located at the right-hand side is at the location ξ=1 while the one on the left-hand side is at ξ=-1 (see Figure 2). In terms of the electrostatic forces, the system is defined by (2.1a) jointly with appropriate boundary conditions [3, 4]:(3.1a)d2ψdξ2=λ2sinhψ,(3.1b)ψ(-1)=ψw,ψ(+1)=ψw,dψdξ|ξ=0=0.The HAM explicit analytical solution (2.31) depends on the auxiliary parameter h. Let Rh denote the set of all values of h which ensure the convergence of the series (2.31). According to Liao's proof [16], all of these series solutions must converge to the solution of the original problem (3.1a)-(3.1b). For a fixed ξ0 in the range of ξ, plot the h-curve, ψ(ξ0) versus h. Because the limit of all convergent series solutions (2.31) is the same for a given ψw, there exists a horizontal line segment above the region h∈Rh. So, by plotting the h-curve at a high enough order approximation, one can find an approximation Rh={h:-2<h<0}, as shown in Figure 1 in the case of ψw=1 and λ=0.24. In Table 1, the value h=-1.5 ensures the convergence of the series (2.31).Table 1
Approximation ofψ(0.2) when ψw=1, λ=0.24, h=-1.5.
Order of approximationψ(0.2)20.869040.894360.961780.9661100.9687120.9687140.9687Figure 1
Geometrical sketch of the capillary channel.Figure 2
Curveψ(0.2) versus h at the 15th order of approximation.We can conclude that, for any given value ofψw, we can always find a proper value of h to ensure the convergence of the series of the solution ψ(ξ).Let us consider the average residual errors of (3.1a):
(3.2)E=12∫-11|ψ′′-λ2sinhψ|dξ.
The average residual errors of the series (2.31) at different orders of approximations when h=-1.5 are listed in Table 2. Obviously, as the order M increases, the average residual errors decrease. This clearly indicates that our series solution (2.31) is indeed the solution of (3.1a)-(3.1b).Table 2
Residual errors of our analytic solution whenh=-1.5 at different orders of approximation.
Order of approximationE22.72×10-253.22×10-3101.43×10-4152.57×10-5208.11×10-6251.40×10-6304.13×10-7In a previous work [22], and for coordinates system as shown in Figure 1, the authors have obtained the approximate analytical solution of the model described by (3.1a)-(3.1b), and it is given by
(3.3)ψ(ξ)=ψwcoshλξcoshλ,
which is only valid in the range -1≤ψ≤+1.The proposed technique does not require small parameters in the equations under study. As a result, the technique completely eliminates the difficulty arising in the classical perturbation methods. The homotopy perturbation method (HPM) [23] was used to obtain an explicit analytical solution to the nonlinear Poisson- Boltzmann problem (3.1a)-(3.1b). The HPM solution is not accurate enough since it depends on an approximated form of (3.1a):
(3.4)d2ψdξ2=λ2(ψ+ψ3),
and it is very difficult to obtain more than the first few terms of the series solution since the constant coefficients are not given by recursive formulas.Figure3 shows the electrostatic potential variation, along the transversal coordinate ξ, predicted by the numerical solution of the complete Poisson-Boltzmann equation by Runge-Kutta method (RK4), Debye-Hückel approximation (3.3), 15th-order HAM approximation with h=-1.5, and fAO approximation [4].Figure 3
ψ versus ξ: comparison of the numerical solutions when Ψw=1, λ=0.24 and h=-1.5. Hollow dots: 15th-order HAM approximation; continued lines: numerical solution; dashed lines: Debye-Hückel approximation; dot-dashed lines: fAO approximation.In particular, Figure3 demonstrates that the use of the proposed approach, considerably improves the prediction of the electrostatic potential, ψ, yielding values nearly identical with the numerical solutions. Obviously, we dare say that the best possible approximation with the numerical solution has been achieved by using HAM solution (2.31).


## 4. Conclusion

Some models show nonlinear sources, such as, electrostatic potential. The general situation is that these models lack analytical or closed form solutions that are quite handy when coupling is required. In this study we developed a simple procedure to obtain an explicit and purely analytic solution for the prediction of the electrostatic potential; that is, the structure of the solution is known and the constant coefficients are given by recursive formulas and it is unnecessary to use any numerical methods to get any coefficients. This new approach embraces excellent results for a wide range of electrostatic potential values and can be highly adaptable to many situations.



---

*Source: 980987-2009-12-10.xml*