A novel family of composite Newton–Traub methods for solving systems of nonlinear equations

Applied Mathematics and Computation 269 (2015) 520–535

Contents lists available at ScienceDirect

Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc

A novel family of composite Newton–Traub methods for solving systems of nonlinear equations Janak Raj Sharma a,∗, Rajni Sharma b, Nitin Kalra b,c a

Department of Mathematics, Sant Longowal Institute of Engineering and Technology, Longowal 148106, Punjab, India Department of Applied Sciences, D.A.V. Institute of Engineering and Technology, Jalandhar 144008, Punjab, India c Research scholar, I.K. Gujral Punjab Technical University, Jalandhar, Punjab, India b

a r t i c l e

i n f o

a b s t r a c t

Keywords: Systems of nonlinear equations Newton’s method Traub’s method Order of convergence Computational efficiency

We present a family of three-step iterative methods of convergence order five for solving systems of nonlinear equations. The methodology is based on the two-step Traub’s method with cubic convergence for solving scalar equations. Computational efficiency of the new methods is considered and compared with some well-known existing methods. Numerical tests are performed on some problems of different nature, which confirm robust and efficient convergence behavior of the proposed methods. Moreover, theoretical results concerning order of convergence and computational efficiency are verified in the numerical problems. Stability of the methods is tested by drawing basins of attraction in a two-dimensional polynomial system. © 2015 Elsevier Inc. All rights reserved.

1. Introduction The construction of fixed point methods for solving nonlinear equations and systems of nonlinear equations is an interesting and challenging task in numerical analysis and many applied scientific branches. A great importance of this topic has led to the development of many numerical methods, most frequently of iterative nature (see [1–5]). With the advancement of computer hardware and software, the problem of solving nonlinear equations by numerical methods has gained an additional importance. In this paper, we consider the problem of finding solution of the system of nonlinear equations F(x) = 0 by iterative methods of a high order of convergence. This problem can be precisely stated as to find a vector r = (r1 , r2 , . . . , rn )T such that F(r) = 0, where F : D ⊆ Rn → Rn is the given nonlinear vector function F(x) = ( f1 (x), f2 (x), . . . , fn (x))T and x = (x1 , x2 , . . . , xn )T . The solution vector r of F(x) = 0 can be obtained as a fixed point of some function φ : Rn → Rn by means of the fixed point iteration

x(k+1) = φ(x(k) ),

k = 0, 1, 2, . . . .

One of the basic procedures for solving systems of nonlinear equations is the quadratically convergent Newton’s method (see [1–3]), which is given as,

x(k+1) = φ1(2) (x(k) ) = x(k) − F (x(k) )−1 F(x(k) ),

(1)

where F (x)−1 is the inverse of first Fréchet derivative F (x) of the function F(x). In terms of computational cost Newton’s method requires the evaluations of one F, one F and one matrix inversion per iteration. Throughout the paper, we shall use the abbrevi( p)

ation φi ∗

to denote an ith iterative function of convergence order p.

Corresponding author. Tel.: +91 1672 253256; fax: +91 1672 280057. E-mail addresses: [email protected] (J.R. Sharma), [email protected] (R. Sharma), [email protected] (N. Kalra).

http://dx.doi.org/10.1016/j.amc.2015.07.092 0096-3003/© 2015 Elsevier Inc. All rights reserved.

J.R. Sharma et al. / Applied Mathematics and Computation 269 (2015) 520–535

521

In order to improve the order of convergence of Newton’s method, many higher order methods have been proposed in literature. For example, Frontini and Sormani [6], Homeier [7], Cordero and Torregrosa [8], Noor and Waseem [9], and Grau et al. [10] have developed third order methods requiring one F, two F and two matrix inversions per iteration. Cordero and Torregrosa have also derived two third-order methods in [11]. One of the methods require one F and three F whereas other requires one F and four F evaluations per iteration. Both the methods also require two matrix inversions in each iteration. Darvishi and Barati [12] have proposed a third order method which uses two F, one F and one matrix inversion. Grau et al. presented a fourth order method in [10] utilizing three F, one F and one matrix inversion. Cordero et al. developed a fourth order method in [13], which uses two F, two F and one matrix inversion. Darvishi and Barati [14] presented a fourth order method requiring two F, three F and two matrix inversions per iteration. Cordero et al. in [15] have implemented fourth order Jarratt’s method [16] for scalar equations to systems of equations which requires one F, two F and two matrix inversions. Sharma et al. [17] developed a fourth order method requiring one F, two F and two matrix inversions. In quest of more fast algorithms, researchers have also proposed fifth and sixth order methods in [10,13,15,18,19]. The fifth order methods by Grau et al. [10] and Cordero et al. [15,18] require four evaluations namely, two F and two F per iteration. The fifth order method by Cordero et al. [13] requires three F and two F . In addition, the fifth order method in [10] requires two matrix inversions, in [13] one matrix inversion and in [15,18] three matrix inversions. One sixth order method by Cordero et al. [15] uses two F and two F while other sixth order method [18] uses three F and two F . The sixth order methods, apart from the mentioned evaluations, also require two matrix inversions per one iteration. Recently, Sharma and Gupta [19] proposed a fifth order method requiring two F, two F and two matrix inversions per one iteration. The main goal of this paper is to develop iterative methods of high computational efficiency, which may assume a high convergence order and low computational cost. To do so, we here propose a one-parameter family of methods with fifth order of convergence by employing the iterative scheme that utilizes the number of function evaluations and inverse operators as minimum as possible. In this way, we attain low computational cost and hence an increased computational efficiency. Moreover, we show that the proposed methods are efficient than existing methods of similar nature in general. The scheme of present contribution consists of three steps of which the first two steps are the generalizations of Traub’s third order two-step scheme ([1], p. 181) for solving scalar equation f (x) = 0, which is given as

f (x ) y k = xk − θ  k , f (xk ) xk+1 = xk −



1 1+ 2θ



1 f  (yk ) − 2θ f  (xk )



f (xk ) , f  (xk )

θ ∈ R \ {0 },

(2)

whereas the third step is based on Newton-like scheme. In terms of computational cost the proposed technique utilizes two F, two F and only one matrix inversion per iteration. Our presentation is organized as follows. Some basic definitions relevant to the present work are given in Section 2. In Section 3, the fifth order scheme is developed and its behavior is analyzed. Computational efficiency of the new methods is studied and then compared with some well-known existing methods in Section 4. In Section 5, we present various numerical examples to confirm the theoretical results and to compare convergence properties of the proposed methods with existing methods. In Section 6, we study dynamics of the methods on a two-dimensional polynomial system. Concluding remarks are given in Section 7. 2. Basic definitions 2.1. Order of convergence Let {x(k) }k ≥ 0 be a sequence in Rn which converges to r. Then, convergence is called of order p, p > 1, if there exists M, M > 0, and k0 such that

||x(k+1) − r||  M||x(k) − r|| p ∀ k  k0 or

||e(k+1) ||  M||e(k) || p ∀ k  k0 , where e(k) = x(k) − r. The convergence is called linear if p = 1 and there exists M such that 0 < M < 1. 2.2. Error equation Let e(k) = x(k) − r be the error in the kth iteration, we call the relation

e(k+1) = L(e(k+1) ) p + O((e(k) ) p+1 ), p

