# Nonlinearities Distribution Homotopy Perturbation Method Applied to Solve Nonlinear Problems: Thomas-Fermi Equation as a Case Study

Authors: U. Filobello-Nino; H. Vazquez-Leal; K. Boubaker; A. Sarmiento-Reyes; A. Perez-Sesma; A. Diaz-Sanchez; V. M. Jimenez-Fernandez; J. Cervantes-Perez; J. Sanchez-Orea; J. Huerta-Chua; L. J. Morales-Mendoza; M. Gonzalez-Lee; C. Hernandez-Mejia; F. J. Gonzalez-Martinez
Journal: Journal of Applied Mathematics
 (2015)

---

## Abstract

We propose an approximate solution of T-F equation, obtained by using the nonlinearities distribution homotopy perturbation method (NDHPM). Besides, we show a table of comparison, between this proposed approximate solution and a numerical of T-F, by establishing the accuracy of the results.


---

## Body

## 1. Introduction

The numerical simulation for the charge distribution and the electric field inside an atom is a very difficult task, especially for complex atoms. For these cases, the Thomas-Fermi (T-F) equation can be employed in order to obtain highly accurate approximate results. T-F method is based on the fact that most of complex atoms, with a large number of electrons, have relatively large quantum numbers, and therefore the semiclassical approach can be employed. Then, it is possible to apply the concept of cells in the phase space to the states of individual electrons [1].The Thomas-Fermi equation describes mathematically an infinite atom, without border; however the model is not applicable to both distances from the nucleus: too small and too large. More precisely, its application domain is limited to a range of distancesR: large relative to the quantity 1
/
Z and small compared to 1, where Z is the atomic number of the atom. However, for complex atoms, most of the electrons are in that interval, for which the study of T-F equation becomes important.One of the important consequences of the Thomas-Fermi theory follows from the above. The atom has an outer boundary for values ofR
≅
1; thus the theory predicts that the dimensions of the atom are independent of Z, and for the same reason the ionization energy of the outer electrons is also independent of the atomic number [1].Problems like the one mentioned give rise to the search of solutions to nonlinear differential equations but, unfortunately, solving this kind of equations is a difficult task. As a matter of fact, most of the time, only an approximate solution to such problems can be got. In order to approach various types of nonlinear differential equations, several methods have been proposed such as those based on variational approaches [2–4], tanh method [5], exp-function [6], Adomian’s decomposition method [7–12], parameter expansion [13], homotopy perturbation method [14–27], homotopy analysis method [25, 28–32], and perturbation method [28, 33, 34], among others. Also, a few exact solutions to nonlinear differential equations have been reported occasionally [35].This study proposes a variation of the homotopy perturbation method (HPM), by using nonlinearities distributions [17], which allows finding an approximate solution of Thomas-Fermi equation [23, 36–41]. On the other hand, it will be seen that it is convenient to introduce an adjusting parameter in order to enlarge the domain of convergence of our approach.This paper is organized as follows. In Section2, we provide the basic concept of nonlinearities distributions homotopy perturbation method (NDHPM). Section 3 presents the application of NDHPM to find an approximate solution of Thomas-Fermi equation. Section 4 discusses the main results obtained. Finally, a brief conclusion is given in Section 5.


## 2. Basic Idea of NDHPM

