New approximate solutions to electrostatic differential equations obtained by using numerical and analytical methods

1 Introduction

Many problems of science and engineering lead to nonlinear differential equations. Except a few number of these equations, most of them cannot be solved analytically using traditional methods. Therefore, these problems are often handled by the most common analytical methods, such as the Adomian decomposition method [3, 8], the homotopy decomposition method [1], the optimal perturbation iteration method [4, 7, 9], the homotopy perturbation method [10], the variational iteration method [2], etc. These methods are able to deal with strongly nonlinear problems, but there is also the problem regarding the convergence region of solutions represented by the series. These regions which depend on the desired solution are generally small. Marinca et al. [17, 19, 18] introduced the Optimal Homotopy Asymptotic Method (OHAM) which is straightforward and effective for the approximate solution of nonlinear problems. The OHAM also provides us with a convenient way to control the convergence of approximation series and adjust convergence regions.

The Poisson–Boltzmann equation is a nonlinear equation that has been used for the description of the distribution of electrostatic potential in colloidal systems. One can compute the force of interactions between particles by knowing the electrostatic potential. Although it has a simple form, the Poisson–Boltzmann equation cannot be solved analytically. There have been many attempts to get numerical solutions of the nonlinear Poisson–Boltzmann equation. One such attempt was made by Hoskins and Levine [15] for interactions of identical spheres under the constant surface potential. Bowen and Sharif [5] tackled the axisymmetric situations of the problem. They used geometrical symmetry to reduce the problem in two space dimensions while it was defined in the three-dimensional space. In their paper, an adaptive finite element method was used for two interacting charged spheres enclosed in a charged cylindrical capillary. Oyanader, Arce and Dzurik [21, 22] also provided important contributions to this subject. In [20], Oyanader and Arce proposed a more effective and accurate solution for the prediction of the electrostatic potential.

For this study, we investigate the nonlinear Poisson–Boltzmann equation written as

where ψ is the electrical potential, ξ is the transversal or radial coordinate, and λ is the inverse Debye length. Note that all these quantities are dimensionless. As seen, due to the hyperbolic sine function (related to the free charge density), these two differential equations are nonlinear. Thus, we need to approximate the hyperbolic sine function in order to obtain an effective analytical solution of the Poisson–Boltzmann differential equation. According to the Debye–Huckel theory, discussed in [11], a general simplification of approximation has the form

By using (1.3), equations (1.1) and (1.2) become linear and can be easily solved. However, these solutions are valid only within the range -1≤ψ≤1. This range is very small for electrostatic potential values. In order to extend this range, one can use the following approximation:

This idea has been useful for some authors in solving the problem. However, using this approximation, they obtained a considerable number of equations, subordinated equations and restrictions, as in [16, 24].

In this paper, we apply the Optimal Homotopy Asymptotic Method (OHAM) to solve the Poisson–Boltzmann equation and to obtain an explicit mathematical expression for the prediction of the electrostatic potential. Additionally, our new results are compared with the solutions obtained by the Conjugate Gradient algorithms (CGA), see [6].

2 Optimal homotopy asymptotic method

In order to give the basic principles of the OHAM, let us consider the nonlinear differential equation

where ψ⁢(ξ) is an unknown function, ξ denotes an independent variable, g⁢(ξ) is a source function, and L,N and B are linear, nonlinear and boundary operators, respectively. First, the deformation equation is constructed as follows:

where ϕ⁢(ξ,p) is an unknown function, p is an embedding parameter, and H⁢(p) is a nonzero auxiliary function for p≠0. Clearly, at p=0 and p=1, we have ϕ⁢(ξ,0)=ψ0⁢(ξ) and ϕ⁢(ξ,1)=ψ⁢(ξ). So, as the embedding parameter p increases from 0 to 1, the solution ϕ⁢(ξ,p) deforms from the initial guess ψ0⁢(ξ) to the exact solution ψ⁢(ξ) of the original nonlinear differential equation. We can define ψ0⁢(ξ) from (2.2) for p=0:

On the other hand, for the sake of simplicity, here H⁢(p) is taken as

where C1,C2,… are the convergence-control parameters. Most recently, Herisanu et al. have proposed the generalized auxiliary function