as the error equation. Here, p is the order of convergence, L is a p-linear function, i.e. L ∈ L(Rn × · · · · ×Rn , Rn ), L denotes the set of bounded linear functions.

522

J.R. Sharma et al. / Applied Mathematics and Computation 269 (2015) 520–535

2.3. Computational order of convergence Let x(k−1) , x(k) and x(k+1) be the three consecutive iterations close to the zero r of F(x). Then, the computational order of convergence can be approximated using the formula (see [20,21])

ρk =

log (||F(x(k) )||/||F(x(k−1) )||) . log (||F(x(k−1) )||/||F(x(k−2) )||)

2.4. Computational efficiency Computational efficiency of an iterative is measured by the efficiency index E = p1/C (see [10]), where p is the order of convergence and C is the computational cost, which is given by

C (μ0 , μ1 , n) = P0 (n)μ0 + P1 (n)μ1 + P (n). Here, P0 (n) represents the number of evaluations of scalar functions ( f1 , f2 , . . . , fn ) used in the evaluation of F, P1 (n) is the ∂f number of evaluations of scalar functions of F , i.e. ∂ x i , 1 ≤ i, j ≤ n, P(n) represents the number of products or quotients needed j

per iteration, and μ0 and μ1 are the ratios between products and evaluations required to express the value of C(μ0 , μ1 , n) in terms of products. 3. Development of the method We consider the Traub’s scheme (2) for solving the nonlinear system F(x) = 0, thereby writing it in the generalized form as

y(k) = x(k) − θ F (x(k) )−1 F(x(k) ), z(k) = x(k) −



1+





1 1  (k) −1  (k) I− F (x ) F (y ) F (x(k) )−1 F(x(k) ), 2θ 2θ

(3a)

where θ ∈ R \ {0}. Now based on the two-step scheme (3a), let us write third step to obtain the approximation x(k+1) to a solution of F(x) = 0 in the following way:

x(k+1) = z(k) −



α I + β F (x(k) )−1 F (y(k) ) F (x(k) )−1 F(z(k) ),

(3b)

where α and β are some parameters and I denotes an n × n identity matrix. It is clear that the third step is a Newton-like step. For this reason we call the scheme (3b) along with (3a) as the composite Newton–Traub scheme. In order to explore the convergence property of this three-step scheme, we need the following result of Taylor’s expansion on vector functions (see [2]). Lemma 1. Let F : D ⊆ Rn → Rn be a q-times Fréchet differentiable in a convex set D ⊆ Rn , then for any x, h ∈ Rn the following expression holds:

F(x + h) = F(x) + F (x)h + where ||Rq || 

1 q!

1  1  1 F (x)h2 + F (x)h3 + · · · + F(q−1) (x)hq−1 + Rq , 2! 3! q!

(4)

sup ||F(q) (x + th)||||h||q and hq = (h, h, . .q. ., h).

0
We are in a position to prove the following theorem: Theorem 1. Let the function F : D ⊆ Rn → Rn be sufficiently differentiable in a convex set D containing the zero r of F(x). Suppose that F (x) is continuous and nonsingular in r. If an initial approximation x(0) is sufficiently close to r, then the local convergence order of the composite Newton–Traub method is at least five, provided α = 1 + θ1 and β = − θ1 . Proof. Let e(k) = x(k) − r. Developing F(x(k) ) in the neighborhood of r and assuming that  = F (r)−1 exists, we write





F(x(k) ) = F (r) e(k) + A2 (e(k) )2 + A3 (e(k) )3 + A4 (e(k) )4 + A5 (e(k) )5 + O((e(k) )6 ) , where Ai = Also,

(5)

1  F(i) (r), F(i) (r) ∈ L(Rn ×, . .i . ., ×Rn , Rn ),  ∈ L(Rn , Rn ) and (e(k) )i = (e(k) ,e(k) , . .i . ., e(k) ) with e(k) ∈ Rn . i!





F (x(k) ) = F (r) I + 2A2 (e(k) ) + 3A3 (e(k) )2 + 4A4 (e(k) )3 + O((e(k) )4 ) . Inversion of F (x(k) ) yields,





F (x(k) )−1 = I + X1 (e(k) ) + X2 (e(k) )2 + X3 (e(k) )3 + O((e(k) )4 )  , where

X1 = −2A2 , X2 = 4A22 − 3A3 , X3 = −(8A32 − 6A2 A3 − 6A3 A2 + 4A4 ).

(6)

(7)

J.R. Sharma et al. / Applied Mathematics and Computation 269 (2015) 520–535

523

Post multiplication of the above equation by F(x(k) ) yields,

F (x(k) )−1 F(x(k) ) = e(k) − A2 (e(k) )2 + 2(A22 − A3 )(e(k) )3 + ( − 4A32 + 4A2 A3 + 3A3 A2 − 3A4 )(e(k) )4 + O((e(k) )5 ).

(8)

Using (8) in the first step of (3a), we obtain

e˜(k) = y(k) − r = (1 − θ )e(k) + θ A2 (e(k) )2 − 2θ (A22 − A3 )(e(k) )3 + θ (4A32 − 4A2 A3 − 3A3 A2 + 3A4 )(e(k) )4 + O((e(k) )5 ). Expanding F (y(k) ) about r and using above result, we have



F (y(k) ) = F (r) I + 2A2 (e˜(k) ) + 3A3 (e˜(k) )2 + 4A4 (e˜(k) )3 + O(e˜(k) )4

(9)





= F (r) I + 2(1 − θ )A2 e(k) + (2θ A22 + 3(1 − θ )2 A3 )(e(k) )2 +

( − 4θ (A32 − A2 A3 ) + 6θ (1 − θ )A3 A2 + 4(1 − θ )3 A4 )(e(k) )3 + O((e(k) )4 ) .

(10)

Then, from (7) and (10), it follows that

F (x(k) )−1 F (y(k) ) = I − 2θ A2 e(k) + θ (6A22 − 3(2 − θ )A3 )(e(k) )2 + θ ((16 − 6θ )A2 A3 + 6(2 − θ )A3 A2 − 16A32 − 4(3 − 3θ + θ

2

)A4 )(e(k) )3 + O((e(k) )4 ).

(11)

Using Eqs. (8) and (11) in the second step of (3a), we obtain



eˆ(k) = z(k) − r =



2A22 +

 

  9θ A3 A2 (e(k) )3 − 9A32 + ( − 6 + 3θ )A2 A3 + −6 + 2  +(3 − 6θ + 2θ 2 )A4 (e(k) )4 + O((e(k) )5 ).

3θ − 1 A3 2

(12)

Simple algebraic calculations by using (11) yield

α I + β F (x(k) )−1 F (y(k) ) = (α + β)I − 2βθ A2 e(k) + βθ (6A22 − 3(2 − θ )A3 )(e(k) )2 + βθ ((16 − 6θ )A2 A3 + 6(2 − θ )A3 A2 − 16A32 − 4(3 − 3θ + θ Using Taylor’s series of

F(z(k) )

2

)A4 )(e(k) )3 + O((e(k) )4 ).

(13)

about r, we have

  F(z(k) ) = F (r) eˆ(k) + O((eˆ(k) )2 ) .

(14)

Employing Eqs. (5), (13) and (14) in (3b) and simplifying, we get