The standard homotopy perturbation method (HPM) was proposed by Ji Huan He. It was introduced like a powerful tool to approach various kinds of nonlinear problems. The homotopy perturbation method can be considered as a combination of the classical perturbation technique and the homotopy (whose origin is in the topology), and it is not restricted to small parameters as traditional perturbation methods. For example, HPM method requires neither small parameter nor linearization but only few iterations to obtain accurate solutions.To figure out how HPM method works, consider a general nonlinear equation in the form(1)
A
u
-
f
r
=
0
,
r
∈
Ω
,
with the following boundary conditions:
(2)
B
u
,
∂
u
∂
n
=
0
,
r
∈
Γ
,
where A is a general differential operator, B is a boundary operator, f
(
r
) a known analytical function, and Γ is the domain boundary for Ω. A can be divided into two operators L and N, where L is linear and N nonlinear; from this last statement, (1) can be rewritten as
(3)
L
(
u
)
+
N
(
u
)
-
f
(
r
)
=
0
.
Generally, a homotopy can be constructed in the form
(4)
H
v
,
p
=
1
-
p
L
v
-
L
u
0
+
p
L
v
+
N
v
-
f
r
=
0
,
llllllllllllllllllll
p
∈
0,1
,
r
∈
Ω
.
Or
(5)
H
v
,
p
=
L
(
v
)
-
L
(
u
0
)
+
p
L
u
0
+
N
v
-
f
r
=
0
,
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
p
∈
0,1
,
r
∈
Ω
,
where p is a homotopy parameter, whose values are within range of 0 and 1; u
0 is the first approximation for the solution of (3) that satisfies the boundary conditions.Assuming that solution for (4) or (5) can be written as a power series of p,
(6)
v
=
v
0
+
v
1
p
+
v
2
p
2
+
⋯
.
Substituting (6) into (5) and equating identical powers of p terms, there can be found values for the sequence v
0
,
v
1
,
v
2
,
….Whenp
→
1, it yields the approximate solution for (1) in the form
(7)
v
=
v
0
+
v
1
+
v
2
+
v
3
+
⋯
.
Another way to build a homotopy, which is relevant for this paper, it is by considering the following general equation:
(8)
L
(
v
)
+
N
(
v
)
=
0
,
where L
(
v
) and N
(
v
) are the linear and nonlinear operators, respectively. It is desired that solution for L
(
v
)
=
0 describes, accurately, the original nonlinear system.By the homotopy technique, a formulation is constructed as follows [16]:
(9)
(
1
-
p
)
L
(
v
)
+
p
[
L
(
v
)
+
N
(
v
)
]
=
0
.
Again, it is assumed that solution for (9) can be written in the form (6); thus, taking the limit when p
→
1 results in the approximate solution of (8). The convergence of solutions obtained by HPM method is discussed in [20, 21, 25, 26].A recent report [17] introduced a modified version of homotopy perturbation method, which eases the solutions searching process for (3). As first step, a homotopy of the following form is introduced:
(10)
H
v
,
p
=
1
-
p
L
v
-
L
u
0
+
p
L
v
+
N
v
,
p
-
f
r
,
p
=
0
.
It can be noticed that the homotopy function (10) is essentially the same as (4), except for the nonlinear operator N and the nonhomogeneous function f, which contain embedded the homotopy parameter p. We propose first to subdivide N and f into a series of terms and multiply p
i by the most nonlinear terms, where i is an integer number greater or equal to zero. The i power is selected according to how much displacement is desired in the interactions for the corresponding nonlinear term of N or f. The introduction of that parameter within the differential equation is a strategy to redistribute the nonlinearities between the successive iterations of the HPM method, and thus it increases the probabilities of finding the sought solution. The rest of the method is exactly the same as the standard procedure for the HPM. Therefore, we establish that
(11)
v
=
v
0
+
v
1
p
+
v
2
p
2
+
⋯
.
When p
→
1, it results in the approximate solution for (3) in the form
(12)
v
=
v
0
+
v
1
+
v
2
+
v
3
+
⋯
.
An advantage of this procedure is that, given the distribution of nonlinearities, from the differential equation, over the successive iterations of (11), less complex analytic approximations may be obtained than those generated by the original standard HPM.Finally, convergence of solutions obtained by NDHPM method is discussed in [17].


## 3. Approximate Solution of Thomas-Fermi Equation

