- Email: [email protected]
On the local convergence of a family of two-step iterative methods for solving nonlinear equations
Journal of Computational and Applied Mathematics 255 (2014) 753–764
Contents lists available at ScienceDirect
Journal of Computational and Applied Mathematics journal homepage: www.elsevier.com/locate/cam
On the local convergence of a family of two-step iterative methods for solving nonlinear equations Miquel Grau-Sánchez ∗ , Miquel Noguera, José L. Diaz-Barrero Technical University of Catalonia, Department of Applied Mathematics II and III, Jordi Girona 1-3, Omega, 08034 Barcelona, Spain
article
info
Article history: Received 27 February 2013 Received in revised form 14 May 2013 MSC: 41A25 65H10 Keywords: Order of convergence Nonlinear equations Iterative methods Efficiency
abstract A local convergence analysis for a generalized family of two step Secant-like methods with frozen operator for solving nonlinear equations is presented. Unifying earlier methods such as Secant’s, Newton, Chebyshev-like, Steffensen and other new variants the family of iterative schemes is built up, where a profound and clear study of the computational efficiency is also carried out. Numerical examples and an application using multiple precision and a stopping criterion are implemented without using any known root. Finally, a study comparing the order, efficiency and elapsed time of the methods suggested supports the theoretical results claimed. © 2013 Elsevier B.V. All rights reserved.
1. Introduction There are a great variety of iterative methods for solving a system of nonlinear equations F (x) = 0, where F : D ⊆ Rm → R , and D is a non-empty open convex subset of Rm that contains a simple root α of F . A classical iterative process for solving nonlinear equations is Chebyshev’s method (see [1–3]) m
x ∈ D, 0 yn = xn − F ′ (xn )−1 F (xn ) xn+1 = yn − 1 F ′ (xn )−1 F ′′ (xn ) (yn − xn )2 , 2
n ≥ 0.
The above one-point iterative scheme depends explicitly on the two first derivatives of F . In [1], Ezquerro and Hernández present some modifications in Chebyshev’s method by reducing in one the number of evaluations of the first derivative and maintaining third-order of convergence. It has the following form:
x ∈ D, 0 zn = xn − a F ′ (xn )−1 F (xn ), xn+1 = xn − 1 F ′ (xn )−1 (a2 + a − 1) F (xn ) + F (zn ) , a2
n ≥ 0.
Using the well-known Secant method [4], in [5] a generalization of it employing the divided difference operator of order one (namely, Bn = [xn−1 , xn ; F ]) that substitutes the derivative of F (F ′ (xn ) ≡ [xn , xn ; F ]) is given. The authors call this family
∗
Corresponding author. Tel.: +34 93 413 79 82. E-mail address: [email protected] (M. Grau-Sánchez).
0377-0427/$ – see front matter © 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.cam.2013.06.043
754
M. Grau-Sánchez et al. / Journal of Computational and Applied Mathematics 255 (2014) 753–764
the Chebyshev–Secant-type method and it is defined by
x−1 , x0 ∈ D, 1 zn = xn − a B− n F (xn ), −1 xn+1 = xn − Bn (b F (xn ) + c F (zn )) ,
n ≥ 0,
where a, b, c are non-negative parameters to be chosen so that the sequence {xn } converges to α with maximum local order of convergence. The work presented in [6] analyzes free-derivative iterative processes considering the operator Cn = [xn , L(xn ); F ] and they are called the Steffensen-type method:
x0 ∈ D , zn = xn − a Cn−1 F (xn ), xn+1 = xn − Cn−1 (b F (xn ) + c F (zn )) , m
n ≥ 0,
m
where L : D ⊆ R → R . Our goal in the present paper is to unify these methods by considering the following iterative family
x −1 , x 0 ∈ D , zn = xn − a Θn−1 F (xn ), xn+1 = xn − Θn−1 (b F (xn ) + c F (zn )) ,
(1) n ≥ 0,
where Θn = [G(xn−1 , xn ), H (xn−1 , xn ); F ] and G, H : D ⊆ Rm × Rm → Rm . As can be observed, the above family includes as particular cases the preceding works [5,6] and other well-known algorithms. Furthermore, (1) also offers us the possibility to suggest new methods as we will see later on. 2. Local order of convergence To obtain explicitly an expression of the inverse operator of Θn , we denote by J indistinctly both G or H. On account that a necessary condition about functions G and H is: G(α, α) = H (α, α) = α , then from Taylor’s formulae we get the error expression
εJ = J (xn−1 , xn ) − α = J1 en−1 + J2 en + 1/2 J11 e2n−1 + J22 e2n + 2J12 en−1 en + · · · , ∂J
∂J
∂ 2J
where J1 = ∂ x (α, α), J2 = ∂ x (α, α), J11 = 2 (α, α), etc. n ∂ xn−1 n−1 The expression of Θn is given in the previous works appeared in [7,8]. Namely,
Θn = [G(xn−1 , xn ), H (xn−1 , xn ); F ] = Γ I + A2 (εG + εH ) + A3 (εG2 + εG εH + εH2 ) + A4 (εG3 + εG2 εH + εG εH2 + εH3 ) + · · · , where Γ = F ′ (α) ∈ L(Rm ), Ak = k1! Γ −1 F (k) (α), Ak ∈ Lk (Rm , Rm ) and Ak is k-symmetric. From the preceding follows
Θn−1 = I − A2 (εG + εH ) − (A3 − A22 ) (εG2 + εG εH + εH2 ) + A22 εG εH + A2 (εG + εH ) · A3 (εG2 + εG εH + εH2 ) + A3 (εG2 + εG εH + εH2 ) · A2 (εG + εH ) − A4 (εG3 + εG2 εH + εG εH2 + εH3 ) Γ −1 , where A22 εG2 = (A2 εG )2 . Setting En = zn − α, en = xn − α and using Taylor’s development yield F (xn ) = Γ en + A2 e2n + o(e2n ) .
Subtracting α from both sides of the first equation of (1), we obtain En = en − a Θn−1 Γ en + A2 e2n + o(e2n )
= (1 − a) en + a A2 (G1 + H1 ) en−1 en 1 2 2 2 +a A2 (G11 + H11 ) + A3 (G1 + H1 + G1 H1 ) − (A2 (G1 + H1 )) e2n−1 en 2
+ a A2 (G2 + H2 − I ) e2n + o(e2n−1 en , e2n ),
(2)
where the following notation was used: A2 (G1 + H1 ) en−1 en = A2 (G1 en−1 + H1 en−1 ) en , A2 (G11 + H11 ) e2n−1 en = A2 (G11 e2n−1 + H11 e2n−1 ) en , A3 (G21 + H12 + G1 H1 ) e2n−1 en = A3 (G1 en−1 )(G1 en−1 ) + (H1 en−1 )(H1 en−1 ) − (G1 en−1 )(H1 en−1 ) en ,
M. Grau-Sánchez et al. / Journal of Computational and Applied Mathematics 255 (2014) 753–764
755
and
(A2 (G1 + H1 ))2 e2n−1 en = (A2 (G1 en−1 + H1 en−1 )) (A2 (G1 en−1 + H1 en−1 )) en . Note that with this notation we may order the terms of the error vector En according to their order of convergence. Recall γ β that if we have the term K en−1 en , its local order of convergence is the unique positive root of the associated polynomial 2 equation t − β t − γ = 0. Then (2) is the expression of the error En up to second order. From Taylor’s formulae F (zn ) = Γ (En + o(En )) and subtracting α from both sides of the second equation of (1), we have en+1 = en − Θn−1 b Γ
en + A2 e2n + o(e2n ) + c Γ (En + o(En ))
= (1 − b − c + ac ) en + (b + c − 2ac ) A2 (G1 + H1 ) en−1 en + o(en−1 en ).
(3)
Annulling the coefficients of the terms of low order of (3) yields
1 − b − c + ac = 0, b + c − 2ac = 0,
and we get b = 2 −
1 a
and c =
1 . a
Substituting these values in (3), we obtain
en+1 = (A2 (G1 + H1 ))2 e2n−1 en + (1 − a) A2 e2n + o(e2n−1 en , e2n ).
(4)
Hereafter, a = 1 and G1 + H1 ≡ 0 will be used in order that the family (1) becomes a super-quadratic iterative method. Namely, from (1) and taking into account the previous considerations, we have
x−1 , x0 ∈ D, zn = xn − Θn−1 F (xn ), xn+1 = zn − Θn−1 F (zn ),
(5) n ≥ 0.
Note that for a = b = c = 1, we obtain an iterative method of two steps very similar ones with the same operator. This fact is the reason why this kind of methods are called frozen iterative methods, since the operator is repeated. If we consider (5), where G1 + H1 ≡ 0, then (2) becomes
En =
1 2
A2 (G11 + H11 ) + A2 G21
+
1
A2 (G111 + H111 ) +
6
e2n−1 en + A2 (G2 + H2 − I ) e2n
1
A3 G1 (G11 − H11 ) e3n−1 en
2
+ (A2 (G12 + H12 ) + A3 G1 (G2 − H2 )) en−1 e2n + o(en−1 e2n ), where A2 (G2 + H2 − I ) e2n = A2 (G2 en + H2 en − en ) en , A2 (G111 + H111 ) e3n−1 en = A2 (G111 e3n−1 + H111 e3n−1 ) en , A3 G1 (G11 − H11 ) e3n−1 en = A3 (G1 en−1 ) G11 e2n−1 − H11 e2n−1 en ,
A2 (G12 + H12 ) en−1 e2n = A2 (G12 en−1 en + H12 en−1 en ) en , and A3 G1 (G2 − H2 ) en−1 e2n = A3 (G1 en−1 )(G2 en − H2 en ) en . Furthermore, (4) becomes
e n +1 =
1 4
(A2 (G11 + H11 )) + A2 A3 (G11 + H11 ) 2
+ A2
1 2
G21
A2 (2(G2 + H2 ) − I ) (G11 + H11 ) +
+
2 A3 G21
A3 G21
e4n−1 en
(2(G2 + H2 ) − I ) e2n−1 e2n + o(e2n−1 e2n ),
where A2 A3 (G11 + H11 ) G21 e4n−1 en =
1 2
A2 (G11 e2n−1 + H11 e2n−1 )A3 (G1 en−1 )2
+
1 2
A3 G21
2
e4n−1 en = A3 (G1 en−1 )2
en
A3 (G1 en−1 )2 A2 (G11 e2n−1 + H11 e2n−1 ) en ,
and
A3 (G1 en−1 )2 en .
(6)
756
M. Grau-Sánchez et al. / Journal of Computational and Applied Mathematics 255 (2014) 753–764
√
The orders of the first two terms in (6) are 1 +
17 /2 = 2.562 and 1 +
√
3 = 2.732 respectively. But if we set G1 ≡ 0
and G11 + H11 ≡ 0 the two terms above described and also the terms of third order in e6n−1 en and e3n−1 e2n are zero. Therefore, the only nonzero term of third order is en+1 = A2 (G2 + H2 ) A2 (G2 + H2 − I ) e3n + o(e3n ).
(7)
Now we are ready to state the following result. Theorem 2.1. If a = b = c = 1, G1 ≡ H1 ≡ 0 and G11 + H11 ≡ 0, then iterative method (1) has R-order of convergence at least 3. Furthermore, the vectorial equation of the error is given in (7). Finally, we conclude this section presenting five particular cases of the preceding results. Case 1. From (4), if b = 2 − 1/ a and c = 1/ a, family (1) is a set of iterative methods with memory of at least second order. An element of this family can be found in [7] where the scheme is called the frozen secant method. Here, it is set as Θ = [G, H ] = [xn−1 , xn ; F ] and it is obtained as
x−1 , x0 ∈ D, zn = xn − [xn−1 , xn ; F ]−1 F (xn ), Φ1 : xn+1 = zn − [xn−1 , xn ; F ]−1 F (zn ),
(8) n ≥ 0,
where a = 1, but G1 + H1 ≡ 1 ̸= 0. Others works with this operator can be found in [9–11]. Case 2. Taking G(xn−1 , xn ) = γ xn−1 + (1 − γ )xn
and H (xn−1 , xn ) = δ xn−1 + (1 − δ)xn ,
(9)
if we impose that G1 + H1 ≡ 0, then we obtain γ +δ = 0, and from (6) we have a one-parametric family with super-quadratic order after adding the condition a = b = c = 1. For γ = 1 (δ = −1), the operator Θ becomes [xn−1 , 2xn − xn−1 ; F ] that is the operator used by Kurchatov [12]. This algorithm, called the frozen Kurchatov method, can be found in [13] and it is defined by
x−1 , x0 ∈ D, zn = xn − [xn−1 , 2xn − xn−1 ; F ]−1 F (xn ), Φ2 : xn+1 = zn − [xn−1 , 2xn − xn−1 ; F ]−1 F (zn ),
(10) n ≥ 0,
where a = b = c = 1, G1 + H1 ≡ 0, G11 + H ≡ 0. Notice that in general, in (9), G1 = γ ̸= 0. Hence, the local √ 11 convergence order of this method is 1 + 17 / 2 ≈ 2.56155. Observe that if (5) is an iterative method without memory, then G1 ≡ H1 ≡ 0. In other words, the iterative methods without memory will be super-quadratic if a = b = c = 1. Case 3. Setting γ = 0 (δ = 0) in (9), the operator considered in Case 2 becomes [xn , xn ; F ] ≡ F ′ (xn ) and all the hypotheses of Theorem 2.1 hold, because now G1 ≡ 0. This scheme is known as the frozen Newton method [1,14]:
x0 ∈ D , zn = xn − [F ′ (xn )]−1 F (xn ), Φ3 : xn+1 = zn − [F ′ (xn )]−1 F (zn ),
(11) n ≥ 0.
As we have set G = H = xn , we get G2 = H2 = I, and substituting in (7) we obtain the well known error equation of the frozen Newton method: en+1 = 2 A22 e3n + o(e3n ). Case 4. A family of iterative methods without memory was presented by many authors (see [15–17] and the references therein). Considering G = xn + λF (xn ) and H = xn + ν F (xn ), where λ, ν ∈ R, and on account that the hypotheses of the theorem hold, yields
x0 ∈ D, zn = xn − [xn + λF (xn ), xn + ν F (xn ); F ]−1 F (xn ), Φ4 : xn+1 = zn − [xn + λF (xn ), xn + ν F (xn ); F ]−1 F (zn ),
(12) n ≥ 0.
In this case, we have G2 = I + λ Γ and H2 = I + ν Γ . The error equation (7) becomes en+1 = A2 2 I + (λ + ν) Γ A2 I +
(λ + ν) Γ e3n + o(e3n ). These methods are called Steffensen-like methods.
Case 5. The last family considered in this paper is a set of algorithms that we call quadratic. The definition of a component (i) (i) (i) (i) (i) (i) function of G is given by Gi (xn−1 , xn ) = xn + r (xn − xn−1 )2 , likewise, H i (xn−1 , xn ) = xn + s (xn − xn−1 )2 , where x =
(x(1) , . . . , x(m) ) ∈ Rm , r , s ∈ R and 1 ≤ i ≤ m. Putting qn = (q(n1) , . . . , q(nm) ) ∈ Rm where q(ni) = (x(ni) − x(ni−) 1 )2 , the resulting new family of iterative methods is
x−1 , x0 ∈ D, zn = xn − [xn + r qn , xn + s qn ; F ]−1 F (xn ), Φ5 : xn+1 = zn − [xn + r qn , xn + s qn ; F ]−1 F (zn ),
(13) n ≥ 0.
M. Grau-Sánchez et al. / Journal of Computational and Applied Mathematics 255 (2014) 753–764
757
Table 1 Local convergence order ρi of the methods Φi for 1 ≤ i ≤ 5.
ρ1
ρ2
2
1+ 17 2
ρ3
ρ4
ρ5
3
3
3
√
Table 2 Cost of the methods Φi for i = 1, 2, 4, 5. Using [y, x; F ](1) defined in (15)
Φ1 Φ2 Φ4 Φ5
m 6 m 6 m 6 m 6
Using [y, x; F ](2) defined in (16)
(2m2 + 6mµ+ 9m + 9ℓm + 6µ+ 9ℓ− 11) (2m2 + 6mµ+ 9m + 9ℓm + 12µ+ 9ℓ− 11) (2m2 + 6mµ+ 9m + 9ℓm + 12µ+ 9ℓ− 11) (2m2 + 6mµ+ 9m + 9ℓm + 18µ+ 9ℓ− 5)
m 6 m 6 m 6 m 6
(2m2 + 12mµ + 9m + 9ℓm + 9ℓ − 11) (2m2 + 12mµ+ 9m + 9ℓm + 6µ+ 9ℓ− 11) (2m2 + 12mµ+ 9m + 9ℓm + 6µ+ 9ℓ− 11) (2m2 + 12mµ+ 9m + 9ℓm + 12µ+ 9ℓ− 5)
Taking into account that G1 ≡ H1 ≡ 0 and imposing that G11 + H11 ≡ 0, we have r + s = 0. Furthermore, from G2 + H2 = 2 I and (7) we obtain en+1 = 2 A22 e3n + o(e3n ). Finally, in Table 1 the local convergence orders of the iterative families Φi , 1 ≤ i ≤ 5 are shown. 3. Computational efficiency index The computational efficiency index (CEI ) and the computational cost per iteration (C ) are defined by 1
CEI (µ, µ1 , m, ℓ) = ρ C (µ,µ1 ,m,ℓ) ,
(14)
where C (µ, µ1 , m, ℓ) = A0 (m) µ + A1 (m) µ1 + P (m, ℓ) being A0 (m) the number of evaluations of the scalar functions (F1 , . . . , Fm ) used in the evaluation of F and [y, x; F ]; A1 (m) is the number of evaluations of scalar functions of F ′ , namely, ∂ Fi , 1 ≤ i, j ≤ m. The function P (m, ℓ) represents the number of products needed per iteration, ℓ is the ratio between prod∂x j
ucts and quotients and µ and µ1 are ratios between products and evaluations required to express the value of C (µ, µ1 , m, ℓ) in terms of products (see [7,18]). To evaluate the divided difference operator we will use two different schemes (see [19] and the references therein). Namely, setting
∆j Fi (u, v) = Fi (u1 , . . . , uj−1 , uj , vj+1 , . . . , vm ) − Fi (u1 , . . . , uj−1 , vj , vj+1 , . . . , vm ) / (uj − vj ), we define for 1 ≤ i, j ≤ m,
[y, x; F ](i 1j ) = ∆j Fi (y, x), [y, x; F ](i 2j ) =
1 2
(15)
∆j Fi (y, x) + ∆j Fi (x, y) .
(16)
We point out that expression (16) presents a definition of a computational divided difference operator that led us to overcome the problem of losing convergence order that may appear when the iterative method is applied to some set of functions where the theoretical order does not agree with the one obtained using (15) as we will see later on. Hereafter we will use Φ4 when we choose λ = 0 and ν = 1 in (12), and the new iterative method Φ5 choosing r = 1 and s = −1 in (13). The cost Ci of methods Φi for i = 1, 2, 4, 5 is in the Table 2 and for the method Φ3 is m 6
(2m2 + 6mµ1 + 9m + 3ℓm + 12µ + 9ℓ − 11).
Note that in Table 2 the computational cost of Φ2 (frozen Kurchatov) and Φ4 (frozen Steffensen) is the same. Since ρ4 > ρ2 , then we conclude that Φ4 is more efficient than Φ2 . In general, to compare the efficiency of the other pairs of iterative methods Φi and Φj , we use the ratio Ri,j =
log CEIi (µ, µ1 , m, ℓ) log CEIj (µ, µ1 , m, ℓ)
=
log(ρi )Cj (µ, µ1 , m, ℓ) log(ρj )Ci (µ, µ1 , m, ℓ)
.
(17)
If Ri,j > 1 then the iterative method Φi is more efficient than Φj . From the preceding equation (17) and taking into account that the boundary between two computational efficiencies is attained by Ri,j = 1, then this boundary is expressed by the equation of µ written as a function of µ1 , ℓ and m. (k) Now, comparing by pairs of indices i, j the values of Ri,j and denoting by CEIi for 1 ≤ i ≤ 5, where k = 1, 2 when (15) and (16) are used, we state the following two results whose proofs are easy and they are left to the reader.
758
M. Grau-Sánchez et al. / Journal of Computational and Applied Mathematics 255 (2014) 753–764
Theorem 3.1. For all µ > 0, m ≥ 2, ℓ ≥ 1 and k = 1, 2 we have (k)
CEI4
> CEI2(k) > CEI1(k) and CEI4(k) > CEI5(k) > CEI1(k) .
This result does not depend on the divided difference scheme. (k)
(k)
(a) For all m ≥ 4, µ > 0 and ℓ ≥ 1 we have CEI5 > CEI2 , for k = 1, 2. (b) For the divided difference scheme [y, x; F ](1) and m = 2, 0 < µ ≤ 2.303791 ℓ − 1.767930, ℓ ≥ 1, or m = 3, 0 < µ ≤ (1) (1) 6.293914 ℓ − 0.300676, ℓ ≥ 1, we obtain CEI5 > CEI2 . (c) For the divided difference scheme [y, x; F ](2) and m = 2, 0 < µ ≤ 4.720435 ℓ− 3.622464, ℓ ≥ 1, or m = 3, µ > 0, ℓ ≥ 1, (2) (2) we get CEI5 > CEI2 . We point out that Φ4 (frozen Steffensen) from a theoretical point of view presents the highest efficiency among the freederivative iterative methods considered in this work. Likewise, the results for Φ3 (frozen Newton) are stated in the following. Theorem 3.2. For ℓ ≥ 1 holds (a) CEI1 > CEI3 , if m µ1 >
1 2Cm2 + 6Aµm + 3(2A + C )ℓm + 9Cm + 6(C − B)µ + 9C ℓ − 11C ,
6A
using [y, x; F ](1) ,
and m µ1 >
1 2Cm2 + 12Aµm + 3(2A + C )ℓm + 9Cm − 12Bµ + 9C ℓ − 11C , 6A
using [y, x; F ](2) ,
where A = log ρ3 , B = log ρ1 and C = log (ρ3 /ρ1 ). (b) CEI2 > CEI3 , if m µ1 >
1 2Em2 + 6Aµm + 3(2A + E )ℓm + 9Em + 12E µ + 9E ℓ − 11E , 6D
using [y, x; F ](1) ,
and m µ1 >
1 2Em2 + 12Aµm + 3(2A + E )ℓm + 9Em + 6(E − D)µ + 9E ℓ − 11E ,
6D
using [y, x; F ](2) ,
where A = log ρ3 , D = log ρ2 and E = log (ρ3 /ρ2 ). (c) CEI4 > CEI3 for µ1 > µ + ℓ using [y, x; F ](1) and mµ1 > 2mµ + mℓ − µ using [y, x; F ](2) . (d) CEI5 > CEI3 for mµ1 > mµ + µ + mℓ + 1 using [y, x; F ](1) and mµ1 > 2mµ + mℓ + 1 using [y, x; F ](2) . The boundary between the regions of two computational efficiencies, Ri,j = 1, is a plane in case (c) with [y, x; F ](1) , a hyperbolic cylinder in case (d) with [y, x; F ](2) and a hyperbolic paraboloid (see [18]) in the other cases. 4. Numerical results In this section numerical experiments are performed to illustrate the results previously presented. We begin with some comparisons between the theoretical divided difference operator and the corresponding first computational operator given in (15). We observe that in some situations a losing of convergence order occurs. More precisely, we will see that iterative methods Φ2 and Φ5 have a lower order. Finally, two numerical examples are also given. 4.1. Theoretical considerations Using the Taylor expression of the theoretical divided difference operator [ x + h, x; F ] we have
[ x + h, x; F ]ij =
1
Dj Fi (x + th) dt
0
= Dj Fi (x) +
m 1
2 k=1
Dkj Fi (x) hk +
m 1
6 k,ℓ=1
Dk ℓj Fi (x) hk hℓ + O(h3 ).
Likewise, for the computational divided difference operator [ x + h, x; F ](1) we get
[ x + h, x; F ](ij1) = Dj Fi (x) +
j −1 k=1
Dkj Fi (x) hk +
1 2
Dj j Fi (x) hj + O(h2 ).
(18)
M. Grau-Sánchez et al. / Journal of Computational and Applied Mathematics 255 (2014) 753–764
From (18) we obtain [ x + h, x; F ](1) = F ′ (x) + M h + O(h2 ), instead of [ x + h, x; F ] = F ′ (x) +
M ̸=
1 2
F (x) because M depends on the function F (see [19] and the references therein).
1 2
759
F ′′ (x) + O(h2 ). In general,
′′
Setting Ψk = [ xn + hk , xn ; F ](1) , 1 ≤ k ≤ 4, we have
en−1 − en , where Ψ1 = [ xn−1 , xn ; F ](1) which is the operator of Φ1 defined in (8), 2(en−1 − en ), where Ψ2 = [ xn−1 , 2xn − xn−1 ; F ](1) which is given by Φ2 in (10), hk = −F (xn ) = −Γ (en + o(en )) = −e˜n + o(en ), where Φ3 = [ xn , xn + F (xn ); F ](1) , −2qn = −2(en−1 − en )T (en−1 − en ) = −2e2n−1 + o(e2n−1 ), if Φ4 = [ xn − qn , xn + qn ; F ](1) . Taking into account that F ′ (xn ) = Γ I + 2 A2 en + O(e2n ) and M = Γ (N + O(en )), from Ψk = Γ I + 2 A2 en + N hk + O(h2k ) , (k)
we get Ψk−1 = (I − 2 A2 en − N hk + o(en−1 )) Γ −1 . The error expression corresponding to the first step, zn − α , defined in (8), (10), (12) and (13) is
En(k) = zn(k) − α = en − Ψk−1 Γ en + O(e2n ) = (2 A2 en + N hk ) en + o(hk en ) = [ O(en−1 en ), O(en−1 en ), O(e2n ), O(e2n−1 en ) ],
for 1 ≤ k ≤ 4.
Moreover, the final error using the inverse operator of Ψk is
(k) en+1 = x(nk) − α = En(k) − Ψk−1 Γ En(k) + O((En(k) )2 ) = (2 A2 en + N hk ) En(k) + o(hk En(k) ) = [ O(e2n−1 en ), O(e2n−1 en ), O(e3n ), O(e4n−1 en ) ],
for 1 ≤ k ≤ 4.
(19)
From the results of (19) we may conclude methods Φ2 (k = 2) and Φ5 (k =4) for some set of functions the that iterative
√
17 /2 to ρ˜ 2 = 2 and from ρ5 = 3 to ρ˜ 5 = 1 +
local convergence order go from ρ2 = 1 +
√
17 /2. Furthermore, we
also observe that Φ1 (k = 1, frozen Secant) and Φ4 (k = 3, frozen Steffensen) keep the order independently of the function tested. Given a function F to obtain the maximum order of convergence it may be necessary to use the operator [y, x; F ](2) given in (16) for Φ2 and Φ5 . This requires a higher computational cost; while for Φ1 and Φ4 we can carry out [y, x; F ](1) defined (1) (2) (1) (2) in (15). Namely, the comparison of CEI1 , CEI2 , CEI4 and CEI5 is done in the following theorem. Theorem 4.1. For all µ > 0, m ≥ 2 and ℓ ≥ 1 we have (1)
CEI4 (2)
(a) CEI5
> CEI1(1) ,
(1)
CEI4
> CEI2(2) and CEI4(1) > CEI5(2) .
> CEI2(2) , (see Theorem 3.1(a) and (c)).
(b) For m ≥ 2, ℓ ≥ 1 and 0 < µ <
1 2m2 +9(1+ℓ)m+9ℓ−11 6 ( B −1)m−1 A
(c) For m ≥ 2, ℓ ≥ 1 and 0 < µ <
(2)
we obtain CEI2
D 2 1 2m +9(1+ℓ)m+9ℓ−5−6 C B 6 −1 (m+1)
> CEI1(1) , where A = log (ρ2 /ρ1 ) and B = log ρ1 . (2)
we obtain CEI5
> CEI1(1) , where C = log(ρ5 /ρ1 ) and D = log ρ5 .
C
For ℓ = 1.731, Fig. 1 illustrates the following four zones
(1): CEI4(1) > CEI1(1) > CEI2(2) > CEI5(2) (2): CEI4(1) > CEI1(1) > CEI5(2) > CEI2(2)
(20)
(3): CEI4(1) > CEI5(2) > CEI1(1) > CEI2(2) (4): CEI4(1) > CEI5(2) > CEI2(2) > CEI1(1) .
Finally, using the operator [y, x; F ](1) for all four methods and recalling that the convergence orders of Φ2 and Φ5 are √ ρ2 = 2 and ρ5 = 1 + 17 /2 respectively, we state the following theorem, where the CEIs with these two orders of (1)
(1)
2 and CEI 5 . local convergence are denoted by CEI Theorem 4.2. For all µ > 0, m ≥ 2 and ℓ ≥ 1 we have (1)
CEI4
(1)
2 > CEI1(1) > CEI
and
(1)
CEI4
(1)
(1)
5 > CEI 2 . > CEI
760
M. Grau-Sánchez et al. / Journal of Computational and Applied Mathematics 255 (2014) 753–764
m Fig. 1. Boundary curves Ri,j = 1 for ℓ = 1.731 and m = 2, 3, 4, 5 (see (20)).
(1)
(1)
5 > CEI for Finally, we obtain CEI 1 (a) (b) (c) (d)
all m ≥ 5 and ℓ ≥ 1; m = 2, ℓ ≥ 1 and 0 < µ < 1.729473 ℓ − 0.115673; m = 3, ℓ ≥ 1 and 0 < µ < 3.745439 ℓ + 1.7888879; m = 4, ℓ ≥ 1 and 0 < µ < 12.459539 ℓ + 11.128903.
The proof as before is left to the reader. 4.2. Numerical examples The numerical computations were performed using the MPFR library of C++ multi-precision arithmetics [20,21] with 4096 digits of mantissa. All algorithms were compiled by g++(4.2.1) for i686-apple-darwin1 with libgmp (v.5.0.2) R and libmpfr (v.3.1.0) libraries in a processor Intel⃝ Xeon E5620, 2.4 GHz (64-bit machine). In this machine the quotient and product ratio is ℓ = 1.731. In each example the starting point is the same for all methods tested. The classical stopping criteria ∥eI ∥ = ∥xI − α∥ > ε and ∥eI +1 ∥ < ε , with ε = 10−ν , where ν = 4096, are replaced by ∥˘eI ∥ > 10−η and ∥˘eI +1 ∥ < 10−η , where η = [ν(ρ − 1)/ρ ] and ∥ · ∥ = ∥ · ∥∞ . Moreover, e˘ n is obtained by
e˘ n =
F r ( xn )
Fr (xn−1 )
.
(21)
1≤r ≤m
Note that this criterion is independent of the knowledge of the root (see [22]). To evaluate the numerical efficiency of each method we compute its factor κ˜ (see [8]). That is, for each iteration k we get ˜ (Dk ), where the precision (number of correct decimals) Dk in time Θ Dk ≈ − log10 ∥ek ∥ ≈ −
ρ log10 ∥˘ek ∥. ρ−1
˜ (Dk )), 1 ≤ k ≤ I, in the least-squares sense by a polynomial of degree one, where the We approximate the pairs (log Dk , Θ slope κ˜ is computed. Namely, κ˜ is a coefficient measuring the time of execution in the function of the approximate number of correct decimals, say ˜ (DI ) = κ( Θ ˜ log DI − log D0 ). defined by Taking into account this slope we compare the time factor TF with the computed time factor TF = TF
˜ (DI ) Θ tp log q
=
κ˜ tp
≈ TF =
1 log CEI
,
where tp is the necessary time spent by one product and q = DI /D0 . In these experiments we calculate the computational order of convergence ρ˘ , called PCLOC, defined in [22] by
ρ˘ =
log ∥F (xI )∥ log ∥F (xI −1 )∥
.
M. Grau-Sánchez et al. / Journal of Computational and Applied Mathematics 255 (2014) 753–764
761
Table 3 Numerical results for the nonlinear system (22).
(1)
Φ1 (1) Φ2 Φ3 (1 ) Φ4 (1 ) Φ5
ρ
TF
I
T
qI
κ˜
TF
rTF
2.0000
3204.73
2.5616 3.0000
2922.52 1509.05
10
229.05
2897
80.31
3318.73
3.56
7 6
193.60 98.80
1631 1852
72.82 36.86
3009.11 1522.95
2.96 0.92
3.0000
2502.21
7
201.93
3212
62.70
2590.82
3.54
3.0000
2988.75
7
232.51
3949
74.32
3071.26
2.76
Table 4 Numerical results for the nonlinear system (23).
(1)
Φ1 (1) Φ2 (2 ) Φ2 Φ3 (1 ) Φ4 (1 ) Φ5 (2 ) Φ5
ρ
TF
I
T
qI
κ˜
TF
rTF
2.0000
4199.31
9
264.86
2113
101.77
4205.51
0.15
2.0000
5182.73
9
327.45
2737
116.56
4816.71
7.06
2.5616 3.0000
5326.98 2018.04
7 6
350.72 131.10
2454 2628
130.19 48.40
5379.92 2000.20
0.99 0.88 0.26
3.0000
3286.60
6
224.42
1740
79.75
3295.27
2.5616
4566.81
7
299.18
2257
11.82
4620.64
1.18
3.0000
5204.29
6
337.70
1950
126.86
5242.17
0.73
If ρ = ρ˘ ± ∆ρ˘ , where ρ is the local order of convergence and ∆ρ˘ is the error of PCLOC, then we get ∆ρ˘ < 10−3 . This fact means that in all computations of PCLOC we obtain at least 3 significant digits and this result is a good test of the local convergence orders of the family of iterative methods worked in this paper. Next we present two examples. The first one does not have any problem when we apply the definition [y, x; F ](1) given in (15). The second example needs the definition (16) for the Φ2 and Φ5 methods to attain the local convergence order that theoretical results state. Tables 3 and 4 show the results obtained for the iterative methods Φk , 1 ≤ k ≤ 5. In each table it is shown the local order of convergence ρ , the measure time factor TF defined as TF = 1/ log CEIk , the necessary iteration number I, the necessary time T in milliseconds to reach the I-iteration, the number of correct decimals in xI , say qI , the slope κ and the computed . Furthermore, the last column shows the percentage of relative error rTF between TF and TF . time factor TF In the first example it is considered a system F (x1 , x2 , x3 ) = 0 with three exponential equations defined by
x2 + x3 − e−x1 , F (x1 , x2 , x3 ) = x1 + x3 − e−x2 , x + x − e−x3 , 1 2
(22)
where (m, µ, µ1 = µ/3) = (3, 76.4, 25.5). Setting x−1 = (−0.83, 1.14, 1.14)t and x0 = (−0.84, 1.15, 1.15)t the following approximation of the root α ≈ (−.83202504 . . . , 1.1489838 . . . , 1.1489838 . . .)t with the prefixed precision is obtained. Notice that the preceding results agree with the theoretical ones exposed previously. Observing the computed time factor of the different methods we conclude that they also agree with the computational efficiency index CEI or the time factor TF TF . Moreover, the inequalities (1)
CEI3 > CEI4
> CEI2(1) > CEI5(1) > CEI1(1) ,
are in concordance with the ones in Theorems 3.1 and 3.2. In Fig. 2 the different slopes κ˜ for each method are drawn. Note that the minimum slope and time correspond to the Newton method Φ3 followed by the Steffensen method Φ4 . A second example where the theoretical local convergence orders are different from the ones studied in Section 2 and presented in Table 1 is x1 − cos(x1 − (x2 + x3 )), x2 − cos(x2 − (x1 + x3 )), x3 − cos(x3 − (x1 + x2 )),
F (x1 , x2 , x3 ) =
(23)
where (m, µ, µ1 ) = (3, 100.8, 34.5). Using the initial values x−1 = (0.544, 0.544, 0.996)t and x0 = (0.543, 0.543, 0.995)t we have the following approximation of the root
α ≈ (0.5438500415, 0.5438500415, 0.9957781534)t . Recall that, in this case, it is necessary to use the second computational definition of the divided difference operator. Table 4 shows the results of the system (23). Notice that there are more rows than in Table 3 because we have studied the iterative methods using the two divided difference operators for the iterative functions Φ2 and Φ5 . Also note that CEI , TF
762
M. Grau-Sánchez et al. / Journal of Computational and Applied Mathematics 255 (2014) 753–764
Fig. 2. Precisions and computational time of the five methods for example (22).
Fig. 3. Precisions and computational time of the five methods with maximum order for example (23).
are ordered in the same way as all the theorems stated in this paper. Furthermore, we have and TF (1)
CEI3 > CEI4
> CEI1(1) > CEI5(1) > CEI2(1) > CEI5(2) > CEI2(2) .
In Fig. 3 the different slopes κ˜ for each method are drawn. Observe that the minimum slope and time correspond to the (1) Newton method Φ3 , followed by the Steffensen method Φ4 . 5. Application Hereafter, we apply the technique presented previously to the well-known Hammerstein integral equation [23,24]. Namely, x(s) = 1 +
1 2
1
G(s, t ) x(t )3 dt ,
s ∈ [0, 1],
0
where x ∈ C [0, 1], t ∈ [0, 1], and the kernel G is the Green function G(s, t ) =
(1 − s) t , s (1 − t ),
t ≤ s, s < t.
M. Grau-Sánchez et al. / Journal of Computational and Applied Mathematics 255 (2014) 753–764
763
Table 5 Numerical results for the nonlinear system (25).
(1)
Φ1 (1) Φ2 Φ3 (1) Φ4 (1) Φ5
ρ
TF
I
2.0000
4089.12
2.5616 3.0000
3228.73 1241.12
3.0000 3.0000
T
qI
10
192.85
2866
8 7
152.03 63.01
4035 3795
2764.39
6
124.16
1795
2965.60
6
123.23
1402
We begin by writing the preceding integral equation as R(x) = 0, where R : C [0, 1] → C [0, 1] and
[R(x)](s) = x(s) − 1 −
1 2
1
G(s, t ) x(t )3 dt ,
s ∈ [0, 1].
(24)
0
To obtain the discretization of the integral we can use Gauss–Legendre formulae. That is, 1
f (t ) dt ≈ 0
m
ϖj f (tj ),
j =1
where the abscissas tj and the weights ϖj were computed with maximum precision for m = 8 in [25]. Next, we present the discrete version of (24). Indeed, denoting the approximation of x(ti ) by xi (1 ≤ i ≤ 8), we get the following system of nonlinear equations xi = 1 +
8 1
2 j=1
bij x3j ,
where bij =
ϖj tj (1 − ti ) ϖj ti (1 − tj )
if j ≤ i, if j > i,
1 ≤ i ≤ 8,
(25)
with solution
α ≈ (1.005771, 1.027442, 1.055178, 1.074414, 1.074414, 1.055178, 1.027442, 1.005771)t .
(26)
In this case the values of parameters are (m, µ, µ1 ) = (8, 11, 11/8). Taking x−1 = (1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0)t and x0 = (1.1, 1.1, 1.1, 1.1, 1.1, 1.1, 1.1, 1.1)t we show and com(1) pare the convergence of the methods Φ k , 1 ≤ k ≤ 5, towards the root α (26) in Table 5. 6. Concluding remarks In this paper we present a generalized family of two step Secant-like methods with frozen operator for solving nonlinear equations. A local convergence analysis of this family, where five specific cases are considered, is studied in detail. Namely, using earlier methods such as Secant, Kurchatov, Newton, Steffensen and other new algorithms, the family of iterative schemes is built up and their computational efficiency is also given. Numerical examples using multiple precision and a stopping criterion are implemented without using any known root. Finally, a study comparing the order, efficiency and elapsed time of the methods suggested supports the theoretical results claimed. We point out that the real efficiency of each scheme is related to the slope of the lines shown in the figures presented. Acknowledgement The research was supported partially by the project MTM2011-28636-C02-01. References [1] J.A. Ezquerro, M.A. Hernández, An optimization of Chebyshev’s method, J. Complexity 25 (2009) 343–361. [2] M. Grau, J.L. Díaz-Barrero, An improvement of the Euler–Chebyshev iterative method, J. Math. Anal. Appl. 315 (2006) 1–7. [3] M. Grau-Sánchez, J.M. Gutiérrez, Some variants of the Chebyshev–Halley family of methods with fifth order of convergence, Int. J. Comput. Math. 87 (2010) 818–833. [4] J.A. Ezquerro, M.A. Hernández, M.J. Rubio, Secant-like methods for solving nonlinear integral equations of the Hammerstein type, J. Comput. Appl. Math. 115 (1–2) (2000) 245–254. [5] I.K. Argyros, J.A. Ezquerro, J.M. Gutiérrez, M.A. Hernández, S. Hilout, On the semilocal convergence of efficient Chebyshev–Secant-type methods, J. Comput. Appl. Math. 235 (2011) 3195–3206. [6] I.K. Argyros, H. Ren, Efficient Steffensen-type algorithms for solving nonlinear equations, Int. J. Comput. Math. 90 (2013) 691–704. [7] M. Grau-Sánchez, A. Grau, M. Noguera, Frozen divided difference scheme for solving systems of nonlinear equations, J. Comput. Appl. Math. 235 (2011) 1739–1743. [8] M. Grau-Sánchez, M. Noguera, A technique to choose the most efficient method between secant method and some variants, Appl. Math. Comput. 218 (2012) 6415–6426.
764
M. Grau-Sánchez et al. / Journal of Computational and Applied Mathematics 255 (2014) 753–764
[9] J.A. Ezquerro, A. Grau, M. Grau-Sánchez, M.A. Hernández, M. Noguera, Analysing the efficiency of some modifications of the secant method, Comput. Math. Appl. 64 (2012) 2066–2073. [10] J.A. Ezquerro, M. Grau-Sánchez, M.A. Hernández, Solving non-differentiable equations by a new one-point iterative method with memory, J. Complexity 28 (2012) 48–58. [11] J.A. Ezquerro, A. Grau, M. Grau-Sánchez, M.A. Hernández, Construction of derivative-free iterative methods from Chebyshev’s method, Anal. Appl. 11 (2013) 1350009 (16 pp.). [12] V.A. Kurchatov, On a method of linear interpolation for the solution of functional equations, Dokl. Akad. Nauk SSSR 198 (3) (1971) 524–526 (in Russian); translation in Sov. Math. Dokl. 12 (1971) 835–838. [13] J.A. Ezquerro, M. Grau-Sánchez, M.A. Hernández, M. Noguera, Semilocal convergence of secant-like methods for differentiable and nondifferentiable operator equations, J. Math. Anal. Appl. 398 (2013) 100–112. [14] S. Amat, S. Busquier, J.M. Gutiérrez, Third-order iterative methods with applications to Hammerstein equations: a unified approach, J. Comput. Appl. Math. 235 (2011) 2936–2943. [15] I.K. Argyros, On the Newton–Kantorovich hypothesis for solving equations, J. Comput. Appl. Math. 169 (2004) 315–332. [16] S. Amat, S. Busquier, V. Candela, A class of quasi-Newton generalized Steffensen methods on Banach spaces, J. Comput. Appl. Math. 149 (2002) 397–406. [17] V. Alarcón, S. Amat, S. Busquier, D.J. López, A Steffensen’s type method in Banach spaces with applications on boundary-value problems, J. Comput. Appl. Math. 216 (2008) 243–250. [18] M. Grau-Sánchez, A. Grau, M. Noguera, On the computational efficiency index and some iterative methods for solving systems of nonlinear equations, J. Comput. Appl. Math. 236 (2011) 1259–1266. [19] M. Grau-Sánchez, M. Noguera, S. Amat, On the approximation of derivatives using divided difference operators preserving the local convergence order of iterative methods, J. Comput. Appl. Math. 237 (2013) 363–372. [20] L. Fousse, G. Hanrot, V. Lefèvre, P. Pélissier, P. Zimmermann, MPFR: a multiple-precision binary floating-point library with correct rounding, ACM Trans. Math. Software 33 (2) (2007) 15. Art. 13. [21] The GNU MPFR library 3.1.0. Available in http://www.mpfr.org. [22] M. Grau-Sánchez, A. Grau, M. Noguera, J.R. Herrero, A study on new computational local orders of convergence, Appl. Math. Lett. 25 (2012) 2023–2030. [23] J.M. Ortega, W.C. Rheinboldt, Iterative Solution of Nonlinear Equations in Several Variables, Academic Press, New York, 1970. [24] A.D. Polyanin, A.V. Manzhirov, Handbook of Integral Equations, CRC Press, Boca Raton, 1998. [25] J.A. Ezquerro, M. Grau-Sánchez, A. Grau, M.A. Hernández, M. Noguera, N. Romero, On iterative methods with accelerated convergence for solving systems of nonlinear equations, J. Optim. Theory Appl. 151 (2011) 163–174.
Contents lists available at ScienceDirect
Journal of Computational and Applied Mathematics journal homepage: www.elsevier.com/locate/cam
On the local convergence of a family of two-step iterative methods for solving nonlinear equations Miquel Grau-Sánchez ∗ , Miquel Noguera, José L. Diaz-Barrero Technical University of Catalonia, Department of Applied Mathematics II and III, Jordi Girona 1-3, Omega, 08034 Barcelona, Spain
article
info
Article history: Received 27 February 2013 Received in revised form 14 May 2013 MSC: 41A25 65H10 Keywords: Order of convergence Nonlinear equations Iterative methods Efficiency
abstract A local convergence analysis for a generalized family of two step Secant-like methods with frozen operator for solving nonlinear equations is presented. Unifying earlier methods such as Secant’s, Newton, Chebyshev-like, Steffensen and other new variants the family of iterative schemes is built up, where a profound and clear study of the computational efficiency is also carried out. Numerical examples and an application using multiple precision and a stopping criterion are implemented without using any known root. Finally, a study comparing the order, efficiency and elapsed time of the methods suggested supports the theoretical results claimed. © 2013 Elsevier B.V. All rights reserved.
1. Introduction There are a great variety of iterative methods for solving a system of nonlinear equations F (x) = 0, where F : D ⊆ Rm → R , and D is a non-empty open convex subset of Rm that contains a simple root α of F . A classical iterative process for solving nonlinear equations is Chebyshev’s method (see [1–3]) m
x ∈ D, 0 yn = xn − F ′ (xn )−1 F (xn ) xn+1 = yn − 1 F ′ (xn )−1 F ′′ (xn ) (yn − xn )2 , 2
n ≥ 0.
The above one-point iterative scheme depends explicitly on the two first derivatives of F . In [1], Ezquerro and Hernández present some modifications in Chebyshev’s method by reducing in one the number of evaluations of the first derivative and maintaining third-order of convergence. It has the following form:
x ∈ D, 0 zn = xn − a F ′ (xn )−1 F (xn ), xn+1 = xn − 1 F ′ (xn )−1 (a2 + a − 1) F (xn ) + F (zn ) , a2
n ≥ 0.
Using the well-known Secant method [4], in [5] a generalization of it employing the divided difference operator of order one (namely, Bn = [xn−1 , xn ; F ]) that substitutes the derivative of F (F ′ (xn ) ≡ [xn , xn ; F ]) is given. The authors call this family
∗
Corresponding author. Tel.: +34 93 413 79 82. E-mail address: [email protected] (M. Grau-Sánchez).
0377-0427/$ – see front matter © 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.cam.2013.06.043
754
M. Grau-Sánchez et al. / Journal of Computational and Applied Mathematics 255 (2014) 753–764
the Chebyshev–Secant-type method and it is defined by
x−1 , x0 ∈ D, 1 zn = xn − a B− n F (xn ), −1 xn+1 = xn − Bn (b F (xn ) + c F (zn )) ,
n ≥ 0,
where a, b, c are non-negative parameters to be chosen so that the sequence {xn } converges to α with maximum local order of convergence. The work presented in [6] analyzes free-derivative iterative processes considering the operator Cn = [xn , L(xn ); F ] and they are called the Steffensen-type method:
x0 ∈ D , zn = xn − a Cn−1 F (xn ), xn+1 = xn − Cn−1 (b F (xn ) + c F (zn )) , m
n ≥ 0,
m
where L : D ⊆ R → R . Our goal in the present paper is to unify these methods by considering the following iterative family
x −1 , x 0 ∈ D , zn = xn − a Θn−1 F (xn ), xn+1 = xn − Θn−1 (b F (xn ) + c F (zn )) ,
(1) n ≥ 0,
where Θn = [G(xn−1 , xn ), H (xn−1 , xn ); F ] and G, H : D ⊆ Rm × Rm → Rm . As can be observed, the above family includes as particular cases the preceding works [5,6] and other well-known algorithms. Furthermore, (1) also offers us the possibility to suggest new methods as we will see later on. 2. Local order of convergence To obtain explicitly an expression of the inverse operator of Θn , we denote by J indistinctly both G or H. On account that a necessary condition about functions G and H is: G(α, α) = H (α, α) = α , then from Taylor’s formulae we get the error expression
εJ = J (xn−1 , xn ) − α = J1 en−1 + J2 en + 1/2 J11 e2n−1 + J22 e2n + 2J12 en−1 en + · · · , ∂J
∂J
∂ 2J
where J1 = ∂ x (α, α), J2 = ∂ x (α, α), J11 = 2 (α, α), etc. n ∂ xn−1 n−1 The expression of Θn is given in the previous works appeared in [7,8]. Namely,
Θn = [G(xn−1 , xn ), H (xn−1 , xn ); F ] = Γ I + A2 (εG + εH ) + A3 (εG2 + εG εH + εH2 ) + A4 (εG3 + εG2 εH + εG εH2 + εH3 ) + · · · , where Γ = F ′ (α) ∈ L(Rm ), Ak = k1! Γ −1 F (k) (α), Ak ∈ Lk (Rm , Rm ) and Ak is k-symmetric. From the preceding follows
Θn−1 = I − A2 (εG + εH ) − (A3 − A22 ) (εG2 + εG εH + εH2 ) + A22 εG εH + A2 (εG + εH ) · A3 (εG2 + εG εH + εH2 ) + A3 (εG2 + εG εH + εH2 ) · A2 (εG + εH ) − A4 (εG3 + εG2 εH + εG εH2 + εH3 ) Γ −1 , where A22 εG2 = (A2 εG )2 . Setting En = zn − α, en = xn − α and using Taylor’s development yield F (xn ) = Γ en + A2 e2n + o(e2n ) .
Subtracting α from both sides of the first equation of (1), we obtain En = en − a Θn−1 Γ en + A2 e2n + o(e2n )
= (1 − a) en + a A2 (G1 + H1 ) en−1 en 1 2 2 2 +a A2 (G11 + H11 ) + A3 (G1 + H1 + G1 H1 ) − (A2 (G1 + H1 )) e2n−1 en 2
+ a A2 (G2 + H2 − I ) e2n + o(e2n−1 en , e2n ),
(2)
where the following notation was used: A2 (G1 + H1 ) en−1 en = A2 (G1 en−1 + H1 en−1 ) en , A2 (G11 + H11 ) e2n−1 en = A2 (G11 e2n−1 + H11 e2n−1 ) en , A3 (G21 + H12 + G1 H1 ) e2n−1 en = A3 (G1 en−1 )(G1 en−1 ) + (H1 en−1 )(H1 en−1 ) − (G1 en−1 )(H1 en−1 ) en ,
M. Grau-Sánchez et al. / Journal of Computational and Applied Mathematics 255 (2014) 753–764
755
and
(A2 (G1 + H1 ))2 e2n−1 en = (A2 (G1 en−1 + H1 en−1 )) (A2 (G1 en−1 + H1 en−1 )) en . Note that with this notation we may order the terms of the error vector En according to their order of convergence. Recall γ β that if we have the term K en−1 en , its local order of convergence is the unique positive root of the associated polynomial 2 equation t − β t − γ = 0. Then (2) is the expression of the error En up to second order. From Taylor’s formulae F (zn ) = Γ (En + o(En )) and subtracting α from both sides of the second equation of (1), we have en+1 = en − Θn−1 b Γ
en + A2 e2n + o(e2n ) + c Γ (En + o(En ))
= (1 − b − c + ac ) en + (b + c − 2ac ) A2 (G1 + H1 ) en−1 en + o(en−1 en ).
(3)
Annulling the coefficients of the terms of low order of (3) yields
1 − b − c + ac = 0, b + c − 2ac = 0,
and we get b = 2 −
1 a
and c =
1 . a
Substituting these values in (3), we obtain
en+1 = (A2 (G1 + H1 ))2 e2n−1 en + (1 − a) A2 e2n + o(e2n−1 en , e2n ).
(4)
Hereafter, a = 1 and G1 + H1 ≡ 0 will be used in order that the family (1) becomes a super-quadratic iterative method. Namely, from (1) and taking into account the previous considerations, we have
x−1 , x0 ∈ D, zn = xn − Θn−1 F (xn ), xn+1 = zn − Θn−1 F (zn ),
(5) n ≥ 0.
Note that for a = b = c = 1, we obtain an iterative method of two steps very similar ones with the same operator. This fact is the reason why this kind of methods are called frozen iterative methods, since the operator is repeated. If we consider (5), where G1 + H1 ≡ 0, then (2) becomes
En =
1 2
A2 (G11 + H11 ) + A2 G21
+
1
A2 (G111 + H111 ) +
6
e2n−1 en + A2 (G2 + H2 − I ) e2n
1
A3 G1 (G11 − H11 ) e3n−1 en
2
+ (A2 (G12 + H12 ) + A3 G1 (G2 − H2 )) en−1 e2n + o(en−1 e2n ), where A2 (G2 + H2 − I ) e2n = A2 (G2 en + H2 en − en ) en , A2 (G111 + H111 ) e3n−1 en = A2 (G111 e3n−1 + H111 e3n−1 ) en , A3 G1 (G11 − H11 ) e3n−1 en = A3 (G1 en−1 ) G11 e2n−1 − H11 e2n−1 en ,
A2 (G12 + H12 ) en−1 e2n = A2 (G12 en−1 en + H12 en−1 en ) en , and A3 G1 (G2 − H2 ) en−1 e2n = A3 (G1 en−1 )(G2 en − H2 en ) en . Furthermore, (4) becomes
e n +1 =
1 4
(A2 (G11 + H11 )) + A2 A3 (G11 + H11 ) 2
+ A2
1 2
G21
A2 (2(G2 + H2 ) − I ) (G11 + H11 ) +
+
2 A3 G21
A3 G21
e4n−1 en
(2(G2 + H2 ) − I ) e2n−1 e2n + o(e2n−1 e2n ),
where A2 A3 (G11 + H11 ) G21 e4n−1 en =
1 2
A2 (G11 e2n−1 + H11 e2n−1 )A3 (G1 en−1 )2
+
1 2
A3 G21
2
e4n−1 en = A3 (G1 en−1 )2
en
A3 (G1 en−1 )2 A2 (G11 e2n−1 + H11 e2n−1 ) en ,
and
A3 (G1 en−1 )2 en .
(6)
756
M. Grau-Sánchez et al. / Journal of Computational and Applied Mathematics 255 (2014) 753–764
√
The orders of the first two terms in (6) are 1 +
17 /2 = 2.562 and 1 +
√
3 = 2.732 respectively. But if we set G1 ≡ 0
and G11 + H11 ≡ 0 the two terms above described and also the terms of third order in e6n−1 en and e3n−1 e2n are zero. Therefore, the only nonzero term of third order is en+1 = A2 (G2 + H2 ) A2 (G2 + H2 − I ) e3n + o(e3n ).
(7)
Now we are ready to state the following result. Theorem 2.1. If a = b = c = 1, G1 ≡ H1 ≡ 0 and G11 + H11 ≡ 0, then iterative method (1) has R-order of convergence at least 3. Furthermore, the vectorial equation of the error is given in (7). Finally, we conclude this section presenting five particular cases of the preceding results. Case 1. From (4), if b = 2 − 1/ a and c = 1/ a, family (1) is a set of iterative methods with memory of at least second order. An element of this family can be found in [7] where the scheme is called the frozen secant method. Here, it is set as Θ = [G, H ] = [xn−1 , xn ; F ] and it is obtained as
x−1 , x0 ∈ D, zn = xn − [xn−1 , xn ; F ]−1 F (xn ), Φ1 : xn+1 = zn − [xn−1 , xn ; F ]−1 F (zn ),
(8) n ≥ 0,
where a = 1, but G1 + H1 ≡ 1 ̸= 0. Others works with this operator can be found in [9–11]. Case 2. Taking G(xn−1 , xn ) = γ xn−1 + (1 − γ )xn
and H (xn−1 , xn ) = δ xn−1 + (1 − δ)xn ,
(9)
if we impose that G1 + H1 ≡ 0, then we obtain γ +δ = 0, and from (6) we have a one-parametric family with super-quadratic order after adding the condition a = b = c = 1. For γ = 1 (δ = −1), the operator Θ becomes [xn−1 , 2xn − xn−1 ; F ] that is the operator used by Kurchatov [12]. This algorithm, called the frozen Kurchatov method, can be found in [13] and it is defined by
x−1 , x0 ∈ D, zn = xn − [xn−1 , 2xn − xn−1 ; F ]−1 F (xn ), Φ2 : xn+1 = zn − [xn−1 , 2xn − xn−1 ; F ]−1 F (zn ),
(10) n ≥ 0,
where a = b = c = 1, G1 + H1 ≡ 0, G11 + H ≡ 0. Notice that in general, in (9), G1 = γ ̸= 0. Hence, the local √ 11 convergence order of this method is 1 + 17 / 2 ≈ 2.56155. Observe that if (5) is an iterative method without memory, then G1 ≡ H1 ≡ 0. In other words, the iterative methods without memory will be super-quadratic if a = b = c = 1. Case 3. Setting γ = 0 (δ = 0) in (9), the operator considered in Case 2 becomes [xn , xn ; F ] ≡ F ′ (xn ) and all the hypotheses of Theorem 2.1 hold, because now G1 ≡ 0. This scheme is known as the frozen Newton method [1,14]:
x0 ∈ D , zn = xn − [F ′ (xn )]−1 F (xn ), Φ3 : xn+1 = zn − [F ′ (xn )]−1 F (zn ),
(11) n ≥ 0.
As we have set G = H = xn , we get G2 = H2 = I, and substituting in (7) we obtain the well known error equation of the frozen Newton method: en+1 = 2 A22 e3n + o(e3n ). Case 4. A family of iterative methods without memory was presented by many authors (see [15–17] and the references therein). Considering G = xn + λF (xn ) and H = xn + ν F (xn ), where λ, ν ∈ R, and on account that the hypotheses of the theorem hold, yields
x0 ∈ D, zn = xn − [xn + λF (xn ), xn + ν F (xn ); F ]−1 F (xn ), Φ4 : xn+1 = zn − [xn + λF (xn ), xn + ν F (xn ); F ]−1 F (zn ),
(12) n ≥ 0.
In this case, we have G2 = I + λ Γ and H2 = I + ν Γ . The error equation (7) becomes en+1 = A2 2 I + (λ + ν) Γ A2 I +
(λ + ν) Γ e3n + o(e3n ). These methods are called Steffensen-like methods.
Case 5. The last family considered in this paper is a set of algorithms that we call quadratic. The definition of a component (i) (i) (i) (i) (i) (i) function of G is given by Gi (xn−1 , xn ) = xn + r (xn − xn−1 )2 , likewise, H i (xn−1 , xn ) = xn + s (xn − xn−1 )2 , where x =
(x(1) , . . . , x(m) ) ∈ Rm , r , s ∈ R and 1 ≤ i ≤ m. Putting qn = (q(n1) , . . . , q(nm) ) ∈ Rm where q(ni) = (x(ni) − x(ni−) 1 )2 , the resulting new family of iterative methods is
x−1 , x0 ∈ D, zn = xn − [xn + r qn , xn + s qn ; F ]−1 F (xn ), Φ5 : xn+1 = zn − [xn + r qn , xn + s qn ; F ]−1 F (zn ),
(13) n ≥ 0.
M. Grau-Sánchez et al. / Journal of Computational and Applied Mathematics 255 (2014) 753–764
757
Table 1 Local convergence order ρi of the methods Φi for 1 ≤ i ≤ 5.
ρ1
ρ2
2
1+ 17 2
ρ3
ρ4
ρ5
3
3
3
√
Table 2 Cost of the methods Φi for i = 1, 2, 4, 5. Using [y, x; F ](1) defined in (15)
Φ1 Φ2 Φ4 Φ5
m 6 m 6 m 6 m 6
Using [y, x; F ](2) defined in (16)
(2m2 + 6mµ+ 9m + 9ℓm + 6µ+ 9ℓ− 11) (2m2 + 6mµ+ 9m + 9ℓm + 12µ+ 9ℓ− 11) (2m2 + 6mµ+ 9m + 9ℓm + 12µ+ 9ℓ− 11) (2m2 + 6mµ+ 9m + 9ℓm + 18µ+ 9ℓ− 5)
m 6 m 6 m 6 m 6
(2m2 + 12mµ + 9m + 9ℓm + 9ℓ − 11) (2m2 + 12mµ+ 9m + 9ℓm + 6µ+ 9ℓ− 11) (2m2 + 12mµ+ 9m + 9ℓm + 6µ+ 9ℓ− 11) (2m2 + 12mµ+ 9m + 9ℓm + 12µ+ 9ℓ− 5)
Taking into account that G1 ≡ H1 ≡ 0 and imposing that G11 + H11 ≡ 0, we have r + s = 0. Furthermore, from G2 + H2 = 2 I and (7) we obtain en+1 = 2 A22 e3n + o(e3n ). Finally, in Table 1 the local convergence orders of the iterative families Φi , 1 ≤ i ≤ 5 are shown. 3. Computational efficiency index The computational efficiency index (CEI ) and the computational cost per iteration (C ) are defined by 1
CEI (µ, µ1 , m, ℓ) = ρ C (µ,µ1 ,m,ℓ) ,
(14)
where C (µ, µ1 , m, ℓ) = A0 (m) µ + A1 (m) µ1 + P (m, ℓ) being A0 (m) the number of evaluations of the scalar functions (F1 , . . . , Fm ) used in the evaluation of F and [y, x; F ]; A1 (m) is the number of evaluations of scalar functions of F ′ , namely, ∂ Fi , 1 ≤ i, j ≤ m. The function P (m, ℓ) represents the number of products needed per iteration, ℓ is the ratio between prod∂x j
ucts and quotients and µ and µ1 are ratios between products and evaluations required to express the value of C (µ, µ1 , m, ℓ) in terms of products (see [7,18]). To evaluate the divided difference operator we will use two different schemes (see [19] and the references therein). Namely, setting
∆j Fi (u, v) = Fi (u1 , . . . , uj−1 , uj , vj+1 , . . . , vm ) − Fi (u1 , . . . , uj−1 , vj , vj+1 , . . . , vm ) / (uj − vj ), we define for 1 ≤ i, j ≤ m,
[y, x; F ](i 1j ) = ∆j Fi (y, x), [y, x; F ](i 2j ) =
1 2
(15)
∆j Fi (y, x) + ∆j Fi (x, y) .
(16)
We point out that expression (16) presents a definition of a computational divided difference operator that led us to overcome the problem of losing convergence order that may appear when the iterative method is applied to some set of functions where the theoretical order does not agree with the one obtained using (15) as we will see later on. Hereafter we will use Φ4 when we choose λ = 0 and ν = 1 in (12), and the new iterative method Φ5 choosing r = 1 and s = −1 in (13). The cost Ci of methods Φi for i = 1, 2, 4, 5 is in the Table 2 and for the method Φ3 is m 6
(2m2 + 6mµ1 + 9m + 3ℓm + 12µ + 9ℓ − 11).
Note that in Table 2 the computational cost of Φ2 (frozen Kurchatov) and Φ4 (frozen Steffensen) is the same. Since ρ4 > ρ2 , then we conclude that Φ4 is more efficient than Φ2 . In general, to compare the efficiency of the other pairs of iterative methods Φi and Φj , we use the ratio Ri,j =
log CEIi (µ, µ1 , m, ℓ) log CEIj (µ, µ1 , m, ℓ)
=
log(ρi )Cj (µ, µ1 , m, ℓ) log(ρj )Ci (µ, µ1 , m, ℓ)
.
(17)
If Ri,j > 1 then the iterative method Φi is more efficient than Φj . From the preceding equation (17) and taking into account that the boundary between two computational efficiencies is attained by Ri,j = 1, then this boundary is expressed by the equation of µ written as a function of µ1 , ℓ and m. (k) Now, comparing by pairs of indices i, j the values of Ri,j and denoting by CEIi for 1 ≤ i ≤ 5, where k = 1, 2 when (15) and (16) are used, we state the following two results whose proofs are easy and they are left to the reader.
758
M. Grau-Sánchez et al. / Journal of Computational and Applied Mathematics 255 (2014) 753–764
Theorem 3.1. For all µ > 0, m ≥ 2, ℓ ≥ 1 and k = 1, 2 we have (k)
CEI4
> CEI2(k) > CEI1(k) and CEI4(k) > CEI5(k) > CEI1(k) .
This result does not depend on the divided difference scheme. (k)
(k)
(a) For all m ≥ 4, µ > 0 and ℓ ≥ 1 we have CEI5 > CEI2 , for k = 1, 2. (b) For the divided difference scheme [y, x; F ](1) and m = 2, 0 < µ ≤ 2.303791 ℓ − 1.767930, ℓ ≥ 1, or m = 3, 0 < µ ≤ (1) (1) 6.293914 ℓ − 0.300676, ℓ ≥ 1, we obtain CEI5 > CEI2 . (c) For the divided difference scheme [y, x; F ](2) and m = 2, 0 < µ ≤ 4.720435 ℓ− 3.622464, ℓ ≥ 1, or m = 3, µ > 0, ℓ ≥ 1, (2) (2) we get CEI5 > CEI2 . We point out that Φ4 (frozen Steffensen) from a theoretical point of view presents the highest efficiency among the freederivative iterative methods considered in this work. Likewise, the results for Φ3 (frozen Newton) are stated in the following. Theorem 3.2. For ℓ ≥ 1 holds (a) CEI1 > CEI3 , if m µ1 >
1 2Cm2 + 6Aµm + 3(2A + C )ℓm + 9Cm + 6(C − B)µ + 9C ℓ − 11C ,
6A
using [y, x; F ](1) ,
and m µ1 >
1 2Cm2 + 12Aµm + 3(2A + C )ℓm + 9Cm − 12Bµ + 9C ℓ − 11C , 6A
using [y, x; F ](2) ,
where A = log ρ3 , B = log ρ1 and C = log (ρ3 /ρ1 ). (b) CEI2 > CEI3 , if m µ1 >
1 2Em2 + 6Aµm + 3(2A + E )ℓm + 9Em + 12E µ + 9E ℓ − 11E , 6D
using [y, x; F ](1) ,
and m µ1 >
1 2Em2 + 12Aµm + 3(2A + E )ℓm + 9Em + 6(E − D)µ + 9E ℓ − 11E ,
6D
using [y, x; F ](2) ,
where A = log ρ3 , D = log ρ2 and E = log (ρ3 /ρ2 ). (c) CEI4 > CEI3 for µ1 > µ + ℓ using [y, x; F ](1) and mµ1 > 2mµ + mℓ − µ using [y, x; F ](2) . (d) CEI5 > CEI3 for mµ1 > mµ + µ + mℓ + 1 using [y, x; F ](1) and mµ1 > 2mµ + mℓ + 1 using [y, x; F ](2) . The boundary between the regions of two computational efficiencies, Ri,j = 1, is a plane in case (c) with [y, x; F ](1) , a hyperbolic cylinder in case (d) with [y, x; F ](2) and a hyperbolic paraboloid (see [18]) in the other cases. 4. Numerical results In this section numerical experiments are performed to illustrate the results previously presented. We begin with some comparisons between the theoretical divided difference operator and the corresponding first computational operator given in (15). We observe that in some situations a losing of convergence order occurs. More precisely, we will see that iterative methods Φ2 and Φ5 have a lower order. Finally, two numerical examples are also given. 4.1. Theoretical considerations Using the Taylor expression of the theoretical divided difference operator [ x + h, x; F ] we have
[ x + h, x; F ]ij =
1
Dj Fi (x + th) dt
0
= Dj Fi (x) +
m 1
2 k=1
Dkj Fi (x) hk +
m 1
6 k,ℓ=1
Dk ℓj Fi (x) hk hℓ + O(h3 ).
Likewise, for the computational divided difference operator [ x + h, x; F ](1) we get
[ x + h, x; F ](ij1) = Dj Fi (x) +
j −1 k=1
Dkj Fi (x) hk +
1 2
Dj j Fi (x) hj + O(h2 ).
(18)
M. Grau-Sánchez et al. / Journal of Computational and Applied Mathematics 255 (2014) 753–764
From (18) we obtain [ x + h, x; F ](1) = F ′ (x) + M h + O(h2 ), instead of [ x + h, x; F ] = F ′ (x) +
M ̸=
1 2
F (x) because M depends on the function F (see [19] and the references therein).
1 2
759
F ′′ (x) + O(h2 ). In general,
′′
Setting Ψk = [ xn + hk , xn ; F ](1) , 1 ≤ k ≤ 4, we have
en−1 − en , where Ψ1 = [ xn−1 , xn ; F ](1) which is the operator of Φ1 defined in (8), 2(en−1 − en ), where Ψ2 = [ xn−1 , 2xn − xn−1 ; F ](1) which is given by Φ2 in (10), hk = −F (xn ) = −Γ (en + o(en )) = −e˜n + o(en ), where Φ3 = [ xn , xn + F (xn ); F ](1) , −2qn = −2(en−1 − en )T (en−1 − en ) = −2e2n−1 + o(e2n−1 ), if Φ4 = [ xn − qn , xn + qn ; F ](1) . Taking into account that F ′ (xn ) = Γ I + 2 A2 en + O(e2n ) and M = Γ (N + O(en )), from Ψk = Γ I + 2 A2 en + N hk + O(h2k ) , (k)
we get Ψk−1 = (I − 2 A2 en − N hk + o(en−1 )) Γ −1 . The error expression corresponding to the first step, zn − α , defined in (8), (10), (12) and (13) is
En(k) = zn(k) − α = en − Ψk−1 Γ en + O(e2n ) = (2 A2 en + N hk ) en + o(hk en ) = [ O(en−1 en ), O(en−1 en ), O(e2n ), O(e2n−1 en ) ],
for 1 ≤ k ≤ 4.
Moreover, the final error using the inverse operator of Ψk is
(k) en+1 = x(nk) − α = En(k) − Ψk−1 Γ En(k) + O((En(k) )2 ) = (2 A2 en + N hk ) En(k) + o(hk En(k) ) = [ O(e2n−1 en ), O(e2n−1 en ), O(e3n ), O(e4n−1 en ) ],
for 1 ≤ k ≤ 4.
(19)
From the results of (19) we may conclude methods Φ2 (k = 2) and Φ5 (k =4) for some set of functions the that iterative
√
17 /2 to ρ˜ 2 = 2 and from ρ5 = 3 to ρ˜ 5 = 1 +
local convergence order go from ρ2 = 1 +
√
17 /2. Furthermore, we
also observe that Φ1 (k = 1, frozen Secant) and Φ4 (k = 3, frozen Steffensen) keep the order independently of the function tested. Given a function F to obtain the maximum order of convergence it may be necessary to use the operator [y, x; F ](2) given in (16) for Φ2 and Φ5 . This requires a higher computational cost; while for Φ1 and Φ4 we can carry out [y, x; F ](1) defined (1) (2) (1) (2) in (15). Namely, the comparison of CEI1 , CEI2 , CEI4 and CEI5 is done in the following theorem. Theorem 4.1. For all µ > 0, m ≥ 2 and ℓ ≥ 1 we have (1)
CEI4 (2)
(a) CEI5
> CEI1(1) ,
(1)
CEI4
> CEI2(2) and CEI4(1) > CEI5(2) .
> CEI2(2) , (see Theorem 3.1(a) and (c)).
(b) For m ≥ 2, ℓ ≥ 1 and 0 < µ <
1 2m2 +9(1+ℓ)m+9ℓ−11 6 ( B −1)m−1 A
(c) For m ≥ 2, ℓ ≥ 1 and 0 < µ <
(2)
we obtain CEI2
D 2 1 2m +9(1+ℓ)m+9ℓ−5−6 C B 6 −1 (m+1)
> CEI1(1) , where A = log (ρ2 /ρ1 ) and B = log ρ1 . (2)
we obtain CEI5
> CEI1(1) , where C = log(ρ5 /ρ1 ) and D = log ρ5 .
C
For ℓ = 1.731, Fig. 1 illustrates the following four zones
(1): CEI4(1) > CEI1(1) > CEI2(2) > CEI5(2) (2): CEI4(1) > CEI1(1) > CEI5(2) > CEI2(2)
(20)
(3): CEI4(1) > CEI5(2) > CEI1(1) > CEI2(2) (4): CEI4(1) > CEI5(2) > CEI2(2) > CEI1(1) .
Finally, using the operator [y, x; F ](1) for all four methods and recalling that the convergence orders of Φ2 and Φ5 are √ ρ2 = 2 and ρ5 = 1 + 17 /2 respectively, we state the following theorem, where the CEIs with these two orders of (1)
(1)
2 and CEI 5 . local convergence are denoted by CEI Theorem 4.2. For all µ > 0, m ≥ 2 and ℓ ≥ 1 we have (1)
CEI4
(1)
2 > CEI1(1) > CEI
and
(1)
CEI4
(1)
(1)
5 > CEI 2 . > CEI
760
M. Grau-Sánchez et al. / Journal of Computational and Applied Mathematics 255 (2014) 753–764
m Fig. 1. Boundary curves Ri,j = 1 for ℓ = 1.731 and m = 2, 3, 4, 5 (see (20)).
(1)
(1)
5 > CEI for Finally, we obtain CEI 1 (a) (b) (c) (d)
all m ≥ 5 and ℓ ≥ 1; m = 2, ℓ ≥ 1 and 0 < µ < 1.729473 ℓ − 0.115673; m = 3, ℓ ≥ 1 and 0 < µ < 3.745439 ℓ + 1.7888879; m = 4, ℓ ≥ 1 and 0 < µ < 12.459539 ℓ + 11.128903.
The proof as before is left to the reader. 4.2. Numerical examples The numerical computations were performed using the MPFR library of C++ multi-precision arithmetics [20,21] with 4096 digits of mantissa. All algorithms were compiled by g++(4.2.1) for i686-apple-darwin1 with libgmp (v.5.0.2) R and libmpfr (v.3.1.0) libraries in a processor Intel⃝ Xeon E5620, 2.4 GHz (64-bit machine). In this machine the quotient and product ratio is ℓ = 1.731. In each example the starting point is the same for all methods tested. The classical stopping criteria ∥eI ∥ = ∥xI − α∥ > ε and ∥eI +1 ∥ < ε , with ε = 10−ν , where ν = 4096, are replaced by ∥˘eI ∥ > 10−η and ∥˘eI +1 ∥ < 10−η , where η = [ν(ρ − 1)/ρ ] and ∥ · ∥ = ∥ · ∥∞ . Moreover, e˘ n is obtained by
e˘ n =
F r ( xn )
Fr (xn−1 )
.
(21)
1≤r ≤m
Note that this criterion is independent of the knowledge of the root (see [22]). To evaluate the numerical efficiency of each method we compute its factor κ˜ (see [8]). That is, for each iteration k we get ˜ (Dk ), where the precision (number of correct decimals) Dk in time Θ Dk ≈ − log10 ∥ek ∥ ≈ −
ρ log10 ∥˘ek ∥. ρ−1
˜ (Dk )), 1 ≤ k ≤ I, in the least-squares sense by a polynomial of degree one, where the We approximate the pairs (log Dk , Θ slope κ˜ is computed. Namely, κ˜ is a coefficient measuring the time of execution in the function of the approximate number of correct decimals, say ˜ (DI ) = κ( Θ ˜ log DI − log D0 ). defined by Taking into account this slope we compare the time factor TF with the computed time factor TF = TF
˜ (DI ) Θ tp log q
=
κ˜ tp
≈ TF =
1 log CEI
,
where tp is the necessary time spent by one product and q = DI /D0 . In these experiments we calculate the computational order of convergence ρ˘ , called PCLOC, defined in [22] by
ρ˘ =
log ∥F (xI )∥ log ∥F (xI −1 )∥
.
M. Grau-Sánchez et al. / Journal of Computational and Applied Mathematics 255 (2014) 753–764
761
Table 3 Numerical results for the nonlinear system (22).
(1)
Φ1 (1) Φ2 Φ3 (1 ) Φ4 (1 ) Φ5
ρ
TF
I
T
qI
κ˜
TF
rTF
2.0000
3204.73
2.5616 3.0000
2922.52 1509.05
10
229.05
2897
80.31
3318.73
3.56
7 6
193.60 98.80
1631 1852
72.82 36.86
3009.11 1522.95
2.96 0.92
3.0000
2502.21
7
201.93
3212
62.70
2590.82
3.54
3.0000
2988.75
7
232.51
3949
74.32
3071.26
2.76
Table 4 Numerical results for the nonlinear system (23).
(1)
Φ1 (1) Φ2 (2 ) Φ2 Φ3 (1 ) Φ4 (1 ) Φ5 (2 ) Φ5
ρ
TF
I
T
qI
κ˜
TF
rTF
2.0000
4199.31
9
264.86
2113
101.77
4205.51
0.15
2.0000
5182.73
9
327.45
2737
116.56
4816.71
7.06
2.5616 3.0000
5326.98 2018.04
7 6
350.72 131.10
2454 2628
130.19 48.40
5379.92 2000.20
0.99 0.88 0.26
3.0000
3286.60
6
224.42
1740
79.75
3295.27
2.5616
4566.81
7
299.18
2257
11.82
4620.64
1.18
3.0000
5204.29
6
337.70
1950
126.86
5242.17
0.73
If ρ = ρ˘ ± ∆ρ˘ , where ρ is the local order of convergence and ∆ρ˘ is the error of PCLOC, then we get ∆ρ˘ < 10−3 . This fact means that in all computations of PCLOC we obtain at least 3 significant digits and this result is a good test of the local convergence orders of the family of iterative methods worked in this paper. Next we present two examples. The first one does not have any problem when we apply the definition [y, x; F ](1) given in (15). The second example needs the definition (16) for the Φ2 and Φ5 methods to attain the local convergence order that theoretical results state. Tables 3 and 4 show the results obtained for the iterative methods Φk , 1 ≤ k ≤ 5. In each table it is shown the local order of convergence ρ , the measure time factor TF defined as TF = 1/ log CEIk , the necessary iteration number I, the necessary time T in milliseconds to reach the I-iteration, the number of correct decimals in xI , say qI , the slope κ and the computed . Furthermore, the last column shows the percentage of relative error rTF between TF and TF . time factor TF In the first example it is considered a system F (x1 , x2 , x3 ) = 0 with three exponential equations defined by
x2 + x3 − e−x1 , F (x1 , x2 , x3 ) = x1 + x3 − e−x2 , x + x − e−x3 , 1 2
(22)
where (m, µ, µ1 = µ/3) = (3, 76.4, 25.5). Setting x−1 = (−0.83, 1.14, 1.14)t and x0 = (−0.84, 1.15, 1.15)t the following approximation of the root α ≈ (−.83202504 . . . , 1.1489838 . . . , 1.1489838 . . .)t with the prefixed precision is obtained. Notice that the preceding results agree with the theoretical ones exposed previously. Observing the computed time factor of the different methods we conclude that they also agree with the computational efficiency index CEI or the time factor TF TF . Moreover, the inequalities (1)
CEI3 > CEI4
> CEI2(1) > CEI5(1) > CEI1(1) ,
are in concordance with the ones in Theorems 3.1 and 3.2. In Fig. 2 the different slopes κ˜ for each method are drawn. Note that the minimum slope and time correspond to the Newton method Φ3 followed by the Steffensen method Φ4 . A second example where the theoretical local convergence orders are different from the ones studied in Section 2 and presented in Table 1 is x1 − cos(x1 − (x2 + x3 )), x2 − cos(x2 − (x1 + x3 )), x3 − cos(x3 − (x1 + x2 )),
F (x1 , x2 , x3 ) =
(23)
where (m, µ, µ1 ) = (3, 100.8, 34.5). Using the initial values x−1 = (0.544, 0.544, 0.996)t and x0 = (0.543, 0.543, 0.995)t we have the following approximation of the root
α ≈ (0.5438500415, 0.5438500415, 0.9957781534)t . Recall that, in this case, it is necessary to use the second computational definition of the divided difference operator. Table 4 shows the results of the system (23). Notice that there are more rows than in Table 3 because we have studied the iterative methods using the two divided difference operators for the iterative functions Φ2 and Φ5 . Also note that CEI , TF
762
M. Grau-Sánchez et al. / Journal of Computational and Applied Mathematics 255 (2014) 753–764
Fig. 2. Precisions and computational time of the five methods for example (22).
Fig. 3. Precisions and computational time of the five methods with maximum order for example (23).
are ordered in the same way as all the theorems stated in this paper. Furthermore, we have and TF (1)
CEI3 > CEI4
> CEI1(1) > CEI5(1) > CEI2(1) > CEI5(2) > CEI2(2) .
In Fig. 3 the different slopes κ˜ for each method are drawn. Observe that the minimum slope and time correspond to the (1) Newton method Φ3 , followed by the Steffensen method Φ4 . 5. Application Hereafter, we apply the technique presented previously to the well-known Hammerstein integral equation [23,24]. Namely, x(s) = 1 +
1 2
1
G(s, t ) x(t )3 dt ,
s ∈ [0, 1],
0
where x ∈ C [0, 1], t ∈ [0, 1], and the kernel G is the Green function G(s, t ) =
(1 − s) t , s (1 − t ),
t ≤ s, s < t.
M. Grau-Sánchez et al. / Journal of Computational and Applied Mathematics 255 (2014) 753–764
763
Table 5 Numerical results for the nonlinear system (25).
(1)
Φ1 (1) Φ2 Φ3 (1) Φ4 (1) Φ5
ρ
TF
I
2.0000
4089.12
2.5616 3.0000
3228.73 1241.12
3.0000 3.0000
T
qI
10
192.85
2866
8 7
152.03 63.01
4035 3795
2764.39
6
124.16
1795
2965.60
6
123.23
1402
We begin by writing the preceding integral equation as R(x) = 0, where R : C [0, 1] → C [0, 1] and
[R(x)](s) = x(s) − 1 −
1 2
1
G(s, t ) x(t )3 dt ,
s ∈ [0, 1].
(24)
0
To obtain the discretization of the integral we can use Gauss–Legendre formulae. That is, 1
f (t ) dt ≈ 0
m
ϖj f (tj ),
j =1
where the abscissas tj and the weights ϖj were computed with maximum precision for m = 8 in [25]. Next, we present the discrete version of (24). Indeed, denoting the approximation of x(ti ) by xi (1 ≤ i ≤ 8), we get the following system of nonlinear equations xi = 1 +
8 1
2 j=1
bij x3j ,
where bij =
ϖj tj (1 − ti ) ϖj ti (1 − tj )
if j ≤ i, if j > i,
1 ≤ i ≤ 8,
(25)
with solution
α ≈ (1.005771, 1.027442, 1.055178, 1.074414, 1.074414, 1.055178, 1.027442, 1.005771)t .
(26)
In this case the values of parameters are (m, µ, µ1 ) = (8, 11, 11/8). Taking x−1 = (1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0)t and x0 = (1.1, 1.1, 1.1, 1.1, 1.1, 1.1, 1.1, 1.1)t we show and com(1) pare the convergence of the methods Φ k , 1 ≤ k ≤ 5, towards the root α (26) in Table 5. 6. Concluding remarks In this paper we present a generalized family of two step Secant-like methods with frozen operator for solving nonlinear equations. A local convergence analysis of this family, where five specific cases are considered, is studied in detail. Namely, using earlier methods such as Secant, Kurchatov, Newton, Steffensen and other new algorithms, the family of iterative schemes is built up and their computational efficiency is also given. Numerical examples using multiple precision and a stopping criterion are implemented without using any known root. Finally, a study comparing the order, efficiency and elapsed time of the methods suggested supports the theoretical results claimed. We point out that the real efficiency of each scheme is related to the slope of the lines shown in the figures presented. Acknowledgement The research was supported partially by the project MTM2011-28636-C02-01. References [1] J.A. Ezquerro, M.A. Hernández, An optimization of Chebyshev’s method, J. Complexity 25 (2009) 343–361. [2] M. Grau, J.L. Díaz-Barrero, An improvement of the Euler–Chebyshev iterative method, J. Math. Anal. Appl. 315 (2006) 1–7. [3] M. Grau-Sánchez, J.M. Gutiérrez, Some variants of the Chebyshev–Halley family of methods with fifth order of convergence, Int. J. Comput. Math. 87 (2010) 818–833. [4] J.A. Ezquerro, M.A. Hernández, M.J. Rubio, Secant-like methods for solving nonlinear integral equations of the Hammerstein type, J. Comput. Appl. Math. 115 (1–2) (2000) 245–254. [5] I.K. Argyros, J.A. Ezquerro, J.M. Gutiérrez, M.A. Hernández, S. Hilout, On the semilocal convergence of efficient Chebyshev–Secant-type methods, J. Comput. Appl. Math. 235 (2011) 3195–3206. [6] I.K. Argyros, H. Ren, Efficient Steffensen-type algorithms for solving nonlinear equations, Int. J. Comput. Math. 90 (2013) 691–704. [7] M. Grau-Sánchez, A. Grau, M. Noguera, Frozen divided difference scheme for solving systems of nonlinear equations, J. Comput. Appl. Math. 235 (2011) 1739–1743. [8] M. Grau-Sánchez, M. Noguera, A technique to choose the most efficient method between secant method and some variants, Appl. Math. Comput. 218 (2012) 6415–6426.
764
M. Grau-Sánchez et al. / Journal of Computational and Applied Mathematics 255 (2014) 753–764
[9] J.A. Ezquerro, A. Grau, M. Grau-Sánchez, M.A. Hernández, M. Noguera, Analysing the efficiency of some modifications of the secant method, Comput. Math. Appl. 64 (2012) 2066–2073. [10] J.A. Ezquerro, M. Grau-Sánchez, M.A. Hernández, Solving non-differentiable equations by a new one-point iterative method with memory, J. Complexity 28 (2012) 48–58. [11] J.A. Ezquerro, A. Grau, M. Grau-Sánchez, M.A. Hernández, Construction of derivative-free iterative methods from Chebyshev’s method, Anal. Appl. 11 (2013) 1350009 (16 pp.). [12] V.A. Kurchatov, On a method of linear interpolation for the solution of functional equations, Dokl. Akad. Nauk SSSR 198 (3) (1971) 524–526 (in Russian); translation in Sov. Math. Dokl. 12 (1971) 835–838. [13] J.A. Ezquerro, M. Grau-Sánchez, M.A. Hernández, M. Noguera, Semilocal convergence of secant-like methods for differentiable and nondifferentiable operator equations, J. Math. Anal. Appl. 398 (2013) 100–112. [14] S. Amat, S. Busquier, J.M. Gutiérrez, Third-order iterative methods with applications to Hammerstein equations: a unified approach, J. Comput. Appl. Math. 235 (2011) 2936–2943. [15] I.K. Argyros, On the Newton–Kantorovich hypothesis for solving equations, J. Comput. Appl. Math. 169 (2004) 315–332. [16] S. Amat, S. Busquier, V. Candela, A class of quasi-Newton generalized Steffensen methods on Banach spaces, J. Comput. Appl. Math. 149 (2002) 397–406. [17] V. Alarcón, S. Amat, S. Busquier, D.J. López, A Steffensen’s type method in Banach spaces with applications on boundary-value problems, J. Comput. Appl. Math. 216 (2008) 243–250. [18] M. Grau-Sánchez, A. Grau, M. Noguera, On the computational efficiency index and some iterative methods for solving systems of nonlinear equations, J. Comput. Appl. Math. 236 (2011) 1259–1266. [19] M. Grau-Sánchez, M. Noguera, S. Amat, On the approximation of derivatives using divided difference operators preserving the local convergence order of iterative methods, J. Comput. Appl. Math. 237 (2013) 363–372. [20] L. Fousse, G. Hanrot, V. Lefèvre, P. Pélissier, P. Zimmermann, MPFR: a multiple-precision binary floating-point library with correct rounding, ACM Trans. Math. Software 33 (2) (2007) 15. Art. 13. [21] The GNU MPFR library 3.1.0. Available in http://www.mpfr.org. [22] M. Grau-Sánchez, A. Grau, M. Noguera, J.R. Herrero, A study on new computational local orders of convergence, Appl. Math. Lett. 25 (2012) 2023–2030. [23] J.M. Ortega, W.C. Rheinboldt, Iterative Solution of Nonlinear Equations in Several Variables, Academic Press, New York, 1970. [24] A.D. Polyanin, A.V. Manzhirov, Handbook of Integral Equations, CRC Press, Boca Raton, 1998. [25] J.A. Ezquerro, M. Grau-Sánchez, A. Grau, M.A. Hernández, M. Noguera, N. Romero, On iterative methods with accelerated convergence for solving systems of nonlinear equations, J. Optim. Theory Appl. 151 (2011) 163–174.