e(k+1) = (1 − α − β)eˆ(k) + 2(α + β + βθ )A2 e(k) eˆ(k) − ((α + β)(4A22 − 3A3 ) + 4βθ A22 + β(6θ A22 + 3(θ 2 − 2θ )A3 ))(e(k) )2 eˆ(k) + O((e(k) )6 ).

(15)

Our aim is to find the values of parameters α and β in such a way that the proposed iterative scheme may produce order of convergence as high as possible. Thus, if we take α + β = 1 and β = − θ1 , the above equation yields



e(k+1) = − 6A22 + 3(1 − θ )A3 (e(k) )2 eˆ(k) + O((e(k) )6 ). Combining (12) and (16),

e(k+1) =



−

6A22



+ 3(1 − θ )A3

2A22

+

 3θ 2

  − 1 A3 (e(k) )5 + O((e(k) )6 ).

(16)

(17)

This shows the fifth order convergence. Hence, the proof is completed.  Finally, we present the proposed composite Newton–Traub scheme in the form

y(k) = x(k) − θ F (x(k) )−1 F(x(k) ),







1 1  (k) −1  (k) I− F (x ) F (y ) F (x(k) )−1 F(x(k) ), 2θ 2θ  1   1 x(k+1) = z(k) − 1 + I − F (x(k) )−1 F (y(k) ) F (x(k) )−1 F(z(k) ).

z(k) = x(k) −

1+

θ

θ

(18)

The scheme (18) defines a one-parameter family of three-step fifth order methods with the first two steps as that of Traub’s method (3a) and third is Newton-like step. It is clear that this scheme requires the evaluations of two functions, two first derivatives and only one matrix inversion per one iteration. Some special cases of the family of composite Newton–Traub methods (18) are:

524

J.R. Sharma et al. / Applied Mathematics and Computation 269 (2015) 520–535

(5)

Case 1. For θ = 1, the special case of (18), that we denote by φ1 , is given as

y(k) = x(k) − F (x(k) )−1 F(x(k) ),

1 3I − F (x(k) )−1 F (y(k) ) F (x(k) )−1 F(x(k) ), 2  x(k+1) = z(k) − 2I − F (x(k) )−1 F (y(k) ) F (x(k) )−1 F(z(k) ).

z(k) = x(k) −

(19)

(5)

Case 2. For θ = −1, the special case (denoted by φ2 ) is defined as

y(k) = x(k) + F (x(k) )−1 F(x(k) ),

1 I + F (x(k) )−1 F (y(k) ) F (x(k) )−1 F(x(k) ), 2 x(k+1) = z(k) − F (x(k) )−1 F (y(k) )F (x(k) )−1 F(z(k) ). z(k) = x(k) −

(20)

1 (5) Case 3. For θ = − , we obtain the special case (denoted by φ3 ) given as 2

1  (k) −1 F (x ) F(x(k) ), 2 z(k) = x(k) − F (x(k) )−1 F (y(k) )F (x(k) )−1 F(x(k) ),

y(k) = x(k) +





x(k+1) = z(k) + I − 2F (x(k) )−1 F (y(k) ) F (x(k) )−1 F(z(k) ).

(21)

4. Computational efficiency In order to assess the computational efficiency of presented methods, we will consider all possible number of evaluations that contribute to the total cost of computation. For example, to compute F in any iterative method we evaluate n scalar functions, whereas the number of scalar evaluations is n2 for any new derivative F . In addition, we must include the amount of computational work required to evaluate the inverse of a matrix. Instead of computing the inverse operator we solve a linear system, where we have n(n − 1)(2n − 1)/6 products and n(n − 1)/2 quotients in the LU decomposition, and n(n − 1) products and n quotients in the resolution of two triangular linear systems. Moreover, we must add n2 products for the multiplication of a matrix with a vector or of a matrix by a scalar and n products for the multiplication of a vector by a scalar. We suppose that a quotient is equivalent to l products. (5) Computational efficiency of the presented Newton–Traub methods φi , (i = 1, 2, 3) is compared with existing fourth order method by Cordero et al. [13], fourth order generalized Jarratt’s method [16], fifth order method by Cordero et al. [18] and fifth order Sharma–Gupta method [19]. In addition, the presented methods are also compared with each other. The existing mentioned methods are given as follows: (4) Method by Cordero et al. (φ1 ):

y(k) = φ1(2) (x(k) ),





x(k+1) = y(k) − 2F (x(k) )−1 − F (x(k) )−1 F (y(k) )F (x(k) )−1 F(y(k) ). (4)

Generalized Jarratt method (φ2 ):

2  (k) −1 F (x ) F(x(k) ), 3  −1   1 3F (y(k) ) − F (x(k) ) 3F (y(k) ) + F (x(k) ) F (x(k) )−1 F(x(k) ). x(k+1) = x(k) − 2

y(k) = x(k) −

(5)

Method by Cordero et al. (φ4 ):

y(k) = φ1(2) (x(k) ),



−1

z(k) = x(k) − 2 F (y(k) ) + F (x(k) ) x(k+1) = z(k) − F (y(k) )−1 F(z(k) ). (5)

Sharma–Gupta method (φ5 ):

1  (k) −1 F (x ) F(x(k) ), 2 z(k) = x(k) − F (y(k) )−1 F(x(k) ),

y(k) = x(k) −

F(x(k) ),

J.R. Sharma et al. / Applied Mathematics and Computation 269 (2015) 520–535

525

x(k+1) = z(k) − [2F (y(k) )−1 − F (x(k) )−1 ]F(z(k) ). ( p)

( p)

Let us denote efficiency indices of φi by Ei account all possible evaluations, we have

C1(4) = 2nμ0 + 2n2 μ1 + C2(4) = nμ0 + 2n2 μ1 +

( p)

and computational cost by Ci

n 2 2n + 21n − 11 + 3l (n + 5) 6

n 2 2n + 9n + 1 + 3l (n + 1) 3

and

and

. Then using the definition 2.4 while taking into (4)

E1(4) = 41/C1 .

(22)

(4)

E2(4) = 41/C2 .

(23)

C1(5) = 2nμ0 + 2n2 μ1 +

n 2 2n + 33n − 5 + 3l (n + 7) 6

C2(5) = 2nμ0 + 2n2 μ1 +

n 2 2n + 33n − 17 + 3l (n + 7) 6

and E2(5) = 51/C2 .

C3(5) = 2nμ0 + 2n2 μ1 +

n 2 2n + 33n − 11 + 3l (n + 7) 6

and E3(5) = 51/C3 .

C4(5) = 2nμ0 + 2n2 μ1 +

n 2 2n + 3n − 3 + 3l (n + 1) 2

and E4(5) = 51/C4 .

C5(5) = 2nμ0 + 2n2 μ1 +

n 2 2n + 9n − 5 + 3l (n + 3) 3

and E5(5) = 51/C5 .

and

(5)

E1(5) = 51/C1 .

(24)

(5)

(25)

(5)

(26)

(5)

(27)

(5)

(28)

4.1. Efficiency comparison ( p)

To compare the computational efficiencies of the iterative methods, say φi

(q)

and φ j , we consider the ratio

C (j q) log ( p) = , (q) Ci( p) log (q) log E j ( p)

log E i

Ri,p,qj =

(29) ( p)