In order to facilitate understanding of the NDHPM method, we will solve approximately the equation(13)
y
′′
=
y
3
/
2
x
1
/
2
,
y
(
0
)
=
1
,
y
′
(
0
)
=
-
L
,
where prime denotes differentiation with respect to x.y
(
x
) is a function related to the electrostatic potential inside the atom, x is a variable proportional to distance from the nucleus (
r
), and the most accurate value of y
′
(
0
) is L
=
1.588071022611375313 [36].The above equation is a modified version of the original problem of Thomas-Fermi, where (13) is subject to the boundary conditions y
(
0
)
=
1, y
(
∞
)
=
0 [28, 41] although there are other conditions that can be applied [40].Following references [16, 24] instead of defining a linear and nonlinear part in the above equation, we add and subtract α
y
′
+
β
y as it is shown, in order to add a term to linear operator to facilitate the search for convergent solutions; thus (13) can be rewritten as
(14)
y
′′
+
α
y
′
+
β
y
-
y
3
/
2
x
1
/
2
+
α
y
′
+
β
y
=
0
,
where α and β are constant parameters to determine.The linear part is identified as(15)
L
y
=
y
′′
+
α
y
′
+
β
y
,
and the nonlinear is
(16)
N
(
y
)
=
-
y
3
/
2
x
1
/
2
+
α
y
′
+
β
y
.
In order to simplify the calculations, we will redistribute the nonlinear term by using NDHPM method, starting from (9), in the form
(17)
y
′′
+
α
y
′
+
β
y
-
p
y
3
/
2
x
1
/
2
+
p
α
y
′
+
p
β
y
=
0
,
by substituting
(18)
y
=
y
0
+
y
1
p
+
y
2
p
2
+
⋯
into (17), and equating terms having identical powers of p we obtain
(19)
y
0
′′
+
α
y
0
′
+
β
y
0
=
0
,
y
0
0
=
1
,
y
0
′
0
=
-
L
,
(20)
y
1
′′
+
α
y
1
′
+
β
y
1
-
y
0
3
/
2
x
1
/
2
=
0
,
y
1
(
0
)
=
0
,
y
1
′
(
0
)
=
0
,
⋮In this paper we study the first order approximation, (19)-(20). We will see that the nonlinear term of N
(
y
) (see (16)) contains the sufficient information to obtain a good analytical approximation for (13).The solution of (19) that satisfies the initial conditions is given by
(21)
y
0
=
B
-
L
B
-
A
e
-
A
x
+
L
-
A
B
-
A
e
-
B
x
,
where A and B are constants related to α and β.By substituting (21) into (20) we obtain
(22)
y
1
′′
+
α
y
1
′
+
β
y
1
-
1
x
1
/
2
B
-
L
B
-
A
e
-
A
x
+
L
-
A
B
-
A
e
-
B
x
3
/
2
=
0
.
In order to simplify (22) and obtain a handy approximation for (13), we will assume that there is an adequate choice of A and B
(
A
≪
B
) such that e
-
A
x
≫
e
-
B
x, valid for 0
<
x
<
∞; in such a way we can rewrite (22) approximately as
(23)
y
1
′′
+
α
y
1
′
+
β
y
1
=
k
1
3
/
2
e
-
3
A
x
/
2
x
1
/
2
1
+
3
k
2
2
k
1
e
-
(
B
-
A
)
x
,
where we have defined
(24)
k
1
=
B
-
L
B
-
A
,
k
2
=
L
-
A
B
-
A
.
This assumption is justified later.To solve (23), we employ the method of variation of parameters [42] which requires evaluating the following integrals:
(25)
u
1
=
-
∫
f
x
e
-
B
x
d
x
W
,
u
2
=
∫
f
x
e
-
A
x
d
x
W
,
where y
H
1
=
e
-
A
x and y
H
2
=
e
-
B
x are the solutions of the homogeneous differential equation
(26)
y
1
′′
+
α
y
1
′
+
β
y
1
=
0
,
W is the Wronskian of these two functions, which is given by
(27)
W
y
1
,
y
2
=
A
-
B
e
-
A
+
B
x
A
≠
B
,
and f
(
x
) is the right-hand side of (23).Substitutingf
(
x
) and (27) into (25) leads to
(28)
u
1
=
-
1
A
-
B
k
1
3
/
2
lim
⁡
ε
→
0
∫
ε
x
e
-
A
y
/
2
d
y
y
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
+
3
2
k
2
k
1
1
/
2
lim
⁡
ε
→
0
∫
ε
x
e
-
(
B
-
A
/
2
)
y
d
y
y
,
u
2
=
1
A
-
B
k
1
3
/
2
lim
⁡
ε
→
0
∫
ε
x
e
-
(
3
A
/
2
-
B
)
y
d
y
y
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
+
3
2
k
2
k
1
1
/
2
lim
⁡
ε
→
0
∫
ε
x
e
-
A
y
/
2
d
y
y
.
By substituting y
=
u
2, the above equations take the form
(29)
u
1
=
-
1
A
-
B
2
k
1
3
/
2
lim
⁡
ε
→
0
∫
ε
x
e
-
u
2
A
/
2
d
u
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
+
3
k
2
k
1
1
/
2
lim
⁡
ε
→
0
∫
ε
x
e
-
(
B
-
A
/
2
)
u
2
d
u
,
(30)
u
2
=
1
A
-
B
2
k
1
3
/
2
lim
⁡
ε
→
0
∫
ε
x
e
(
B
-
3
A
/
2
)
u
2
d
u
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
+
3
k
2
k
1
1
/
2
lim
⁡
ε
→
0
∫
ε
x
e
-
A
u
2
/
2
d
u
.
Therefore, the solution of (23) is written according to the method of variation of parameters as [42]
(31)
y
1
=
f
e
-
A
x
+
g
e
-
B
x
+
u
1
e
-
A
x
+
u
2
e
-
B
x
,
where f and g are constants and u
1 and u
2 are given by (29) and (30), respectively.Applying the initial conditiony
1
(
0
)
=
0 to (31) leads to
(32)
y
1
0
=
f
+
g
-
1
A
-
B
2
k
1
3
/
2
lim
⁡
ε
→
0
x
→
0
∫
ε
x
e
-
A
u
2
/
2
d
u
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
+
3
k
2
k
1
1
/
2
lim
⁡
ε
→
0
x
→
0
∫
ε
x
e
-
(
B
-
A
/
2
)
u
2
d
u
+
1
A
-
B
2
k
1
3
/
2
lim
⁡
ε
→
0
x
→
0
∫
ε
x
e
(
B
-
3
A
/
2
)
u
2
d
u
+
3
k
2
k
1
1
/
2
lim
⁡
ε
→
0
x
→
0
∫
ε
x
e
-
A
u
2
2
d
u
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
+
3
k
2
k
1
1
/
2
lim
⁡
ε
→
0
x
→
0
∫
ε
x
e
-
A
u
2
/
2
d
u
=
0
.
Equation (32) is simplified to
(33)
f
+
g
=
0
.
On the other hand, to apply the condition y
1
′
(
0
)
=
0, we differentiate (31) to obtain
(34)
y
1
′
=
-
A
f
e
-
A
x
-
B
g
e
-
B
x
-
e
-
A
x
A
-
B
k
1
3
/
2
e
-
A
x
/
2
x
+
3
k
2
k
1
1
/
2
e
-
(
B
-
A
/
2
)
x
2
x
+
e
-
B
x
A
-
B
k
1
3
/
2
e
-
3
A
/
2
-
B
x
x
+
3
k
2
k
1
1
/
2
e
-
A
/
2
x
2
x
-
A
e
-
A
x
u
1
-
B
e
-
B
x
u
2
.
After performing algebraic simplifications to (34), we obtain
(35)
y
1
′
=
-
A
f
e
-
A
x
-
B
g
e
-
B
x
-
A
e
-
A
x
u
1
-
B
e
-
B
x
u
2
.
Applying the condition y
1
′
(
0
)
=
0 to (35) the following is obtained:
(36)
A
f
+
B
g
=
0
,
since lim
⁡
x
→
0
u
1
=
0 and lim
⁡
x
→
0
u
2
=
0 (see (29) and (30)).From (33) and (36) we obtain f
=
0 and g
=
0 since A
≠
B; therefore (31) becomes
(37)
y
1
=
u
1
e
-
A
x
+
u
2
e
-
B
x
.
By substituting (21) and (37) into (12), we obtain the following first order approximation for the solution of (13):
(38)
y
=
k
1
e
-
A
x
+
k
2
e
-
B
x
+
u
1
e
-
A
x
+
u
2
e
-
B
x
.
We will show later that it is possible to achieve a good fit and enlarge the domain of convergence for (38) (keeping at the same time the handy character of the proposed solution), by introducing in (38) an adjusting parameter (consistent with the initial conditions), as it is shown in
(39)
y
=
k
1
e
-
A
x
+
k
2
e
-
B
x
+
δ
u
1
e
-
A
x
+
δ
u
2
e
-
B
x
.
The inclusion of parameter, δ, is motivated because the solution y
0
(
x
) for the homogeneous equation (19) is undetermined by a constant factor α (without considering the initial values) and from the following argument.LetY
0
(
x
)
=
α
y
0
(
x
), since Y
0
(
0
)
=
α, y
0
′
(
0
)
=
-
L
α; the first order solution for a function, say Y
(
x
), which defines a slightly different problem from the original, is given by Y
(
x
)
=
Y
0
(
x
)
+
Y
1
(
x
) or Y
(
x
)
=
α
y
0
(
x
)
+
Y
1
(
x
), where y
0
(
0
)
=
1, y
0
′
(
0
)
=
-
L, Y
1
(
0
)
=
0, and Y
1
′
(
0
)
=
0. Thus, after dividing by α, y
(
x
)
=
y
0
(
x
)
+
δ
y
1
(
x
) is obtained (where δ
=
f
(
α
)
/
α
,
y
(
x
)
=
Y
(
x
)
/
α, and we have supposed that for some problems Y
1
(
x
)
=
f
(
α
)
y
1
(
x
). As a matter of fact, following (21)–(30) we conclude that for this case f
(
α
)
=
α
3
/
2; therefore δ
=
α
1
/
2
).Substituting (29) and (30) into (39) we obtain
(40)
y
=
k
1
e
-
A
x
+
k
2
e
-
B
x
-
δ
e
-
A
x
A
-
B
2
k
1
3
/
2
lim
⁡
ε
→
0
∫
ε
x
e
-
A
u
2
/
2
d
u
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
+
3
k
2
k
1
1
/
2
lim
⁡
ε
→
0
∫
ε
x
e
-
(
B
-
A
/
2
)
u
2
d
u
+
δ
e
-
B
x
A
-
B
2
k
1
3
/
2
lim
⁡
ε
→
0
∫
ε
x
e
(
B
-
3
A
/
2
)
u
2
d
u
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
+
3
k
2
k
1
1
/
2
lim
⁡
ε
→
0
∫
ε
x
e
-
A
u
2
/
2
d
u
.
In [14] a high accurate approximation of normal distribution integral was reported, which let us write the above integrals as
(41)
lim
⁡
ε
→
0
∫
ε
x
e
-
A
u
2
/
2
d
u
=
1
2
2
π
A
tanh
39
2
A
x
2
π
-
111
2
tan
-
1
35
111
A
x
2
π
,
lim
⁡
ε
→
0
∫
ε
x
e
-
(
B
-
A
/
2
)
u
2
d
u
=
1
2
2
π
2
B
-
A
tanh
39
2
2
B
-
A
x
2
π
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
-
111
2
tan
-
1
35
111
2
B
-
A
x
2
π
.
Note that
(42)
∫
0
x
e
(
B
-
3
A
/
2
)
u
2
d
u
is not a Gaussian integral on the assumption that e
-
A
x
≫
e
-
B
x for 0
<
x
<
∞, since it implies that A
≪
B.Equations (41) allow us to write (40) as
(43)
y
=
k
1
e
-
A
x
+
k
2
e
-
B
x
-
δ
e
-
A
x
A
-
B
·
35
111
2
B
-
A
x
2
π
k
1
3
/
2
2
π
A
tanh
39
2
A
x
2
π
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
-
111
2
tan
-
1
35
111
A
x
2
π
+
3
2
k
2
k
1
1
/
2
2
π
2
B
-
A
·
tanh
35
111
2
B
-
A
x
2
π
39
2
2
B
-
A
x
2
π
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
-
111
2
tan
-
1
35
111
2
B
-
A
x
2
π
+
δ
e
-
B
x
A
-
B