where Hi⁢(t,Ci), i=1,2,…, are auxiliary functions. Some examples of such generalized auxiliary functions are presented in [14, 13].

Let us now consider the Taylor expansion of the solution of equation (1.2) with respect to p, that is,

Substituting (2.3) and (2.4) into (2.2) and equating the coefficients of identical powers of p equal to zero gives the linear equations

where Nm⁢(ψ0⁢(ξ),ψ1⁢(ξ),…,ψm⁢(ξ)) is the coefficient of pm in the expansion of p:

As shown previously, the convergence of series (2.4) depends on the parameters C1,C2,…. If series (2.4) is convergent at p=1, then one has

Generally speaking, the solution of equation (2.1) can be determined approximately in the form

Substituting (2.5) into (2.2), the general problem results in the following residual:

Obviously, when R⁢(ξ,Ci)=0, the approximation ψ(m)⁢(ξ,Ci) will be an exact solution. The numbers a and b are chosen to determine C1,C2,…, so that the optimal values of C1,C2,… are obtained using the method of least squares:

where R=L⁢(u(m))+g⁢(ξ)+N⁢(u(m)) is the residual and

With these constants known, one can get an approximate solution of order m. The convergence control parameters C1,C2,… can also be calculated by

where ki∈(a,b).

3 Application of OHAM to the Poisson–Boltzmann equation

In this section, we solve the nonlinear Poisson–Boltzmann equation in two different forms by using the OHAM. It should again be noted that we use the approximation

to get the solution of equation (1.4).

3.1 Capillary system in planar geometry

We first consider the capillary system in planar geometry. In order to take advantage of the symmetry, the coordinates have been placed at the vertical center of the capillary domain. In the dimensionless coordinate, the wall located on the left-hand side is at the location L=-1, while the one on the right-hand side is at L=1, as shown in Figure 1. Approximation (3.1) leads to the differential model

ψ′′=λ2⁢(ψ+ψ36),ψ⁢(1)=ψ⁢(-1)=L,ψ′⁢(0)=0.

Figure 1

Geometrical sketch of the capillary channeland the coordinate system used in Example 1.

Example 1.

The OHAM formulation of the above problem is

L⁢(ψ⁢(ξ))=ψ′′⁢(ξ)-λ2⁢ψ⁢(ξ),g⁢(ξ)=0,N⁢(ψ⁢(ξ))=-λ2⁢ψ3⁢(ξ)6,

which satisfies

(1-p)⁢[(ψ0⁢(ξ)+p⁢ψ1⁢(ξ)+p2⁢ψ2⁢(ξ)+⋯)′′-λ2⁢(ψ0⁢(ξ)+p⁢ψ1⁢(ξ)+p2⁢ψ2⁢(ξ)+⋯)]

=(pC1+p2C2+⋯)[(ψ0(ξ)+pψ1(ξ)+p2ψ2(ξ)+⋯)′′-λ2(ψ0(ξ)+pψ1(ξ)+p2ψ2(ξ)+⋯)

   -λ2(ψ0⁢(ξ)+p⁢ψ1⁢(ξ)+p2⁢ψ2⁢(ξ)+⋯)36].

By equating the coefficients of identical powers of p, we get the problems

(3.2)O⁢(p0): ⁢ψ0′′⁢(ξ)-λ2⁢ψ0⁢(ξ)=0,

(3.3)O⁢(p1): ⁢ψ1′′⁢(ξ)-λ2⁢ψ1⁢(ξ)=-λ2⁢C1⁢ψ036.

Table 1

Comparison of the absolute errors of Example 1 for different values of λ.

	OHAM solutions			CGA solutions