It is clear that if Ri, j > 1, the iterative method φi p,q

two computational efficiencies is given by

p,q Ri, j

(q)

is more efficient than φ j . Taking into account that the boundary between

= 1, this boundary is given by the equation μ0 written as a function of μ1 , n and

l; (μ1 , μ0 ) ∈ (0, +∞) × (0, +∞), n is a positive integer ≥ 2 and l ≥ 1. In the sequel, we have r = log (5) and t = log (2).

φ1(5) versus φ1(4) case:

For this case the boundary R5,4 = 1, expressed in μ0 as function of μ1 , n and l, is given by 1,1

μ0 =

12(2t − r)nμ1 + 2(2t − r)n2 + 3(22t − 7r)n + (11r − 10t ) + 3(2t − r)ln + (42t − 15r)l . 12(r − 2t ) (5)

(4)

Here, μ0 is always negative for n ≥ 27, so the boundary is out of the admissible region and we conclude that E1 > E1 for μ1 > 0, n ≥ 27 and l ≥ 1. For 2 ≤ n ≤ 26, we present some boundary lines Li (i=1, 2, 3, 4) in the (μ1 , μ0 )-plane corresponding to n = 2, 5, 10 and 15 taking l = 2.66 in each case. Reason for selecting the value 2.66 for l will be clear in the next section. These (5) (4) boundaries are the straight lines with negative slopes, where E1 > E1 on the right and E14 > E15 on the left of each line (see Fig. 1).

φ1(5) versus φ2(4) case:

In this case the boundary R5,4 = 1 is expressed by 1,2

μ0 =

3(2r − 4t )nμ1 + 2(r − t )n2 + 3(3r − 11t )n + (r + 5t ) + 3(r − t )ln + 3(r − 7t )l . 3(4t − r) (4)

In order to compare the efficiencies E2

(5)

and E1 , we present some boundaries Li (i=1, 2, 3, 4) corresponding to n = 2, 5, 10 and (5)

15 (taking l = 2.66) in the (μ1 , μ0 )-plane. Such boundaries are the straight lines with positive slopes, where E1 (4)

below and E2

(5)

> E1

(5)

For this case it is sufficient to compare the corresponding values of C15 and C25 in (24) and (25). Let C1 that E2

(5)

> E1

on the

on the above of each line (see Fig. 2).

φ1(5) versus φ2(5) case: (5)

(4)

> E2

(5)

> C2

which implies

for all l ≥ 1 and n ≥ 2.

φ1(5) versus φ3(5) case:

(5)

In this case comparison of the corresponding values of C15 and C35 in (24) and (26) gives E3

(5)

> E1

for all l ≥ 1 and n ≥ 2.

526

J.R. Sharma et al. / Applied Mathematics and Computation 269 (2015) 520–535

20 L1 L2

15

µ0

L3 L4

10

5

0

0

5

10

15

20

µ1 Fig. 1. Boundary lines in the (μ1 , μ0 )-plane for E14 and E15 .

200

µ0

150

100

L1 L2 L3

50

L4

0

0

50

100

150

200

µ1 Fig. 2. Boundary lines in the (μ1 , μ0 )-plane for E24 and E15 .

φ1(5) versus φ4(5) case:

(5)

Similarly, here also we compare the corresponding values of C15 and C45 in (24) and (27). Thus we find that E1 l ≥ 1 and n ≥ 5.

φ1(5) versus φ5(5) case:

(5)

Comparison of the corresponding values of C15 and C55 in (24) and (28) gives E1 (5)

(4)

φ2 versus φ1 case:

(5)

> E5

for all l ≥ 1 and n ≥ 5.

(5)

> E4

for all

J.R. Sharma et al. / Applied Mathematics and Computation 269 (2015) 520–535

527

20 L1 L2

15

µ0

L3 L4

10

5

0

0

5

10

15

20

µ1 Fig. 3. Boundary lines in the (μ1 , μ0 )-plane for E14 and E25 .

The boundary R5,4 = 1 is given by 2,1

μ0 =

12(2t − r)nμ1 + 2(2t − r)n2 + 3(22t − 7r)n + (11r − 34t ) + 3(2t − r)ln + (42t − 15r)l . 12(r − 2t )

It is easy to check that μ0 is always negative for n ≥ 26, so the boundary is out of the admissible region which implies that (5) (4) E2 > E1 for μ1 > 0, n ≥ 26 and l ≥ 1. Otherwise, for 2 ≤ n ≤ 25 the comparison depends on μ0 , μ1 and l. We present some boundaries Li (i=1, 2, 3, 4) in the (μ1 , μ0 )-plane corresponding to n = 2, 5, 10 and 15 taking l = 2.66. These boundaries are the (5) (4) (4) (5) straight lines with negative slopes, where E2 > E1 on the right and E1 > E2 on the left of each line (see Fig. 3).

φ2(5) versus φ2(4) case:

For this case the boundary R5,4 = 1 is given as 2,2

μ0 =

3(2r − 4t )nμ1 + 2(r − t )n2 + 3(3r − 11t )n + (r + 17t ) + 3(r − t )ln + 3(r − 7t )l . 3(4t − r)

We draw some particular boundaries Li (i=1, 2, 3, 4) in the (μ1 , μ0 )-plane using the same set of values of n and l as in the previous (5) (4) (5) (4) case for comparing the efficiencies E2 and E2 . These boundaries are the straight lines with positive slopes, where E2 > E2 (4)

on the below and E2

(5)

> E2

on the above of each line (see Fig. 4).

φ2(5) versus φ3(5) case:

(5)

In this case, it can be easily seen that comparison of the corresponding values of C2 for all l ≥ 1 and n ≥ 2.

φ2(5) versus φ4(5) case:

(5)

Comparing the corresponding values of C2 (5)

(5)

and C4

(5)

in (25) and (27), we have E2

(5)

and C3

(5)

> E4

(5)

in (25) and (26) yields E2

(5)

> E3

for all l ≥ 1 and n ≥ 4.

(5)

φ2 versus φ5 case:

(5)

Here also we compare the corresponding values of C2 and n ≥ 4.

(5)

and C5

φ3(5) versus φ1(4) case:

The boundary R5,4 = 1, expressed in μ0 as function of μ1 , n and l, is given by 3,1

μ0 =

(5)

in (25) and (28) which implies that E2

12(2t − r)nμ1 + 2(2t − r)n2 + 3(22t − 7r)n + 11(r − 2t ) + 3(2t − r)ln + (42t − 15r)l . 12(r − 2t )

(5)

> E5

for all l ≥ 1

528

J.R. Sharma et al. / Applied Mathematics and Computation 269 (2015) 520–535

200

µ0

150

100

L1 L2 L3

50

L4

0

0

50

100

150

200

µ1 Fig. 4. Boundary lines in the (μ1 , μ0 )-plane for E24 and E25 .

20 L1 L2

15

µ0

L3 L4

10

5

0

0

5

10

15

20

µ1 Fig. 5. Boundary lines in the (μ1 , μ0 )-plane for E14 and E35 .

(5)

(4)