35
111
A
x
2
π
2
k
1
3
/
2
lim
⁡
ε
→
0
∫
ε
x
e
B
-
3
A
/
2
u
2
d
u
llllllllllllllllll
+
3
2
k
2
k
1
1
/
2
2
π
A
lllllllllllllllllll
·
tanh
39
2
A
x
2
π
-
111
2
tan
-
1
35
111
A
x
2
π


.
To transform (43) into an analytical expression, we evaluate (42) keeping terms up to the ninth power, to obtain the following compact fractional power function:
(44)
∫
0
x
e
(
B
-
3
A
/
2
)
u
2
d
u
=
x
+
1
3
B
-
3
A
x
3
/
2
+
1
10
B
-
3
A
2
x
5
/
2
+
1
42
B
-
3
A
3
x
7
/
2
+
1
216
B
-
3
A
4
x
9
/
2
,
where the main contribution of this approximation to (43) is in the range of x
∈
(
0,1
] because it is multiplied by a negative exponential (see (43) and Section 4).Therefore, substituting (44) into (43), the following is obtained:
(45)
y
=
k
1
e
-
A
x
+
k
2
e
-
B
x
-
δ
e
-
A
x
A
-
B
·
35
111
2
B
-
A
x
2
π
k
1
3
/
2
2
π
A
tanh
39
2
A
x
2
π
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
-
111
2
tan
-
1
35
111
A
x
2
π
+
3
2
k
2
k
1
1
/
2
2
π
2
B
-
A
·
tanh
39
2
2
B
-
A
x
2
π
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
-
111
2
tan
-
1
35
111
(
2
B
-
A
)
x
2
π
+
δ
e
-
B
x
A
-
B

