Filomat 30:6 (2016), 1633–1650
DOI 10.2298/FIL1606633T
Published by Faculty of Sciences and Mathematics,
University of Ni
ˇ
s, Serbia
Available at: http://www.pmf.ni.ac.rs/filomat
An Effective Approach for Evaluation of the Optimal Convergence
Control Parameter in the Homotopy Analysis Method
Mustafa Turkyilmazoglu
a
a
Department of Mathematics, Hacettepe University, 06532-Beytepe, Ankara, Turkey
Abstract. A rapid and effective way of working out the optimum parameter of convergence control in
the homotopy analysis method (HAM) is introduced in this paper. As compared with the already known
ways of evaluating the convergence control parameter in HAM either through the classical h−level curves
with h being the convergence control parameter or from the classical squared residual formula as adopted
in the HAM society, an elegant way of calculating the convergence control parameter yielding the same
optimum values is offered. In most cases, the new method is shown to perform quicker and better against
the residual error method when integrations are much harder to evaluate or even by numerical means.
Examples originating from real life applications selected from the literature demonstrate the validity and
usefulness of the introduced technique.
1. Introduction
Ever since the development of the HAM by Liao in 1992 [4], more than a couple of thousand nonlinear
mathematical models have been treated by researchers using HAM. A good collection of such problems,
the power of the technique and its relation to already known approximate analytic methods were contained
in the books [5] and [8].
It is now well-understood that the success of the HAM is constantly attributed to the so-called convergence
control parameter, whose optimum value can be worked out making use of the traditional squared residual
error definition after the work of [7], see also [3, 12, 13]. Desirable progress was also achieved on the
convergence of HAM [9].
The present paper proposes a new way of finding out the optimum value of convergence control
parameter used to ensure the convergence of the HAM series in a fastest manner, an idea first introduced in
[9] (see chapter 5 in [9]). The proposed approach constitutes an alternative to both the classical h-level curves
method and the squared residual error approach for determination of optimal value of the convergence
control parameter. It is demonstrated through examples from the open literature that the squared residual
error method and the present one generate nearly identical optimal convergence control parameter values,
even though less computational cost is the advantage of the presented scheme. Not only the optimum value
of the convergence control parameter is obtained via the new method, but also the interval of convergence
can be gained and the guarantee of quick convergence can be given.
2010 Mathematics Subject Classification. 34Axx; 35Axx; 37Axx; 39Axx; 41Axx; 45xx; 76Axx
Keywords. Homotopy analysis method, Convergence control parameter, h-curve, Squared residual approach, Ratio approach
Received: 11 May 2014; Accepted: 13 December 2014
Communicated by Dragan S. Djordjevi
´
c
Email address: turkyilm@hacettepe.edu.tr (Mustafa Turkyilmazoglu)
M. Turkyilmazoglu / Filomat 30:6 (2016), 1633–1650 1634
2. Traditional Ways and the Presently Proposed Approach
The layout of the HAM is now well-documented in the literature and hence will be omitted here for the
conciseness, anyway, the interested readers may refer to the very recent book [9]. To be more precise, it is
intended to construct a homotopy of the form
(1 − p)L[u(t, p) − u
0
(t)] − p h H(t)N[u(t, p)] = 0, (1)
to a nonlinear problem N(u, t) = 0, where p ∈ [0, 1] is the homotopy embedding parameter, h is the convergence
control parameter adjusting the convergence, Lis a selected linear operator, u
0
(t) is an initial guess for wanted
solution u(t), H(t) is an auxiliary function. Upon adequate number of successive differentiations of (1), the
following homotopy series is reached
u(t) = u
0
(t) +
∞
X
k=1
u
k
(t), (2)
which was proved by Liao [5] to correspond to the exact solution desired.
By truncating the homotopy series (2), eventually an approximate analytic solution of Mth-order for a
physical problem of interest can be represented by
u
M
(t) = u
0
(t) +
M
X
k=1
u
k
(t), (3)
whose limiting value
lim
M→∞
u
M
(t)
leads to the the exact solution u(t).
Based on a sufficient theorem to ensure the convergence of homotopy series (2), a novel approach was
outlined [9] (see chapter 5 in [9]), which is restated in the following Corollary
Corollary 2.1. For a preassigned value of h, for convergence of the homotopy series (2) it is sufficient to keep track of
magnitudes of the ratio β defined by
β =
ku
k+1
(t)k
ku
k
(t)k
, (4)
and to check whether it remains less than unity for increasing values of k. An optimal value for the convergence
control parameter h could also be determined from (4) by requiring the ratio β to be as close to zero as possible, so that
for such a value the rate of convergence of homotopy series (2) will be the fastest, since then the remainder of the series
will most rapidly decay.
Moreover, in the case of a preassigned value h of the convergence control parameter, when the absolute
value norm is chosen, the region of t for the wanted solution can also be identified. Additionally, the ratio
given in (4) can also plot the constant h-curves [5] correctly. Furthermore, on the condition that one is
concerned with the zeros of equation f (x) = 0 of the algebraic kind, the inequality
|x
k+1
|
|x
k
|
< 1 (5)
for large k produces the convergence control parameter interval.
In place of the constant h−level curves, the norm defining the residual
Res(h) = kN[
M
X
k=0
u
k
(t)]k,
M. Turkyilmazoglu / Filomat 30:6 (2016), 1633–1650 1635
where the norm is in the sense L
p
, with p = 2 in general [16], the subsequent squared residual error is often
employed to determine an optimal h
Res(h) =
Z
A
N
M
X
k=0
u
k
(r)
2
dr, (6)
for which the physical problem takes place over the set A. In the case of a positive integrand, the residual
error in (6) can be replaced by
Res(h) =
Z
A
N
M
X
k=0
u
k
(r)
dr, (7)
based on L
1
norm to gain more computational time. We should mention here that following [7], the discrete
forms corresponding to (6) and (7) may be suggested
Res(h) ≈
1
N + 1
N
X
j=0
N
M
X
k=0
u
k
(t
j
)
2
, (8)
Res(h) ≈
1
N + 1
N
X
j=0
N
M
X
k=0
u
k
(t
j
)
. (9)
with equally distributed N discrete points. By requiring from (6), (7), 8 or (9) that
dRes(h)
dh
= 0, (10)
optimal values of the convergence control parameter h for a nonlinear physical problem can be determined.
The main disadvantages of the residual are that exact integration from (6) and (7) may not be clear or
computational cost of discrete forms (8) or (9) may not be at the expense of desire, particularly for the
unbounded physical problems. Owing to such shortcomings, by means of the ratio given in equation (4), a
novel and easy way of finding h was suggested in [9] by letting the ratio β in (4) to be sufficiently small. If
possible then the optimums of the convergence control parameter h might be as a consequence of
dβ
dh
= 0,
or at worst, graphs of constant β−curves will play the role of determining an optimal value of h in (4).
Therefore, in line with [9],
β =
R
A
u
p
k+1
(r)dr
R
A
u
p
k
(r)dr
, (11)
or its discrete counterpart
β ≈
P
N
j=0
[u
k+1
(t
j
)]
p
P
N
j=0
[u
k
(t
j
)]
p
, (12)
will serve good to get the optimum convergence control parameter h.
It is noted that whenever exact solution or numerical estimation u
e
(t) is at our disposal for a considered
problem, the absolute error
err =
Z
A
|u
e
(t) − u(t)|dt (13)
may be employed to check the accuracy of homotopy solution u(t) as obtained from (3).
M. Turkyilmazoglu / Filomat 30:6 (2016), 1633–1650 1636
M=6, 10, 16, 22
-1.0
-0.8
-0.6
-0.4
-0.2
0.0
0.0
0.2
0.4
0.6
0.8
h
u
Figure 1: h−level curves for equation (14).
3. Examples
In order to test the new approach and its validity through the utilities (6–12) (refer also to chapter 5
in [9] for further analysis), some physical problems are considered here collected from various homotopy
analysis research in the open literature.
3.1. A transcendental equation
Our first example is the transcendental equation
f (u) = ue
u
− 1 = 0, (14)
whose numerical solution up to ten significant digits is simply u = 0.5671432904. The choices
u
0
= 0, L(u) = f (u) − f (u
0
)
helps us to construct the homotopy series solution via the homotopy approach (1) from the residual and
absolute errors
Res(h) = ue
u
− 1, (15)
err = u − ProductLog[1]. (16)
It is computed from (15) that the 22th-order approximation h = −0.44 yields the minimum value for the
residual.
Figure 1 clearly demonstrate that the interval of convergence for h is h ∈ [−0.6, 0). Indeed, from (5) by
analytically solving the inequality |u
22
/u
21
| < 1, we find −0.55721 < h < 0. At the above calculated value of
h = −0.44, Table 1 gives the root and the absolute error (15) and also the result using the classical algorithm
of modified Newton iteration (see [8])
u
k
= u
k−1
+ h
f (u
k−1
)
f
0
(u
k−1
)
. (17)
From Table 1 it is understood that h = −0.44 is good enough for the homotopy analysis method,
as verified also from Table 2, even surprisingly better convergent HAM solutions are obtained than the
Newton method. Tables 2-3 and Figures 2 (a–b) are clear evidences that as M tends to infinity, the optimums
computed from the residual using (15) and the ratio with (5) will be equal.
The ratio (5) can also be assessed from Figure 3 and Table 4 corresponding to the optimum convergence
control parameter h = −0.423. The insurance of convergence is due to the less than unity value of β which
actually limits to 0.57634. It is also observed that to find the minimum values of (15) for M = 151, 201,
251 and 351, respectively 40, 70, 110 and 926 seconds are needed, whereas only 5, 14, 36 and 82 seconds
are sufficient for evaluating minimums from (5). This obviously points to advantageous CPU time for the
present approach.
M. Turkyilmazoglu / Filomat 30:6 (2016), 1633–1650 1637
M 5 10 15 20
u
a
0.5612802859 0.5668842143 0.5671363654 0.5671427617
err
a
5.8630 × 10
−3
2.5908 × 10
−4
6.9250 × 10
−6
5.2870 × 10
−7
u
b
0.5566418620 0.5665732963 0.5671119237 0.5671415630
err
b
1.0501 × 10
−2
5.6999 × 10
−4
3.1367 × 10
−5
1.7274 × 10
−6
Table 1: Zeros of (14) and absolute errors for some M at h = −0.44.
a
Solutions from homotopy (3) and
b
Solutions from Newton
iteration (17).
M 10 20 30 50 100 200 300 350
h -0.4501 -0.4392 -0.4350 -0.4312 -0.4280 -0.4262 -0.4255 0.4253
Table 2: The resulting values of optimum for h using the ratio (5) and residual (15).
M 21 51 101 151 201 251 301 351
h
a
−0.4455 −0.4341 −0.4295 −0.4278 −0.4270 −0.4264 −0.4260 −0.4257
h
b
−0.4074 −0.4162 −0.4196 −0.4209 −0.4216 −0.4220 −0.4223 −0.4225
β 0.57160 0.57523 0.57604 0.57621 0.57626 0.57628 0.57629 0.57629
Table 3: The values of optimal h and ratio β.
a
Equation (15) and
b
Equation (5).
M 100 150 200 300 400 460 490 500
β 0.59136 0.58438 0.58082 0.57771 0.57667 0.57643 0.57636 0.57634
Table 4: The values of ratio β with M = 500 and h = −0.423.
3.2. A nonlinear differential-difference equation of Volterra type
The next example is the famous nonlinear Volterra differential-difference initial value problem
u
0
n
(t) = u
n
(t)(u
n−1
(t) − u
n+1
(t) + u
n−2
(t) − u
n+2
(t)), u
n
(0) = n, 0 ≤ t ≤ 1, (18)
whose physical importance was highlighted in [8]. The exact solution was given in the study [18] as
u
n
(t) =
n
1+6t
.
By means of selecting, respectively,
H(t) = 1, L =
d
dt
, u
n,0
(t) = n − t,
the inequality outlined in (5) for u
0
10
(0) gives rise to the the interval [−2,0] h, which is also depicted from the
h-level curves in Figure 4.
Carrying out the numerical integration for the residual (with 500 equal points) and exact integration for
the ratio at n = 10, it is seen from Figures 5 (a–b) and Table 5 that the residual (8) and the ratio (11) (with
p = 2) generate nearly the same region of convergence, both limiting towards the optimal value h = −0.245.
The fast convergence via the present ratio approach is also apparent, see also plot of the ratio β in Figures 6
(a–c).