Note that μ0 is always negative for n ≥ 26, so the boundary is out of the admissible region which shows that E3 > E1 for μ1 > 0, n ≥ 26 and l ≥ 1. On the other hand, for 2 ≤ n ≤ 25 the efficiency comparison depends on μ0 , μ1 and l. We draw some boundaries Li (i=1, 2, 3, 4) in the (μ1 , μ0 )-plane corresponding to n = 2, 5, 10 and 15 taking l = 2.66 in each case. These boundaries are the (5) (4) (4) (5) straight lines with negative slopes, where E3 > E1 on the right and E1 > E3 on the left of each line (see Fig. 5).

φ3(5) versus φ2(4) case:

For this case the boundary R5,4 = 1 is given by 3,2

μ0 =

6(r − 2t )nμ1 + 2(r − t )n2 + 3(3r − 11t )n + (r + 11t ) + 3(r − t )ln + 3(r − 7t )l . 3(4t − r)

J.R. Sharma et al. / Applied Mathematics and Computation 269 (2015) 520–535

529

200

150

µ0

100

L1 L2 L3

50

L4

0

0

50

100

150

200

µ1 Fig. 6. Boundary lines in the (μ1 , μ0 )-plane for E24 and E35 .

We draw the boundaries Li (i=1, 2, 3, 4) in the (μ1 , μ0 )-plane using the same set of values of n and l as in the previous case in (5) (4) (5) (4) order to compare the efficiencies E3 and E2 . These boundaries are the straight lines with positive slopes, where E3 > E2 on (4)

the below and E2 (5)

(5)

> E3

on the above of each line (see Fig. 6).

(5)

φ3 versus φ4 case:

(5)

It is enough to compare the corresponding expressions of C3

(5)

E3

(5)

> E4

(5)

and C4

in (26) and (27). Simple calculations show that

for all l ≥ 1 and n ≥ 4.

(5)

φ3 versus φ5(5) case:

(5)

(5)

On comparing the corresponding values of C3 and C5 Below we summarize the above proved results: Theorem 2. (a) For all μ0 > 0, μ1 > 0 and l ≥ 1, we have: (5) (4) (i) E1 > E1 for n ≥ 27. (5)

(5)

(ii) E1

> E4

for n ≥ 5,

(iii) E1

> E5

for n ≥ 5,

(iv) E2

> E1

for n ≥ 26.

(v) E2

> E1

for n ≥ 2.

(vi) E2

> E3

for n ≥ 2,

(5) (5) (5) (5) (5)

(5) (4) (5) (5) (5)

(vii) E2

> E4

for n ≥ 4,

(viii) E2

> E5

for n ≥ 4,

(ix) E3

> E1

for n ≥ 26.

(x) E3

> E1

for n ≥ 2,

(xi) E3

> E4

for n ≥ 4,

(5) (5) (5) (5) (5)

(5) (4) (5) (5) (5)

(xii) E3 > E5 for n ≥ 5. (b) For μ0 > 0, μ1 > 0, l ≥ 1 and 2 ≤ n ≤ 26, we have: (5) (4) (i) E1 > E1 for μ0 > m1 , where m1 =

12(2t−r)nμ1 +2(2t−r)n2 +3(22t−7r)n+(11r−10t )+3(2t−r)ln+(42t−15r)l . 12(r−2t )

(c) For μ0 > 0, μ1 > 0, l ≥ 1 and 2 ≤ n ≤ 25, we have: (5) (4) (i) E2 > E1 for μ0 > m2 ,

(5)

in (26) and (28), it follows that E3

(5)

> E5

for all l ≥ 1 and n ≥ 5.

530

J.R. Sharma et al. / Applied Mathematics and Computation 269 (2015) 520–535

(5)

(ii) E3

(4)

for μ0 > m3 , 12(2t−r)nμ1 +2(2t−r)n2 +3(22t−7r)n+(11r−34t )+3(2t−r)ln+(42t−15r)l , 12(r−2t ) 12(2t−r)nμ1 +2(2t−r)n2 +3(22t−7r)n+11(r−2t )+3(2t−r)ln+(42t−15r)l . 12(r−2t ) > E1

where m2 = m3 =

(d) For μ0 > 0, μ1 > 0, l ≥ 1 and n ≥ 2, we have: (5) (4) (i) E1 > E2 for μ0 < m4 , (5)

(ii) E2

(5)

(iii) E3

(4)

> E2

(4)

for μ0 < m5 ,

for μ0 < m6 , 6(r−2t )nμ1 +2(r−t )n2 +3(3r−11t )n+(r+5t )+3(r−t )ln+3(r−7t )l , 3(4t−r) 6(r−2t )nμ1 +2(r−t )n2 +3(3r−11t )n+(r+17t )+3(r−t )ln+3(r−7t )l , 3(4t−r) 6(r−2t )nμ1 +2(r−t )n2 +3(3r−11t )n+(r+11t )+3(r−t )ln+3(r−7t )l . 3(4t−r) > E2

where m4 = m5 = m6 =

Remark 1. It is clear from Figs. 1–6 that the efficiency regions of the proposed composite Newton–Traub methods increase in size with the increasing value of n as compared with the efficiency regions of existing fourth order methods which, however, (5) (4) decrease in size. That means Ei > E j , (i = 1, 2, 3; j = 1, 2) in a wide region of (μ1 , μ0 )-plane with increasing n. This shows that the proposed methods are more efficient in case of the systems with large dimensions. Remark 2. The convergence order of the composite Newton–Traub methods is increased from 3 of the basic Traub’s oneparameter family of methods to 5 at the cost of one additional function evaluation. This clearly points to a high computational efficiency of Newton–Traub fifth order methods than Traub’s third order methods. For this reason we have not included generalized Traub’s third order methods in the comparison of computational efficiencies. 5. Numerical results (5)

To illustrate the convergence behavior and computational efficiency of the new fifth order schemes φi (4)

, (i = 1, 2, 3) we

consider some numerical examples and compare the performance with existing fourth order φ j , ( j = 1, 2) and fifth or(5)

der φk , (k = 4, 5) schemes. All computations are performed in the programming package Mathematica [22] using multipleprecision arithmetic with 4096 digits. For every method, we analyze the number of iterations (k) needed to converge to the solution such that the stopping criterion ||x(k+1) − x(k) || + ||F(x(k) )|| < 10−100 is satisfied. In numerical results, we also include CPU time used in the execution of program which is computed by the Mathematica command “TimeUsed[ ]”. The results of Theorem 2 are also verified through numerical examples. In order to do this, we need an estimation of the factors μ0 and μ1 . To claim this estimation, we express the cost of the evaluation of elementary functions in terms of products, which depends on the computer, the software and the arithmetics used (see, for example [23]). In Table 1, the elapsed CPU time (measured in milliseconds) in the computation of elementary functions and an estimation of the cost of the elementary functions in product units are displayed. The programs are performed in the processor Intel (R) Core (TM) i5-520M CPU @ 2.40 GHz (64-bit Machine) Microsoft Windows 7 Home Premium 2009 and are complied by computational software program Mathematica using multiple-precision arithmetic. It can be observed from Table 1 that for this hardware and the software, the computational cost of division with respect to multiplication is, l = 2.66. For numerical tests we consider the following problems: Problem 1. Consider the system of two equations (selected from [19]):



x1 + ex2 − cos x2 = 0, 3x1 − sin x1 − x2 = 0,