35
111
A
x
2
π
2
k
1
3
/
2
x
+
B
-
3
A
x
3
/
2
3
+
B
-
3
A
2
x
5
/
2
10
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
+
B
-
3
A
3
x
7
/
2
42
+
B
-
3
A
4
x
9
/
2
216
lllllllllllllllllll
+
3
2
k
2
k
1
1
/
2
2
π
A
llllllllllllllllll
·
tanh

39
2
A
x
2
π
-
111
2
tan
-
1
35
111
A
x
2
π


.
In order to obtain a good approximation from (45), we optimize the values of the aforementioned parameters as follows: A
=
0.398521712348839, B
=
6, and δ
=
-
0.623239813 using the procedure reported in [14, 17, 27] (note that with these values, the assumption A
≪
B is satisfied, besides the above results allow to know α and β whose values are not needed in the search for the solution of the system (19) and (20)).Figure1 and Table 1 show the comparison between approximate solution (45) and the exact solution. The accuracy of (45) is clear as an approximate solution for (13).Table 1
Numerical comparison of Thomas-Fermi’s numerical solution [39] and approximations [15, 38] and (45).
x
Numerical [39]
[38]
(45)
[15]
0.00
—
1.000000000
1.000000000
1.000000000
0.25
0.755880759
0.776191000
0.708502932
0.680650028
0.50
0.606700008
0.615917000
0.570491745
0.459455663
0.75
0.502964042
0.505380000
0.485018011
0.307042537
1.00
0.424333179
0.423772000
0.421167227
0.202655526
1.25
0.363227937
0.362935000
0.369542163
0.131667603
1.50
0.314660642
0.314490000
0.326298526
0.083800106
1.75
0.275233848
0.275154000
0.289315779
0.051853222
2.00
0.242678587
0.242718000
0.257253699
0.030801947
2.25
0.215439334
0.215630000
0.229210094
0.017153465
2.50
0.192406328
0.192795000
0.204540381
0.008491277
2.75
0.172758691
0.173364000
0.182755264
0.003152503
3.00
0.155871862
0.156719000
0.163464474
0
3.25
0.141260504
0.142371000
0.146346052
-
0.001738122
3.50
0.128541381
0.129937000
0.131128799
-
0.002581206
3.75
0.117408054
0.119108000
0.117581337
-
0.002874924
4.00
0.107612958
0.109632000
0.105504623
-
0.00284628
4.25
0.098954329
0.101303000
0.094726427
-
0.002641801
4.50
0.091266456
0.093950400
0.085097045
-
0.002353931
4.75
0.084412289
0.087432000
0.076485875
-
0.002039168
5.00
0.078277758
0.081629600
0.068778618
-
0.001730443Figure 1
Comparison of Thomas-Fermi’s numerical solution [39] and approximation (45).