ξ	λ=0.5	λ=0.75	λ=1	λ=0.5	λ=0.75	λ=1
0.1	1.0523e-6	8.0569e-4	1.2154e-3	9.0863e-5	9.1669e-4	4.6952e-3
0.2	5.0526e-5	9.2134e-4	5.1028e-3	2.0547e-4	9.9904e-4	8.1245e-3
0.3	7.0879e-5	9.8362e-4	7.5216e-3	2.9056e-4	1.1365e-3	9.0453e-3
0.4	7.9608e-5	9.9901e-4	9.0542e-3	5.0536e-4	3.6589e-3	0.013691
0.5	8.0014e-5	2.0444e-3	0.015368	7.0835e-4	0.023697	0.023674
0.6	8.1086e-5	2.6304e-3	0.036314	7.1143e-4	0.029866	0.056681
0.7	9.3555e-5	7.0569e-3	0.039257	4.6637e-3	0.044125	0.059996
0.8	1.8014e-4	8.7847e-3	0.047361	5.1146e-3	0.063014	0.080114
0.9	3.0965e-4	0.023544	0.054258	5.9335e-3	0.069982	0.088962
1	6.1258e-4	0.050147	0.071236	7.2014e-3	0.080259	0.090481

By using the boundary conditions ψ0⁢(-1)=ψ0⁢(1)=L,ψ0′⁢(0)=0, the solution of the so-called zeroth-order problem (3.2) is

(3.4)ψ0⁢(ξ)=L⁢cosh⁡λ⁢ξcosh⁡λ.

This result is valid only within the range -1≤ψ≤1. Similar solutions can be found in [23], but they may be obtained for a different location of the coordinate system. In practice, researchers usually use this zero-order approximate solution for the sake of simplicity. Nevertheless,one can proceed as follows.

Substituting equation (3.4) into (3.3), we get

(3.5)ψ1′′⁢(ξ)-λ2⁢ψ1⁢(ξ)=-λ2⁢C1⁢L36⁢cosh3⁡λ⁢(cosh⁡λ⁢ξ)3.

The solution of (3.5) with the boundary conditions ψ1⁢(-1)=ψ1⁢(1)=ψ1′⁢(0)=0 is

ψ1(ξ)=-C1⁢L3⁢e-3⁢λ⁢ξ-2⁢λ⁢sec⁡h3⁢λ13824⁢(e2⁢λ+1)(12λe2⁢λ⁢ξ+2⁢λ-12λe2⁢λ⁢ξ+2⁢λ+12λe4⁢λ⁢ξ+2⁢λ-12λe4⁢λ⁢ξ+4⁢λ

-12⁢λ⁢ξ⁢e2⁢λ⁢ξ+2⁢λ-12⁢λ⁢ξ⁢e2⁢λ⁢ξ+4⁢λ+12⁢λ⁢ξ⁢e4⁢λ⁢ξ+2⁢λ+12⁢λ⁢ξ⁢e4⁢λ⁢ξ+4⁢λ-e2⁢λ⁢ξ

-e4⁢λ⁢ξ-e2⁢λ⁢ξ+6⁢λ-e4⁢λ⁢ξ+6⁢λ+e6⁢λ⁢ξ+2⁢λ+e6⁢λ⁢ξ+4⁢λ+e2⁢a+e4⁢a).

With ψ0⁢(ξ) and ψ1⁢(ξ) known, the first-order approximation can be written in the form

ψ(ξ)=ψ0(ξ)+ψ1(ξ)=Lcosh⁡λ⁢ξcosh⁡λ-C1⁢L3⁢e-3⁢λ⁢ξ-2⁢λ⁢sec⁡h3⁢λ13824⁢(e2⁢λ+1)(12λe2⁢λ⁢ξ+2⁢λ-12λe2⁢λ⁢ξ+2⁢λ+12λe4⁢λ⁢ξ+2⁢λ

-12⁢λ⁢e4⁢λ⁢ξ+4⁢λ-12⁢λ⁢ξ⁢e2⁢λ⁢ξ+2⁢λ-12⁢λ⁢ξ⁢e2⁢λ⁢ξ+4⁢λ+12⁢λ⁢ξ⁢e4⁢λ⁢ξ+2⁢λ+12⁢λ⁢ξ⁢e4⁢λ⁢ξ+4⁢λ-e2⁢λ⁢ξ

-e4⁢λ⁢ξ-e2⁢λ⁢ξ+6⁢λ-e4⁢λ⁢ξ+6⁢λ+e6⁢λ⁢ξ+2⁢λ+e6⁢λ⁢ξ+4⁢λ+e2⁢a+e4⁢a).

Following the procedure mentioned in Section 2, the parameter can be easily obtained as

C1=1.20114427441