with initial value x(0) = {−1, 1}T towards the solution: r = {0, 0}T . For this problem the corresponding values of parameters μ0 and μ1 (calculated by using the Table 1) are 136.93 and 68.21 which we use in (22)–(28) for computing computational costs and efficiency indices. Note that the other parameters are (n, l ) = (2, 2.66). Problem 2. The boundary value problem (see [10]):

y + y3 = 0, y(0) = 0, y(1) = 1, Table 1 √ √ CPU time and estimation of computational cost of elementary functions, where x = 3 − 1 and y = 5. √ Functions x∗y x/y ln(x) sin(x) cos(x) arccos(x) x ex CPU time Cost

0.0495 1

0.1317 2.66

0.0601 1.22

3.8934 78.74

3.7347 75.52

4.8079 97.23

4.7908 96.88

7.8906 159.57

arctan(x) 7.6487 154.68

J.R. Sharma et al. / Applied Mathematics and Computation 269 (2015) 520–535

531

Table 2 Abscissas and weights of the Gauss–Legendre quadrature formula for n = 8. j

tj

ϖj

1 2 3 4 5 6 7 8

0.01985507175123188415821957… 0.10166676129318663020422303… 0.23723379504183550709113047… 0.40828267875217509753026193… 0.59171732124782490246973807… 0.76276620495816449290886952… 0.89833323870681336979577696… 0.98014492824876811584178043…

0.05061426814518812957626567… 0.11119051722668723527217800… 0.15685332293894364366898110… 0.18134189168918099148257522… 0.18134189168918099148257522… 0.15685332293894364366898110… 0.11119051722668723527217800… 0.05061426814518812957626567…

is studied. Consider the following partitioning of the interval [0, 1]: u0 = 0 < u1 < u2 < · · · < um−1 < um = 1, u j+1 = u j + h, h = 1/m. Let us define y0 = y(u0 ) = 0, y1 = y(u1 ), . . . , ym−1 = y(um−1 ), ym = y(um ) = 1. If we discretize the problem by using the numerical formula for second derivative

yk =

yk−1 − 2yk + yk+1 , h2

k = 1, 2, 3, . . . , m − 1,

we obtain a system of m − 1 nonlinear equations in m − 1 variables:

yk−1 − 2yk + yk+1 + h2 y3k = 0, k = 1, 2, 3, . . . m − 1.

(30)

In particular, we solve this problem for m = 5 so that n = 4 by selecting y(0) = {0.5, 0.5, 0.5, 0.5}T as the initial value. The solution

of this problem is, r = {0.21054188948074775 . . . , 0.42071046387616439 . . . , 0.62790045371805633 . . . , 0.82518822786851363 . . .}T and concrete values of the parameters (n, l, μ0 , μ1 ) are (4, 2.66, 4, 0.75). Problem 3. Consider the quadratic integral equation related with Chandrasekhar’s work [24]:

x(s) = f (s) + λx(s)



1 0

κ(s, t )x(t )dt,

(31)

which arises in the study of the radiative transfer theory, the transport of neutrons and the kinetic theory of the gases. It is studied in [25]. We consider the max-norm, the kernel κ (s, t) as a continuous function with s, t ∈ [0, 1] such that 0 < κ (s, t) < 1 and κ(s, t ) + κ(t, s) = 1. Moreover, we suppose that f(s) ∈ C[0, 1] is a given function and λ is a real number. Notice that finding a solution of (31) is equivalent to solving the equation F(x) = 0, where F: C[0, 1] → C[0, 1] and

F(x)(s) = x(s) − f (s) − λx(s) In particular, we consider

F(x)(s) = x(s) − 1 −

x(s) 4



1 0



1 0

κ(s, t )x(t )dt, x ∈ C[0, 1], s ∈ [0, 1].

s x(t )dt, s+t

x ∈ C[0, 1], s ∈ [0, 1].

(32)

Transforming the above equation into a finite-dimensional problem by using Gauss–Legendre quadrature formula given as



1 0

f (t )dt ≈

n 

 j f (t j ),

j=1

where the abscissas tj and the weights ϖj are determined for n = 8. Denoting the approximation of x(ti ) by xi (i = 1, 2, . . . , 8), we obtain the system of nonlinear equations

x i − xi

8 

ai j x j − 1 = 0, where ai j =

j=1

ti  j , i = 1, 2, . . . , 8, 4(ti + t j )

wherein the abscissas tj and the weights ϖj are known and given in Table 2 for n = 8. The initial approximation assumed is

x(0) = {−3, −3, · · · · ·, −3}T for obtaining the solution of this problem, r = {1.0217197314617265 . . . , 1.0731863817335818 . . . , 1.1257248936565282. . . , 1.1697533121691148. . . , 1.2030717513053575. . . , 1.2264908746333123. . . , 1.2415246005934998. . . , 1.2494485166934808 . . .}T . In this case the concrete values of parameters are (n, l, μ0 , μ1 ) = (8, 2.66, 12.5, 1.875), which we use in (22)–(28) for the calculation of costs and efficiency indices. Problem 4. Considering the system of thirteen equations (selected from [26]): 13  j=1, j =i

x j − e−xi = 0, 1 ≤ i ≤ 13,

532

J.R. Sharma et al. / Applied Mathematics and Computation 269 (2015) 520–535 Table 3 Comparison of the performances of methods. Methods

C(i p)

E(i p)

k

ρ k ± ρ k

6 6 5 5 5 5 5

4 ± 5.39 × 10−19 4 ± 2.89 × 10−7 5 ± 1.67 × 10−9 5 ± 1.94 × 10−5 5 ± 1.01 × 10−5 5 ± 1.87 × 10−15 5 ± 8.76 × 10−7

1125.02 853.50 1140.34 1136.34 1138.34 1128.34 1134.00

1.0012330 1.0016256 1.0014124 1.0014173 1.0014148 1.0014274 1.0014203

0.40067 0.35367 0.34867 0.34352 0.34824 0.40546 0.37967

5 5 4 4 4 4 4

4 ± 2.52 × 10−5 4 ± 3.07 × 10−5 5 ± 2.44 × 10−4 5 ± 5.75 × 10−4 5 ± 4.71 × 10−4 5 ± 2.80 × 10−4 5 ± 1.47 × 10−4

173.88 185.20 220.52 212.52 216.52 217.80 214.48

1.0080046 1.0075135 1.0073251 1.0076019 1.0074609 1.0074169 1.0075321

0.09834 0.11556 0.12667 0.10304 0.11952 0.12375 0.11336

6 6 5 5 5 5 5

4 ± 4.66 × 10−4 4 ± 3.00 × 10−4 5 ± 1.42 × 10−3 5 ± 1.42 × 10−3 5 ± 1.42 × 10−3 5 ± 9.91 × 10−4 5 ± 1.13 × 10−3

958.32 1067.52 1115.60 1099.60 1107.60 1323.28 1194.08

1.0014476 1.0012995 1.0014437 1.0014647 1.0014541 1.0012170 1.0013488

0.76524 0.87923 0.78112 0.66568 0.67334 0.89453 0.79646

5 5 4 4 4 4 4

4 ± 1.07 × 10−23 4 ± 6.29 × 10−24 5 ± 1.87 × 10−31 5 ± 1.10 × 10−23 5 ± 4.41 × 10−25 5 ± 2.21 × 10−31 5 ± 3.28 × 10−34

5705.73 5531.01 6091.31 6065.31 6078.31 7251.69 6597.79