## 4. Discussion

Table1 shows the comparison between numerical solution given in [39] and approximations [15, 38] and (45). It is clear that our approximation (45) is competitive, considering that it is analytic and handy, unlike other approximations which are expressed in terms of power series. In [37, 38, 41], T-F equation was solved by homotopy analysis method (HAM). The approximation [38] corresponds to an 80th order analytical approximation and therefore is impractical. On the other hand in [37] the results are given in terms of infinite series and show particularly the calculation of the initial slope y
′
(
0
). Moreover, in [41], an accurate approximate solution was obtained, but it corresponds to the 40th order approximation. In [23], (13) was solved by HPM and Padé approximants; in this case a solution of 26 power series terms is provided which is also not handy in practice. Reference [15] presents an analytical approach, but it is not accurate. In this study, we used homotopy perturbation method with nonlinearities distributions to find an accurate solution for (13) (see Figure 1 and Table 1). Although (43) is just almost analytic (due to the presence of a non-Gaussian integral), it can be estimated that contribution of (42) to the solution (43) is well represented by (44) assuming the values are A
=
0.398521712348839, B
=
6, and δ
=
-
0.623239813; thus we obtained the analytical approximation (45). We just kept only five terms in (44) in order to obtain a handy approximation. Besides, the right-hand side of (44) adequately represents the indicated non-Gaussian integral in the range x
∈
[
0,1
]. It implies that the best contribution of the term proportional to e
-
B
x
∫
0
x
e
(
B
-
3
A
/
2
)
u
2
d
u in (43) is just beginning of the interval considered, that is, in a part of the domain where the T-F equation describes correctly the atom (see Section 1). For values x
>
1, the above-mentioned term is dominated for the exponential e
-
B
x. In any case, the value of the integrals of (40) can be evaluated by using numerical methods. An interesting fact is that if more terms had been considered in the development of (12), we would have obtained, as a solution of (13), a series containing terms of hyperbolic tangent functions to the form (41).The strategy to redistribute the nonlinearities of T-F equation, between the successive iterations of the HPM method, was important to obtain an approximate solution of (13) because it allows us to simplify the solution procedure for the successive iterations of the HPM method. Besides, we obtained a soluble coupled system of linear differential equations (19) and (20). In this study, we introduced an adjusting parameter in order to enlarge the domain of convergence of our approach, with good results, by using only the first order approximation (see Table 1 and Figure 1).A relevant fact is that NDHPM requires little iterations to obtain accurate solutions, if we have a first approximation containing as much information as possible for a nonlinear differential equation. For instance, in our case, (21) contains the correct asymptotic character of the exact equation.


## 5. Conclusion

In this work, NDHPM was used to obtain an approximate analytical solution for T-F equation. Our solution is novel because it is expressed in terms of exponential and fractional power functions and above all in terms of normal distribution integrals. It is important to obtain an analytical expression that provides a good description of the solution for the nonlinear differential equations like (13). For instance, the charge distribution and the electrostatic potential inside an atom are adequately described by (43) and (45). It is worth mentioning that, if the initial guess is suitably chosen, it is possible to obtain by this method an accurate approximation like (43), even using the first terms of (12). Finally, in contrast to Runge-Kutta numerical solution, NDHPM method allows to analyze quantitatively and qualitatively the solution.



---

*Source: 405108-2015-01-20.xml*