for λ=1. Table 1 shows the absolute errors of the OHAM and CGA for different values of the inverse Debye length. Note that these errors are obtained only for the first approximate solutions of the OHAM. The CGA values are computed by the algorithm in [6, 12]. All numerical computations were carried out using the Mathematica 9 software.

3.2 Capillary system in cylindrical geometry

In this case, we assume that the unique atom is centered at 0. Then the solution ψ⁢(ξ) to the Poisson–Boltzmann equation can be written as ψ⁢(∥ξ∥) for some function ψ, because of the symmetry of the problem.

The system consists of the cylindrical capillary whose walls are subjected to a dimensionless electrostatic potential equal to L. In dimensionless coordinates, the center of the system is located at ξ=0, while the wall of the capillary is located at ξ=+1 for the positive domain, as seen in Figure 2. Using the Debye–Huckel approximation, equation (1.2) yields the following system of differential equations and boundary conditions:

(3.6)ψ′′+1ξ⁢ψ=λ2⁢(ψ+ψ33!),ψ⁢(+1)=L,ψ′⁢(0)=0.

Figure 2

Geometrical sketch of the capillary channel and the coordinate systemused in Example 2.

Example 2.

Proceeding as in the previous example, the OHAM formulation for problem (3.6) can be constructed as

L⁢(ψ)=ψ′′+ψξ-λ2⁢ψ,g⁢(ξ)=0,N⁢(ψ)=-λ2⁢ψ36,

which satisfies

(1-p)⁢[(ψ0+p⁢ψ1+p2⁢ψ2+⋯)′′+(ψ0+p⁢ψ1+p2⁢ψ2+⋯)ξ-λ2⁢(ψ0+p⁢ψ1+p2⁢ψ2+⋯)]

=(pC1+p2C2+⋯)[(ψ0+pψ1+p2ψ2+⋯)′′+(ψ0+p⁢ψ1+p2⁢ψ2+⋯)ξ

   -λ2(ψ0+pψ1+p2ψ2+⋯)-λ2(ψ0+p⁢ψ1+p2⁢ψ2+⋯)36].

Equating the coefficients yields

(3.7)O⁢(p0): ⁢ψ0′′+ψ0ξ-λ2⁢ψ0=0,

O⁢(p1): ⁢ψ1′′+ψ1ξ-λ2⁢ψ1=-λ2⁢C1⁢ψ036.

By using the boundary conditions ψ⁢(+1)=L, ψ′⁢(0)=0, the solution of (3.7) can be written in the form

(3.8)ψ0⁢(ξ)=L⁢I0⁢(λ⁢ξ)I0⁢(λ),

where I0 is the modified Bessel function of first kind, which is available in the computational software. In the first iteration step, we can get the solution given in [25]. This result is also valid in the range -1≤ψ≤1. One can also proceed to find the first-order solution of (3.6) by using solution (3.8) as in the previous case. Figures 3 and 4 give a comparison between the approximate results of the proposed methods and the Runge method.

Figure 3

Comparison of the CGA solutions (∘) with the Runge–Kutta solutions.

Figure 4

Comparison of the OHAM solutions (∙) with the Runge–Kutta solutions.

4 Results and discussion

The optimal homotopy asymptotic method is applied for the first time to obtain approximate solutions of the electrostatic differential equation which has a simple form, but which is not easy to solve. Therefore, the Debye–Huckel approximation is used to simplify the Poisson–Boltzman equations. Two different geometries are used to illustrate the wide range of validity of this approach. One can easily observe from Example 1 that even the first approximate solution of the OHAM gives more efficient results than the conjugate gradient algorithm. Table 1 displays the absolute errors of these methods. Figure 4 also shows the accuracy of the OHAM solutions when compared to the solutions of the conjugate gradient algorithm in Figure 3.

5 Conclusion

In this paper, we derived more accurate and approximate solutions for the prediction of the electrostatic potential, which are commonly used in electrokinetic research and related applications. The obtained results were also compared with those of other methods in the literature. As a result, more effective solutions were derived for the nonlinear Poisson–Boltzmann equation by means of the Optimal Homotopy Asymptotic Method, and this fact supports the claim that the OHAM is a powerful mathematical tool in dealing with nonlinear equations.