1.0002430 1.0002507 1.0002643 1.0002654 1.0002648 1.0002220 1.0002440

2.86535 2.71433 2.34456 2.31924 2.33210 3.17225 2.60524

5 5 4 4 4 4 4

4 ± 1.86 × 10−22 4 ± 1.24 × 10−25 5 ± 2.03 × 10−33 5 ± 2.74 × 10−32 5 ± 1.49 × 10−32 5 ± 3.92 × 10−34 5 ± 8.47 × 10−34

12459.80 13484.40 13333.00 13293.00 13313.00 18010.60 15488.40

1.0001113 1.0001028 1.0001207 1.0001211 1.0001209 1.0000894 1.0001039

CPU-Time

Problem 1

φ1(4) φ2(4) φ1(5) φ2(5) φ3(5) φ4(5) φ5(5)

Problem 2

φ1(4) φ2(4) φ1(5) φ2(5) φ3(5) φ4(5) φ5(5)

Problem 3

φ1(4) φ2(4) φ1(5) φ2(5) φ3(5) φ4(5) φ5(5)

Problem 4

φ1(4) φ2(4) φ1(5) φ2(5) φ3(5) φ4(5) φ5(5)

Problem 5

φ1(4) φ2(4) φ1(5) φ2(5) φ3(5) φ4(5) φ5(5)

22.7465 23.9976 18.9750 18.6525 18.7885 24.3722 21.2665

with initial value x(0) = {1, 1, · · · · ·, 1}T towards the solution:

r = {0.077146207613064638 . . . , 0.077146207613064638 . . . , · · · · ·, 0.077146207613064638 . . .}T . Values of the parameters (n, l, μ0 , μ1 ) are (13,2.66,78.74,6.057). Problem 5. Lastly, considering the system of twenty equations (selected from [19]):



xi − cos 2xi −

20 

 xj

= 0, 1 ≤ i ≤ 20.

j=1

In this problem a closer choice of initial approximation to the required solution is very much needed since the problem has many solution vectors with the samevalue of each component of magnitude less than one in every solution vector. That means √ 20 (0) = {−0.9, −0.9, · · · · ·, −0.9}T 2 each solution vector satisfies r = i=1 |ri | < 20. We choose the initial approximation x for obtaining the solution,

r = {−0.89797814194212824 . . . , −0.89797814194212824 . . . , · · · · ·, −0.89797814194212824 . . .}T . Here the concrete values of parameters used in (22)–(28) are (n, l, μ0 , μ1 ) = (20, 2.66, 96.88, 4.862). (4)

In Table 3, we exhibit numerical results obtained for the considered Problems 1–5 by implementing the methods φi (5)

, (i =

1, 2) and φ j , ( j = 1, 2, . . . , 5). Displayed in the table are the necessary iterations (k), the approximate computational order of

J.R. Sharma et al. / Applied Mathematics and Computation 269 (2015) 520–535

533

Table 4 Comparison between E1(4) , E2(4) , E1(5) , E2(5) and E3(5) according to Theorem 2(b–d). Problem

μ0

m1

m2

m3

m4

m5

m6

Comparison between efficiencies as per Theorem 2(b–d) case (b-i)

1 2 3 4 5

136.93 4 12.5 78.74 96.88

−120.53 18.48 13.67 −48.58 −78.99

−126.74 12.27 7.46 −54.79 −85.21

−123.64 15.38 10.57 −51.58 −82.10

47.86 2.36 30.94 108.99 235.25

50.25 4.74 33.32 111.38 237.63

49.06 3.55 32.13 110.19 236.44

(5)

E1 E1(5) E1(5) E1(5) E1(5)

(4)

> E1 < E1(4) < E1(4) > E1(4) > E1(4)

case (c-i) (5)

E2 E2(5) E2(5) E2(5) E2(5)

(4)

> E1 < E1(4) > E1(4) > E1(4) > E1(4)

case (c-ii) (5)

E3 E3(5) E3(5) E3(5) E3(5)

(4)

> E1 < E1(4) > E1(4) > E1(4) > E1(4)

case (d-i) (5)

E1 E1(5) E1(5) E1(5) E1(5)

(4)

< E2 < E2(4) > E2(4) > E2(4) > E2(4)

( p)

case (d-ii) (5)

E2 E2(5) E2(5) E2(5) E2(5)

(4)

< E2 > E2(4) > E2(4) > E2(4) > E2(4)

case (d-iii) E3(5) E3(5) E3(5) E3(5) E3(5)

< E2(4) < E2(4) > E2(4) > E2(4) > E2(4)

( p)

convergence (ρ k ± ρ k ), the computational cost (Ci ) in terms of products, the computational efficiency (Ei ) and the mean elapsed CPU time (CPU-Time). Computational cost and efficiency are calculated according to the corresponding expressions given by (22)–(28) using the values of parameters n, μ0 and μ1 as shown in the end of each problem and taking l = 2.66 in each case. The mean elapsed CPU time is calculated by taking the mean of 50 performances of the program, where we use ||x(k+1) − x(k) || + ||F(x(k) )|| < 10−100 as the stopping criterion in single performance of the program. From the numerical results, we can observe that like the existing methods the present methods show consistent convergence behavior. Calculated values of the computational order of convergence displayed in the third column of Table 3 verify the theoretical order of convergence proved in Section 3. Numerical values of the efficiency indices displayed in the penultimate column of the table also verify the results of Theorem 2(a). In order to confirm that the results obtained in Table 3 are also in accordance with the results of Theorem 2(b–d), we find mi , (i = 1, 2, . . . , 6) using the values of l, n, and μ1 for each numerical problem. The results are presented in Table 4. Comparing the efficiency results of Table 3 with that of Table 4, we find that the results are compatible with each other for each numerical problem, and hence Theorem 2(b–d) is also verified. From the numerical values of the ( p) efficiency (Ei ) and elapsed CPU Time (CPU-Time), we can observe that the method with large efficiency uses less computing time than the method with small efficiency. This shows that the efficiency results are in complete agreement with the CPU time utilized in the execution of program. 6. Basins of attraction The study of dynamical behavior of the rational function associated to an iterative method gives important information about convergence and stability of the method, see for example [27,28]. To start with, let us recall some basic dynamical concepts of rational function. Let φ : Rn → Rn be a vectorial rational function. The orbit of a point x(0) ∈ Rn is defined as the set of successive images of x(0) by the vectorial rational function, {x(0) , φ(x(0) ), . . . , φ m (x(0) ), . . .}. The dynamical behavior of the orbit of a point of Rn can be classified depending on its asymptotic behavior. In this way, a point x∗ ∈ Rn is a fixed point of φ if φ(x∗ ) = x∗ . The practical tool for classifying the stability of fixed points is given as (see [27]): Let x∗ be a fixed point of φ. Then, (a) If | (b) If | (c) If |

∂φi (x∗ ) 1 ∗ n ∂ x j |< n for all i, j ∈ {1, . . . , n}, then x ∈ R is attracting. ∂φi (x∗ ) ∗ n ∂ x j |= 0, for all i, j ∈ {1, . . . , n}, then x ∈ R is superattracting. ∂φi (x∗ ) 1 ∗ n ∂ x j |> n for all i, j ∈ {1, . . . , n}, then x ∈ R is unstable and lies at the Julia set,

where φi (x), i = 1, 2, . . . , n, are the coordinate functions of the fixed point multivariate function φ. If is an attracting fixed point of the rational function φ, its basin of attraction A(x∗ ) is defined as the set of pre-images of any order such that x∗

A(x∗ ) = {x(0) ∈ Rn : φ m (x(0) ) → x∗ , m → ∞}. In the same way as in the scalar case [28], the set of points whose orbits tend to an attracting fixed point x∗ is defined as the Fatou set, F (φ). The complementary set, the Julia set J (φ), is the closure of the set consisting of its repelling fixed points, and establishes the borders between the basins of attraction. To study dynamical behavior, we analyze the basins of attraction of the methods in previous sections on the polynomial system (see [27])



x21 − 1 = 0, x22 − 1 = 0,

with roots {{1, 1}T , {1, −1}T , {−1, 1}T , {−1, −1}T }. For generating basins of attraction associated with roots of nonlinear system of equations, we take a square [−2, 2] × [−2, 2] of 1024 × 1024 points, which contains all roots of concerned nonlinear system of equations and we apply the iterative method starting in every point in the square. We assign a color to each point according to the root to which the corresponding orbit of the iterative method, starting from the point, converges. If the corresponding orbit

534

J.R. Sharma et al. / Applied Mathematics and Computation 269 (2015) 520–535

Fig. 7. Basins of attraction for system of equations {x21 − 1 = 0, x22 − 1 = 0} for various methods.

does not reach any root of the polynomial, with tolerance 10−3 in a maximum of 25 iterations, we mark those points with black color. For the given test problem, it can be observed in Fig. 7 that all the roots of the polynomial system have their respective basins of attraction with different colors. Also the Julia set can be seen as black lines of unstable behavior. It can be easily observed that (5) (4) (4) the method φ4 (Fig. 7(f)) seems to produce almost larger basins of attraction than the methods φ1 (Fig. 7(a)), φ2 (Fig. 7(b)),

φ1(5) (Fig. 7(c)), φ2(5) (Fig. 7(d)), φ3(5) (Fig. 7(e)) and φ5(5) (Fig. 7(g)).

J.R. Sharma et al. / Applied Mathematics and Computation 269 (2015) 520–535

535

7. Concluding remarks In the foregoing study, we have developed a one-parameter family of fifth-order iterative methods for solving systems of nonlinear equations. The scheme of the family consists of three steps. Of these three steps the first two are the steps of generalized Traub’s cubically convergent method and third is Newton type step. Hence, the name composite Newton–Traub methods. Direct computation by Taylor’s expansion is used to prove the local convergence order of the methods. The computational efficiency is discussed exhaustively. Then, a comparison between the efficiencies of some particular methods of the new family with some existing methods is performed. It is shown that the proposed methods are more efficient than existing methods, especially when applied for solving the systems of large dimensions. To illustrate the new techniques five numerical examples are presented and completely solved. Computational results have justified robust and efficient convergence behavior of the proposed methods. Similar numerical experimentations, carried out for a number of problems of different type, confirmed the above conclusions to a great extent. Acknowledgments One of the authors, Nitin Kalra, acknowledges the I.K. Gujral Punjab Technical University, Jalandhar for providing research support. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28]

J.F. Traub, Iterative Methods for the Solution of Equations, Prentice-Hall, Englewood Cliffs, New Jersey, 1964. J.M. Ortega, W.C. Rheinboldt, Iterative Solution of Nonlinear Equations in Several Variables, Academic Press, New York, 1970. C.T. Kelley, Solving Nonlinear Equations with Newton’s Method, SIAM, Philadelphia, 2003. I.K. Argyros, Computational Theory of Iterative Methods, in: C.K. Chui, L. Wuytack (Eds.), Studies in Computational Mathematics, 15, Elsevier Publ. Comp., New York, 2007. ´ B. Neta, L.D. Petkovic, ´ J. Džunic, ´ Multipoint Methods for Solving Nonlinear Equations, Elsevier, Boston, 2013. M.S. Petkovic, M. Frontini, E. Sormani, Third-order methods from quadrature formulae for solving systems of nonlinear equations, Appl. Math. Comput. 149 (2004) 771– 782. H.H.H. Homeier, A modified Newton method with cubic convergence: the multivariable case, J. Comput. Appl. Math. 169 (2004) 161–169. A. Cordero, J.R. Torregrosa, Variants of Newton’s method for functions of several variables, Appl. Math. Comput. 183 (2006) 199–208. M.A. Noor, M. Waseem, Some iterative methods for solving a system of nonlinear equations, Comput. Math. Appl. 57 (2009) 101–106. 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. A. Cordero, J.R. Torregrosa, Variants of Newton’s method using fifth-order quadrature formulas, Appl. Math. Comput. 190 (2007) 686–698. M.T. Darvishi, A. Barati, A third-order Newton-type method to solve system of nonlinear equations, Appl. Math. Comput. 187 (2007) 630–635. A. Cordero, E. Martínez, J.R. Torregrosa, Iterative methods of order four and five for systems of nonlinear equations, Appl. Math. Comput. 231 (2009) 541–551. M.T. Darvishi, A. Barati, A fourth-order method from quadrature formulae to solve systems of nonlinear equations, Appl. Math. Comput. 188 (2007) 257–261. A. Cordero, J.L. Hueso, E. Martínez, J.R. Torregrosa, A modified Newton–Jarratt’s composition, Numer. Algor. 55 (2010) 87–99. P. Jarratt, Some fourth order multipoint iterative methods for solving equations, Math. Comput. 20 (1966) 434–437. J.R. Sharma, R.K. Guha, R. Sharma, An efficient fourth order weighted-Newton method for systems of nonlinear equations, Numer. Algor. 62 (2013) 307–323. A. Cordero, J.L. Hueso, E. Martínez, J.R. Torregrosa, Increasing the convergence order of an iterative method for nonlinear systems, Appl. Math. Lett. 25 (2012) 2369–2374. J.R. Sharma, P. Gupta, An efficient fifth order method for solving systems of nonlinear equations, Comput. Math. Appl. 67 (2014) 591–601. ´ On a general class of multipoint root-finding methods of high computational efficiency, SIAM J. Numer. Anal. 49 (2011) 1317–1319. M.S. Petkovic, L.O. Jay, A note on Q-order of convergence, BIT 41 (2001) 422–429. S. Wolfram, The Mathematica Book, fifth edition, Wolfram Media, 2003. L. Fousse, G. Hanrot, V. Lefvre, P. Plissier, P. Zimmermann, MPFR: a multiple-precision binary floating-point library with correct rounding, ACM Trans. Math. Softw. 33 (2) (2007).15. Art. 13 S. Chandrasekhar, Radiative Transfer, Dover, New York, 1960. I.K. Argyros, Quadratic equations and applications to Chandrasekhar’s and related equations, Bull. Aust. Math. Soc. 32 (1985) 275–292. 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. A. Cordero, F. Soleymani, J.R. Torregrosa, Dynamical analysis of iterative methods for nonlinear systems or how to deal with the dimension? Appl. Math. Comput. 244 (2014) 398–412. Á.A. Magreñán, A new tool to study real dynamics: the convergence plane, Appl. Math. Comput. 248 (2014) 215–224.