# On Some Efficient Techniques for Solving Systems of Nonlinear Equations

Authors: Janak Raj Sharma; Puneet Gupta
Journal: Advances in Numerical Analysis
 (2013)

---

## Abstract

We present iterative methods of convergence order three, five, and six for solving systems of nonlinear equations. Third-order method is composed of two steps, namely, Newton iteration as the first step and weighted-Newton iteration as the second step. Fifth and sixth-order methods are composed of three steps of which the first two steps are same as that of the third-order method whereas the third is again a weighted-Newton step. Computational efficiency in its general form is discussed and a comparison between the efficiencies of proposed techniques with existing ones is made. The performance is tested through numerical examples. Moreover, theoretical results concerning order of convergence and computational efficiency are verified in the examples. It is shown that the present methods have an edge over similar existing methods, particularly when applied to large systems of equations.


---

## Body

## 1. Introduction

Solving the systemF(x)=0 of nonlinear equations is a common and important problem in various disciplines of science and engineering [1–4]. This problem is precisely stated as follows: For a given nonlinear function F(x):D⊆Rn→Rn, where F(x)=(f1(x),f2(x),…,fn(x))t and x=(x1,x2,…,xn)t, to find a vector α=(α1,α2,…,αn)t such that F(α)=0. The solution vector α can be obtained as a fixed point of some function G:Rn→Rn by means of fixed point iteration
(1)x(k+1)=G(x(k)),k=0,1,….
One of the basic procedures for solving systems of nonlinear equations is the classical Newton’s method [4, 5] which converges quadratically under the conditions that the function F is continuously differentiable and a good initial approximation x(0) is given. It is defined by
(2)x(k+1)=G2(x(k))=x(k)-F′(x(k))-1F(x(k)),k=0,1,2,…,
where F′(x)-1 is the inverse of first Fréchet derivative F′(x) of the function F(x). Note that this method uses one function, one first derivative, and one matrix inversion evaluations per iteration.In order to improve the order of convergence of Newton’s method, many modifications have been proposed in literature, for example, see [6–18] and references therein. For a system of n equations in n unknowns, the first Fréchet derivative is a matrix with n2 evaluations whereas the second Fréchet derivative has n2(n+1)/2 evaluations. Thus, the methods such as those developed in [6–8] with second derivative are considered less efficient from a computational point of view.In quest of efficient methods without using second Fréchet derivative, a variety of third and higher order methods have been proposed in recent years. For example, Frontini and Sormani in [9] and Homeier in [10] developed some third order methods each requiring the evaluations of one function, two first-order derivatives, and two matrix inversions per iteration. Darvishi and Barati [11] presented a fourth-order method which uses two functions, three first derivatives, and two matrix inversions. Cordero and Torregrosa [12] developed two variants of Newton’s method with third-order convergence; one of the variants requires the evaluations of one function, three first derivatives, and two matrix inversions whereas, the other requires one function, two first derivatives, and two matrix inversions. Noor and Waseem [13] have developed two third-order methods of which one method requires one function, two first derivatives, and two matrix inversions, and the other requires one function, three first derivatives, and two matrix inversions. Cordero et al. [14] presented fourth and fifth-order methods requiring two functions, two first derivatives, and one matrix inversion and three functions, two first derivatives, and one matrix inversion, respectively. Cordero et al. [15] developed a sixth-order method which requires two functions, two first derivatives, and two matrix inversions. Grau-Sánchez et al. [16] proposed the methods with third, fourth, and fifth-order convergence; third-order method with one function, two first derivatives, and two matrix inversions, fourth method with three functions, one first derivative, and one matrix inversion, and fifth-order method with two functions, two first derivatives, and two matrix inversions. Grau-Sánchez et al. [17] have also developed fourth and sixth-order methods requiring the evaluations of one function, one divided difference, one first derivative, and two matrix inversions for the fourth-order method and two functions, one divided difference, one first derivative, and two matrix inversions for the sixth-order method. Sharma et al. [18] proposed a fourth-order method which requires one function, two first derivatives, and two matrix inversions per iteration.In this paper, our aim is to develop second derivative free iterative methods that may satisfy the basic requirements of generating the quality numerical algorithms, that is, the algorithms which have (i) high convergence speed, (ii) minimum computational cost, and (iii) simple structure. Taking into account these considerations, we here devise the methods of third, fifth, and sixth-order of convergence. Third order method is composed of two steps namely, Newton’s and weighted-Newton steps, and per iteration it requires the evaluations of one function, two first derivatives, and one matrix inversion. Fifth and sixth-order methods are composed of three steps of which the first step is Newton’s step and the last two are weighted-Newton steps. Per iteration, both of these schemes require the evaluations of two functions, two first derivatives, and one matrix inversion.We summarize the contents of the paper as follows. In Section2, third-order scheme is developed and its convergence analysis is studied. In Section 3, fifth-order scheme is developed and its convergence analysis is studied. Sixth-order scheme with its convergence analysis is presented in Section 4. In Section 5, the computational efficiency of new methods is discussed and is compared with the methods which lie in the same category. Various numerical examples are considered in Section 6 to show the consistent convergence behavior of the methods and to verify the theoretical results. Section 7 contains the concluding remarks.


## 2. Development of Third-Order Method

Since our objective is to develop a method which accelerates the convergence of Newton’s method with minimum number of evaluations. Thus we consider two-step iterative scheme of the type:(3)y(k)=x(k)-F′(x(k))-1F(x(k)),x(k+1)=x(k)-(a1I+a2F′(x(k))-1F′(y(k)))×F′(x(k))-1F(x(k)),
where a1 and a2 are some parameters and I denotes the n×n identity matrix. In order to explore the convergence property of (3), we recall the following result of Taylor’s expansion on vector functions (see [2]).Lemma 1.
LettingF:D⊆Rn→Rn be a p time Fréchet differentiable in a convex set D⊆Rn, then for any x,h∈Rn, the following expression holds:
(4)F(x+h)=F(x)+F′(x)h+12!F′′(x)h2+13!F′′′(x)h3+⋯+1p!F(p-1)(x)hp-1+Rp,
where
(5)||Rp||⩽1p!sup0⩽t⩽1||F(p)(x+th)||∥h∥p,hp=(h,h,…p,h).The Taylor’s expansion (4) for F(x) about x(k) is
(6)F(x)=F(x(k))+F′(x(k))(x-x(k))+12!F′′(x(k))×(x-x(k))2+13!F′′′(x(k))(x-x(k))3+14!F(iv)(x(k))(x-x(k))4+O(∥x-x(k)∥5).
Let e(k)=x(k)-α. Then setting x=α and using F(α)=0 in (6), we obtain
(7)F(x(k))=F′(α)[e(k)+A2(e(k))2+A3(e(k))3=F′(α)h+A4(e(k))4+O((e(k))5)],
where Ai=(1/i!)ΓF(i)(α)∈Li(Rn,Rn), Γ=F′(α)-1 and (e(k))i=(e(k),e(k),…i,e(k)), e(k)∈Rn.Also,(8)F′(x(k))=F′(α)[I+2A2e(k)+3A3(e(k))2=F′(α)h+4A4(e(k))3+O((e(k))4)].
Then,
(9)F′(x(k))-1=[D(e(k))-1+O((e(k))4)]Γ,
where D(e(k))=I+2A2e(k)+3A3(e(k))2+4A4(e(k))3.Inverse ofD(e(k)) is given by
(10)D(e(k))-1=I+X1e(k)+X2(e(k))2+X3(e(k))3,
where Xi (i=1 to 3) are such that they satisfy the definition
(11)D(e(k))D(e(k))-1=D(e(k))-1D(e(k))=I.
Solving the above system, we have
(12)X1=-2A2,X2=4A22-3A3,X3=-(8A23-6A2A3-6A3A2+4A4).
From the above expressions, we find that
(13)F′(x(k))-1=[I-2A2e(k)+(4A22-3A3)(e(k))2-(8A23-6A2A3-6A3A2+4A4)×(e(k))3+O((e(k))4)]Γ.
From (7) and (13), it follows that
(14)F′(x(k))-1F(x(k))=e(k)-A2(e(k))2+2(A22-A3)×(e(k))3+O((e(k))4).
Letting e-(k)=y(k)-α, then using (14) in the first step of (3), we have
(15)e-(k)=A2(e(k))2-2(A22-A3)(e(k))3+O((e(k))4).
Expanding F′(y(k)) about α, we obtain
(16)F′(y(k))=F′(α)[I+2A2e-(k)+3A3(e-(k))2+O((e-(k))3)].
Then,
(17)F′(x(k))-1F′(y(k))=I-2A2(e(k)-e-(k))+(4A22-3A3)×(e(k))2-4A22e(k)e-(k)-(8A23-6A2A3-6A3A2+4A4)×(e(k))3+O((e(k))4).
Substituting (14) and (17) in the second step of (3), we obtain the error equation as
(18)e(k+1)=(1-a1-a2)e(k)+(a1+3a2)A2(e(k))2-2a2A2e-(k)×e(k)-(2(a1+4a2)A22-(2a1+5a2)A3)×(e(k))3+O((e(k))4).
Our aim is to find values of the parameters a1 and a2 in such a way that the proposed iterative scheme (3) may produce order of convergence as high as possible. This happens when a1=3/2 and a2=-1/2. Then, the error equation (18) on using (15) becomes
(19)e(k+1)=12(4A22+A3)(e(k))3+O((e(k))4).According to the above analysis the following theorem can be formulated.Theorem 2.
Let the functionF:D⊆Rn→Rn be sufficiently differentiable in a neighborhood D of its zero α. If an initial approximation x(0) is sufficiently close to α, then the local order of convergence of the iterative scheme (3) is at least 3, provided a1=3/2 and a2=-1/2.Thus, the proposed scheme (3) finally can be written as
(20)y(k)=x(k)-F′(x(k))-1F(x(k)),x(k+1)=G1,3=x(k)-12(3I-F′(x(k))-1F′(y(k)))×F′(x(k))-1F(x(k)).
Clearly this formula uses one function, two derivatives, and one matrix inversion evaluations per iteration.


## 3. The Fifth-Order Method

Based on the two-step scheme (20), we propose the following three-step scheme:
(21)y(k)=x(k)-F′(x(k))-1F(x(k)),z(k)=x(k)-12(3I-F′(x(k))-1F′(y(k)))×F′(x(k))-1F(x(k)),x(k+1)=z(k)-(b1I+b2F′(x(k))-1F′(y(k)))×F′(x(k))-1F(z(k)),
where b1 and b2 are some parameters that are to be determined.Now to calculate the local order of convergence, lete~(k)=z(k)-α. In view of (19), error in the second step of (21) is given by
(22)e~(k)=12(4A22+A3)(e(k))3+O((e(k))4).
Expanding F(z(k)) about α,
(23)F(z(k))=F′(α)[e~(k)+O((e~(k))2)].
Using (13), (17), (22), and (23) in the third step of (21), it follows that
(24)e(k+1)=(1-b1-b2)e~(k)+2(b1+2b2)A2e(k)e~(k)-2b2A2e-(k)e~(k)-(4(b1+3b2)A22-3(b1+2b2)A3)×(e(k))2e~(k)+O((e(k))6).
It is clear that for the parameter values b1=2 and b2=-1, the terms containing (e(k))3 and (e(k))4  vanish. Thus, from the expressions of e~(k) and e-(k) the error equation (24) yields
(25)e(k+1)=3A22(4A22+A3)(e(k))5+O((e(k))6),
which shows fifth-order of convergence.From the above analysis, we can state the following theorem.Theorem 3.
Let the functionF:D⊆Rn→Rn be sufficiently differentiable in an open neighborhood D of its zero α. If an initial approximation x(0) is sufficiently close to α, then the local order of convergence of method (21) is at least 5, if b1=2 and b2=-1.The proposed algorithm (21) now can be written as
(26)y(k)=x(k)-F′(x(k))-1F(x(k)),z(k)=x(k)-12(3I-F′(x(k))-1F′(y(k)))×F′(x(k))-1F(x(k)),x(k+1)=G1,5=z(k)-(2I-F′(x(k))-1F′(y(k)))×F′(x(k))-1F(z(k)).
It is clear that this formula uses the evaluations of two functions, two derivatives and only one matrix inversion in all.


## 4. The Sixth-Order Method

Here, with the first two steps of the scheme (20) we consider the following three-step scheme:
(27)y(k)=x(k)-F′(x(k))-1F(x(k)),z(k)=x(k)-12(3I-F′(x(k))-1F′(y(k)))×F′(x(k))-1F(x(k)),x(k+1)=z(k)-(c1I+F′(x(k))-1F′(y(k))=z(k)hh-×(c2I+c3F′(x(k))-1F′(y(k))))×F′(x(k))-1F(z(k)),
where c1, c2, and c3 are some arbitrary constants.In the following lines we obtain the local order of convergence of the above proposed scheme. Using (13), (17), (22), and (23) in the third step of (27), it follows that
(28)e(k+1)=(1-c1-c2-c3)e~(k)+2(c1+2c2+3c3)×A2e(k)e~(k)-2(c2+2c3)A2e-(k)e~(k)-(4(c1+3c2+6c3)A22-3(c1+2c2+3c3)A3)×(e(k))2e~(k)+8(c2+3c3)A22e(k)e-(k)e~(k)-(c1+c2+c3)A2(e~(k))2+2(A234(c1+4c2+10c3)A23-3(c1+3c2+6c3)A2A3-3(c1+3c2+6c3)A3A2+2(c1+2c2+3c3)A4A23)×(e(k))3e~(k)+O((e(k))7).
Combining (15), (22), and (28), it is easy to prove that for the parameters c1=7/2, c2=-4, and c3=3/2, error equation (28) produces the maximum error. For this set of values the above equation reduces to
(29)e(k+1)=12(18A23+12A2A3-3A3A2)×(4A22+A3)(e(k))6+O((e(k))7),
which shows sixth-order of convergence. Thus, based on the above discussion, we can formulate the following theorem.Theorem 4.
Let the functionF:D⊆Rn→Rn be sufficiently differentiable in an open neighborhood D of its zero α. If an initial approximation x(0) is sufficiently close to α, then the local order of convergence of method (27) is at least 6, provided c1=7/2, c2=-4 and c3=3/2.Finally, the proposed scheme (27) is given by
(30)y(k)=x(k)-F′(x(k))-1F(x(k)),z(k)=x(k)-12(3I-F′(x(k))-1F′(y(k)))×F′(x(k))-1F(x(k)),x(k+1)=G1,6=z(k)-(72I-F′(x(k))-1F′(y(k))=hG1,6hh=z(k)×(4I-32F′(x(k))-1F′(y(k))))×F′(x(k))-1F(z(k)).
Like fifth-order scheme, this formula also uses two function evaluations, two derivatives evaluations, and only one matrix inversion in all.


## 5. Computational Efficiency

To obtain an assessment of the efficiency of the proposed methods we shall make use of efficiency index, according to which the efficiency of an iterative method is given byE=ρ1/C, where ρ is the order of convergence and C is the computational cost per iteration. In order to do this, we must consider all possible factors which contribute to the total cost of computation. For a system of n nonlinear equations in n unknowns, the computational cost per iteration is given by (see [17])
(31)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 number of evaluations of scalar functions of F′, say ∂fi/∂xj, 1⩽i, j⩽n, P(n) represents the number of products or quotients needed per iteration, and μ0 and μ1 are ratios between products and evaluations required to express the value of C(μ0,μ1,n) in terms of products. We suppose that a quotient is equivalent to l products.To computeF 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 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 n products for the multiplication of a vector by a scalar, while n2 products for the same operation in case of a matrix.We compare the computational efficiency of present third order methodG1,3 with existing third-order methods by Potra-Pták [3], Homeier [10], Noor-Waseem [13], and Cordero-Torregrosa [12]; the fifth-order method G1,5 with fifth-order method given by Grau-Sánchez et al. [16]; the sixth-order method G1,6 with sixth-order method developed by Cordero et al. [15]. These methods are expressed as follows.Method by Potra-Pták (
G
2,3
):
(32)
y
(
k
)
=
x
(
k
)
-
F
′
(
x
(
k
)
)
-
1
F
(
x
(
k
)
)
,
x
(
k
+
1
)
=
y
(
k
)
-
F
′
(
x
(
k
)
)
-
1
F
(
y
(
k
)
)
.Homeier’s method (
G
3,3
):
(33)
y
(
k
)
=
x
(
k
)
-
1
2
F
′
(
x
(
k
)
)
-
1
F
(
x
(
k
)
)
,
x
(
k
+
1
)
=
x
(
k
)
-
F
′
(
y
(
k
)
)
-
1
F
(
x
(
k
)
)
.Method by Noor and Waseem (
G
4,3
):
(34)
y
(
k
)
=
x
(
k
)
-
F
′
(
x
(
k
)
)
-
1
F
(
x
(
k
)
)
,
x
(
k
+
1
)
=
x
(
k
)
-
4
[
F
′
(
x
(
k
)
)
+
3
F
′
(
x
(
k
)
+
2
y
(
k
)
3
)
]
-
1
F
(
x
(
k
)
)
.Method by Cordero and Torregrosa (
G
5,3
):
(35)
y
(
k
)
=
x
(
k
)
-
F
′
(
x
(
k
)
)
-
1
F
(
x
(
k
)
)
,
x
(
k
+
1
)
=
x
(
k
)
-
3
[
2
F
′
(
3
x
(
k
)
+
y
(
k
)
4
)
-
F
′
(
x
(
k
)
+
y
(
k
)
2
)
+
2
F
′
(
x
(
k
)
+
3
y
(
k
)
4
)
]
-
1
F
(
x
(
k
)
)
.Method by Grau-Sánchez et al. (
G
2,5
):
(36)
y
(
k
)
=
x
(
k
)
-
F
′
(
x
(
k
)
)
-
1
F
(
x
(
k
)
)
,
z
(
k
)
=
x
(
k
)
-
1
2
[
F
′
(
x
(
k
)
)
-
1
+
F
′
(
y
(
k
)
)
-
1
]
F
(
x
(
k
)
)
,
x
(
k
+
1
)
=
z
(
k
)
-
F
′
(
y
(
k
)
)
-
1
F
(
z
(
k
)
)
.Method by Cordero et al. (
G
2,6
):
(37)
y
(
k
)
=
x
(
k
)
-
2
3
F
′
(
x
(
k
)
)
-
1
F
(
x
(
k
)
)
,
z
(
k
)
=
x
(
k
)
-
1
2
[
3
F
′
(
y
(
k
)
)
-
F
′
(
x
(
k
)
)
]
-
1
×
(
3
F
′
(
y
(
k
)
)
+
F
′
(
x
(
k
)
)
)
F
′
(
x
(
k
)
)
-
1
F
(
x
(
k
)
)
,
x
(
k
+
1
)
=
z
(
k
)
-
2
[
3
F
′
(
y
(
k
)
)
-
F
′
(
x
(
k
)
)
]
-
1
F
(
z
(
k
)
)
.
In addition, we also compare present methods with each other. Let us denote efficiency indices of Gi,p by Ei,p and computational cost by Ci,p. Then taking into account the above and previous considerations, we have
(38)C1,3=nμ0+2n2μ1+n6(2n2+15n+1+3l(n+3)),E1,3=31/C1,3.(39)C2,3=2nμ0+n2μ1+n6(2n2+9n-11+3l(n+3)),E2,3=31/C2,3.(40)C3,3=nμ0+2n2μ1+n3(2n2+3n-2+3l(n+1)),E3,3=31/C3,3.(41)C4,3=nμ0+2n2μ1+n3(2n2+6n+4+3l(n+1)),E4,3=31/C4,3.(42)C5,3=nμ0+4n2μ1+n3(2n2+9n+13+3l(n+1)),E5,3=31/C5,3.(43)C1,5=2nμ0+2n2μ1+n6(2n2+33n-5+3l(n+7)),E1,5=51/C1,5.(44)C2,5=2nμ0+2n2μ1+n3(2n2+6n-5+3l(n+2)),E2,5=51/C2,5.(45)C1,6=2nμ0+2n2μ1+n6(2n2+45n+1+3l(n+9)),E1,6=61/C1,6.(46)C2,6=2nμ0+2n2μ1+n3(2n2+12n+1+3l(n+2)),E2,6=61/C2,6.


### 5.1. Comparison between the Efficiencies

To compare the computational efficiencies of the iterative methodsGi,p, we consider the ratio
(47)Ri,p;j,q=logEi,plogEj,q=log(p)Cj,qlog(q)Ci,p.
It is clear that if Ri,p;j,q>1, the iterative method Gi,p is more efficient than Gj,q. Taking into account that the border between two computational efficiencies is given by Ri,p;j,q=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.G
1,3  versus  G2,3  Case. Here we study the boundary R1,3;2,3=1 expressed by μ0 as a function of μ1 and n. This boundary, on using (38) and (39) in (47), is given by
(48)μ0=nμ1+n+2.
In order to compare the efficiencies E1,3 and E2,3 in the (μ1,μ0)-plane, we present in Figure 1 some boundary lines Li (i=1,2,3,4) corresponding to the cases n=2,5,10,15. These boundaries are the lines with positive slopes, where E1,3>E2,3 on the above and E2,3>E1,3 on the below of each line.Figure 1
Boundary lines in(μ1,μ0)-plane for E1,3 and E2,3.G
1,3  versus  G3,3  Case. Using the corresponding values of C1,3 and C3,3 from (38) and (40) in (47), it follows that R1,3;3,3>1 for all n⩾5 and l⩾1. Thus, we have E1,3>E3,3 for all n⩾5 and l⩾1. This is shown graphically in Figure 2 for a particular set (μ0,μ1,l)=(1,1,2.8). The reason for taking l=2.8 will be made clear in the next section.Figure 2
Plots forEvalues of G1,3 and G3,3 for n⩾5.G
1,3  versus  G4,3 and G1,3versus  G5,3  Cases. From (47) on using (38), (41) and then (38), (42), it is easy to prove that R1,3;4,3>1 and R1,3;5,3>1 for all n⩾2 and l⩾1. Thus, we conclude that E1,3 is higher than E4,3 and E5,3 for all n⩾2 and l⩾1. This situation is shown graphically in Figure 3 taking the same set of values of (μ0,μ1,l) as in the previous case.Figure 3
Plots forEvalues of G1,3, G4,3, and G5,3 for n⩾2.G
1,3  versus  G1,5  Case. Here we study the boundary R1,5;1,3=1 expressed by μ0 as a function of μ1, l, and n. This boundary is given by
(49)μ0=(2n2(r-s)+3n((r-s)(4μ1+l)+5r-11s)hhh+3l(3r-7s)+r+5s2n2)(6(2s-r))-1,
where r=ln(5) and s=ln(3). In order to compare the efficiencies E1,3 and E1,5 in the (μ1,μ0)-plane, we present in Figure 4 some boundary lines Li (i=1,2,3,4) corresponding to the cases n=2,5,10,15 taking l=2.8 in each case. These boundaries are the straight lines with positive slopes, where E1,3>E1,5 on the above and E1,5>E1,3 on the below of each line.Figure 4
Boundary lines in(μ1,μ0)-plane for E1,3 and E1,5.G
1,3  versus  G1,6  Case. In this case the boundary R1,6;1,3=1 is expressed by
(50)μ0=2n2t+3n(4tμ1+tl+5(t-2s))+9l(t-2s)+t6(s-t),
where s=ln(3) and t=ln(2). Similar to the case above, here we also present the boundaries Li  (i=1,2,3,4) using the same set of  n  and  l  in (μ1,μ0)-plane to compare the efficiencies E1,3 and E1,6. These boundaries are the lines with positive slopes, where E1,3>E1,6 on the above and E1,6>E1,3 on the below of each line (see Figure 5).Figure 5
Boundary lines in(μ1,μ0)-plane for E1,3 and E1,6.G
1,5  versus  G1,6  Case. For this case the boundary R1,6;1,5=1 is given by
(51)μ0=(2n2(u-r)+3n((u-r)(4μ1+l)+11u-15r)μ0h(+3l(7u-9r)-5u-r2n2)(12(r-u))-1,
where r=ln(5) and u=ln(6). Here we also draw boundaries in (μ1,μ0)-plane with the same set of values of n and l. The boundaries Li are the lines with negative slopes, where E1,6>E1,5 on the above and E1,5>E1,6 on the below of each line (see Figure 6).Figure 6
Boundary lines in(μ1,μ0)-plane for E1,5 and E1,6.G
1,5  versus  G2,5  Case. Substituting the corresponding values of C1,5 and C2,5 from (43) and (44) in (47), we see that R1,5;2,5>1 for all n⩾10 and l⩾1, which further implies that E1,5>E2,5 for all n⩾10, l⩾1. This result is verified graphically in Figure 7 by using the values (μ0,μ1,l)=(1,1,2.8).Figure 7
Plots forEvalues of G1,5 and G2,5 for n⩾10.G
1,6  versus  G2,6  Case. In this case substitution of the corresponding values of C1,6 and C2,6 from (45) and (46) in (47) results in R1,6;2,6>1 for all n⩾10 and l⩾1. Thus, we analyze that E1,6>E2,6 for all n⩾10 and l⩾1. This is shown graphically in Figure 8 by using the same values of (μ0,μ1,l) that are used in the previous case.Figure 8
Plots forEvalues of G1,6 and G2,6 for n⩾10.We summarize the above results in following theorem.Theorem 5.
(i) For allμ0>0, μ1>0, and l⩾1 we have:
(a)E1,3>E3,3∀n⩾5,
(b)E1,3>E4,3 and E1,3>E5,3∀n⩾2,
(c)E1,5>E2,5 and E1,6>E2,6∀n⩾10.
Otherwise, the comparison depends onn, μ0, μ1, and l.
(ii) For alln⩾2, we have:
(a)E1,3>E2,3 for μ0>m1,
(b)E1,5>E1,3 for μ0<m2,
(c)E1,6>E1,3 for μ0<m3,
(d)E1,6>E1,5 for μ0>m4,
wherem1=nμ1+n+2, m2=(2n2(r-s)+3n((r-s)(4μ1+l)+5r-11s)+3l(3r-7s)+r+5s)/6(2s-r), m3=(2n2t+3n(4tμ1+tl+5(t-2s))+9l(t-2s)+t)/6(s-t), and m4=(2n2(u-r)+3n((u-r)(4μ1+l)+11u-15r)+3l(7u-9r)-5u-r)/12(r-u).


## 5.1. Comparison between the Efficiencies

To compare the computational efficiencies of the iterative methodsGi,p, we consider the ratio
(47)Ri,p;j,q=logEi,plogEj,q=log(p)Cj,qlog(q)Ci,p.
It is clear that if Ri,p;j,q>1, the iterative method Gi,p is more efficient than Gj,q. Taking into account that the border between two computational efficiencies is given by Ri,p;j,q=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.G
1,3  versus  G2,3  Case. Here we study the boundary R1,3;2,3=1 expressed by μ0 as a function of μ1 and n. This boundary, on using (38) and (39) in (47), is given by
(48)μ0=nμ1+n+2.
In order to compare the efficiencies E1,3 and E2,3 in the (μ1,μ0)-plane, we present in Figure 1 some boundary lines Li (i=1,2,3,4) corresponding to the cases n=2,5,10,15. These boundaries are the lines with positive slopes, where E1,3>E2,3 on the above and E2,3>E1,3 on the below of each line.Figure 1
Boundary lines in(μ1,μ0)-plane for E1,3 and E2,3.G
1,3  versus  G3,3  Case. Using the corresponding values of C1,3 and C3,3 from (38) and (40) in (47), it follows that R1,3;3,3>1 for all n⩾5 and l⩾1. Thus, we have E1,3>E3,3 for all n⩾5 and l⩾1. This is shown graphically in Figure 2 for a particular set (μ0,μ1,l)=(1,1,2.8). The reason for taking l=2.8 will be made clear in the next section.Figure 2
Plots forEvalues of G1,3 and G3,3 for n⩾5.G
1,3  versus  G4,3 and G1,3versus  G5,3  Cases. From (47) on using (38), (41) and then (38), (42), it is easy to prove that R1,3;4,3>1 and R1,3;5,3>1 for all n⩾2 and l⩾1. Thus, we conclude that E1,3 is higher than E4,3 and E5,3 for all n⩾2 and l⩾1. This situation is shown graphically in Figure 3 taking the same set of values of (μ0,μ1,l) as in the previous case.Figure 3
Plots forEvalues of G1,3, G4,3, and G5,3 for n⩾2.G
1,3  versus  G1,5  Case. Here we study the boundary R1,5;1,3=1 expressed by μ0 as a function of μ1, l, and n. This boundary is given by
(49)μ0=(2n2(r-s)+3n((r-s)(4μ1+l)+5r-11s)hhh+3l(3r-7s)+r+5s2n2)(6(2s-r))-1,
where r=ln(5) and s=ln(3). In order to compare the efficiencies E1,3 and E1,5 in the (μ1,μ0)-plane, we present in Figure 4 some boundary lines Li (i=1,2,3,4) corresponding to the cases n=2,5,10,15 taking l=2.8 in each case. These boundaries are the straight lines with positive slopes, where E1,3>E1,5 on the above and E1,5>E1,3 on the below of each line.Figure 4
Boundary lines in(μ1,μ0)-plane for E1,3 and E1,5.G
1,3  versus  G1,6  Case. In this case the boundary R1,6;1,3=1 is expressed by
(50)μ0=2n2t+3n(4tμ1+tl+5(t-2s))+9l(t-2s)+t6(s-t),
where s=ln(3) and t=ln(2). Similar to the case above, here we also present the boundaries Li  (i=1,2,3,4) using the same set of  n  and  l  in (μ1,μ0)-plane to compare the efficiencies E1,3 and E1,6. These boundaries are the lines with positive slopes, where E1,3>E1,6 on the above and E1,6>E1,3 on the below of each line (see Figure 5).Figure 5
Boundary lines in(μ1,μ0)-plane for E1,3 and E1,6.G
1,5  versus  G1,6  Case. For this case the boundary R1,6;1,5=1 is given by
(51)μ0=(2n2(u-r)+3n((u-r)(4μ1+l)+11u-15r)μ0h(+3l(7u-9r)-5u-r2n2)(12(r-u))-1,
where r=ln(5) and u=ln(6). Here we also draw boundaries in (μ1,μ0)-plane with the same set of values of n and l. The boundaries Li are the lines with negative slopes, where E1,6>E1,5 on the above and E1,5>E1,6 on the below of each line (see Figure 6).Figure 6
Boundary lines in(μ1,μ0)-plane for E1,5 and E1,6.G
1,5  versus  G2,5  Case. Substituting the corresponding values of C1,5 and C2,5 from (43) and (44) in (47), we see that R1,5;2,5>1 for all n⩾10 and l⩾1, which further implies that E1,5>E2,5 for all n⩾10, l⩾1. This result is verified graphically in Figure 7 by using the values (μ0,μ1,l)=(1,1,2.8).Figure 7
Plots forEvalues of G1,5 and G2,5 for n⩾10.G
1,6  versus  G2,6  Case. In this case substitution of the corresponding values of C1,6 and C2,6 from (45) and (46) in (47) results in R1,6;2,6>1 for all n⩾10 and l⩾1. Thus, we analyze that E1,6>E2,6 for all n⩾10 and l⩾1. This is shown graphically in Figure 8 by using the same values of (μ0,μ1,l) that are used in the previous case.Figure 8
Plots forEvalues of G1,6 and G2,6 for n⩾10.We summarize the above results in following theorem.Theorem 5.
(i) For allμ0>0, μ1>0, and l⩾1 we have:
(a)E1,3>E3,3∀n⩾5,
(b)E1,3>E4,3 and E1,3>E5,3∀n⩾2,
(c)E1,5>E2,5 and E1,6>E2,6∀n⩾10.
Otherwise, the comparison depends onn, μ0, μ1, and l.
(ii) For alln⩾2, we have:
(a)E1,3>E2,3 for μ0>m1,
(b)E1,5>E1,3 for μ0<m2,
(c)E1,6>E1,3 for μ0<m3,
(d)E1,6>E1,5 for μ0>m4,
wherem1=nμ1+n+2, m2=(2n2(r-s)+3n((r-s)(4μ1+l)+5r-11s)+3l(3r-7s)+r+5s)/6(2s-r), m3=(2n2t+3n(4tμ1+tl+5(t-2s))+9l(t-2s)+t)/6(s-t), and m4=(2n2(u-r)+3n((u-r)(4μ1+l)+11u-15r)+3l(7u-9r)-5u-r)/12(r-u).


## 6. Numerical Results

In this section, some numerical problems are considered to illustrate the convergence behavior and computational efficiency of the proposed methods. The performance is compared with the existing methodsG2,3, G3,3, G4,3, G5,3, G2,5, and G2,6 that we have introduced in previous section. All computations are performed in the programming package Mathematica [19] using multiple-precision arithmetic with 4096 digits. For every method, we analyze the number of iterations (k) needed to converge to the solution such that (||x(k+1)-x(k)||+||F(x(k))||)<10-200. In numerical results, we also include CPU time utilized in the execution of program which is computed by the Mathematica command “TimeUsed[]".  To verify the theoretical order of convergence, we calculate the computational order of convergence (ρc) using the following formula [20]:
(52)ρc=log(||F(x(k))||/||F(x(k-1))||)log(||F(x(k-1))||/||F(x(k-2))||),
taking into consideration the last three approximations in the iterative process.To connect the analysis of computational efficiency with numerical examples, we apply the definition of the computational cost (31), according to which an estimation of the factors μ0 and μ1 is claimed. For this, we express the cost of the evaluation of the elementary functions in terms of products, which depends on the computer, the software, and the arithmetics used (see [21, 22]). In Table 1, an estimation of the cost of the elementary functions in product units is shown, wherein the running time of one product is measured in milliseconds. For the hardware and the software used in the present numerical work, the computational cost of quotient with respect to product is l=2.8 (see Table 1).Table 1
Estimation of computational cost of elementary functions computed with Mathematica7.0 in a processor Intel(R) Core (TM) i5-480M CPU @ 2.67 GHz (64-bit Machine) Microsoft Windows 7 Home Basic 2009, where x=3-1 and y=5.
Digits
x
*
y
l
=
x
/
y
x
exp
(
x
)
ln
(
x
)
sin
(
x
)
cos
(
x
)
arccos
(
x
)
arctan
(
x
)
x
y
4096
0.4665 ms
2.8
1.3
81
78
102
102
166
162
178We consider the following problems for numerical tests.Problem 1.
Consider the system of two equations:(53)x1+ex2-cos(x2)=0,3x1-x2-sin(x2)=0,
with initial approximation x(0)={0.1,0.3}t. The solution is α={0,0}t. The concrete values of parameters in this case are (n,μ0,μ1)=(2,143,71.25), which we use in (38)–(46) for calculating computational costs and efficiency indices of different methods.Problem 2.
Consider the mixed Hammerstein integral equation (see [2, 17]):
(54)x(s)=1+15∫01‍G(s,t)x(t)3dt,
where x∈C[0,1]; s,t∈[0,1] and the kernel G is
(55)G(s,t)={(1-s)t,t⩽s,s(1-t),s⩽t.
We transform the above equation into a finite-dimensional problem by using Gauss-Legendre quadrature formula given as
(56)∫01‍f(t)dt≈∑j=1m‍ϖjf(tj),
where the abscissas tj and the weights ϖj are determined for n=2 by Gauss-Legendre quadrature formula. Denoting the approximation of x(ti) by xi(i=1,2), we obtain the system of nonlinear equations:
(57)xi-1-15∑j=12‍aijxj3=0,
where
(58)aij={ϖjtj(1-ti)ifj⩽i,ϖjti(1-tj)ifi<j.
The initial approximation assumed is x(0)={-1,-1}t and the solution of this problem is α={1.0594598013021009…, 1.0594598013021009…}t. Here the values of (n,μ0,μ1) are (2,5,1.5).Problem 3.
Consider the system of three equations:(59)10x1+sin(x1+x2)-1=0,8x2-cos2(x3-x2)-1=0,12x3+sinx3-1=0,
with initial approximation x(0)={4/5,1/2,1/8}t. The solution correct to 16 places of decimal is α={0.0689783491726666,0.2464424186091830, 0.0769289119875370}t and the concrete values of parameters are (n,μ0,μ1)=(3,103.3,45.56).Problem 4.
Consider the system of eight equations:(60)xi2(xi+xi+1)-1.5=0,i=1,2,…,7x82(x1+x8)-1.5=0.
The initial approximation chosen is x(0)={1.1,1.1,…,1.1}t for the solution α={0.9085602964160698, 0.9085602964160698,…,0.9085602964160698}t. Here the values of (n,μ0,μ1) are (8,2,0.5).Problem 5.
Consider the system of thirteen equations:(61)∑j=1,j≠i13‍xj-e-xi=0,1≤i≤13,
with the initial value x(0)={0.2,0.2,…,0.2}t towards the solution α={0.0771462076130646,0.0771462076130646,…, 0.0771462076130646}t. For this problem the corresponding values of parameters (n,μ0,μ1) are (13,81,6.23).Problem 6.
Consider the system of fifteen equations:(62)xi-cos(2xi-∑j=115‍xj)=0,1≤i≤15,
with initial value x(0)={1.1,1.1,…,1.1}t. The solution of this problem is α={0.9814737799885131,0.9814737799885131,…, 0.9814737799885131}t and (n,μ0,μ1)=(15,102,6.8).Problem 7.
Now consider the system of twenty equations:(63)xi2xi+1-1=0,i=1,2,…,19,x202x1-1=0,
with initial approximation x(0)={1.2,1.2,…,1.2}t and the solution α={1,1,…,1}t. The concrete values of parameters in this case are (n,μ0,μ1)=(20,2,0.15).The numbers of iterations(k), the computational order of convergence (ρc), the computational costs (Ci,p), in terms of products, the computational efficiencies (Ei,p), and the mean CPU time (CPUtime) for each method are displayed in Table 2. Computational cost and efficiency are calculated according to the corresponding expressions given by (38)–(46) by using the values of parameters n, μ0, and μ1 as shown at the end of each problem, while taking l=2.8 in each case. The mean 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-200 as the stopping criterion in single performance of the program.Table 2
Performance of methods.
Methods
k
ρ
c
C
i
,
p
E
i
,
p
CPUtime
Problem1
G1,3
7
3.0000
883.0
1.001245
0.205573
G
2,3
7
3.0000
876.0
1.001255
0.192909
G
3,3
7
3.0000
880.8
1.001248
0.203215
G
4,3
7
3.0000
888.8
1.001237
0.218364
G
5,3
7
5.0000
1468.8
1.000748
0.460909
G
1,5
5
5.0000
1190.2
1.001353
0.176545
G
2,5
5
5.0000
1174.4
1.001370
0.172728
G
1,6
5
6.0000
1205.8
1.001487
0.166418
G
2,6
5
6.0000
1186.4
1.001511
0.162545
Problem2
G
1,3
6
3.0000
49.0
1.022674
0.052327
G
2,3
6
3.0000
45.0
1.024714
0.050854
G
3,3
7
3.0000
46.8
1.023752
0.051545
G
4,3
7
3.0000
54.8
1.020250
0.052173
G
5,3
7
3.0000
76.8
1.014408
0.102091
G
1,5
5
5.0000
80.2
1.020271
0.052734
G
2,5
5
5.0000
64.4
1.025306
0.050727
G
1,6
5
6.0000
95.8
1.018879
0.053545
G
2,6
5
6.0000
76.4
1.023730
0.051182
Problem3
G
1,3
6
3.0000
1187.2
1.000926
0.489455
G
2,3
6
3.0000
1072.0
1.001025
0.398545
G
3,3
6
3.0000
1188.6
1.000925
0.540818
G
4,3
6
3.0000
1203.6
1.000913
0.566636
G
5,3
6
3.0000
2041.7
1.000538
1.023910
G
1,5
5
5.0000
1537.9
1.001047
0.379364
G
2,5
5
5.0000
1512.9
1.001064
0.376345
G
1,6
4
6.0000
1567.3
1.001144
0.365864
G
2,6
4
6.0000
1536.9
1.001167
0.364636
Problem4
G
1,3
7
3.0000
535.2
1.002055
0.340364
G
2,3
7
3.0000
439.2
1.002505
0.228364
G
3,3
7
3.0000
681.6
1.001613
0.370091
G
4,3
7
3.0000
761.6
1.001444
0.387273
G
5,3
7
3.0000
913.6
1.001203
0.393000
G
1,5
5
5.0000
780.0
1.002066
0.237091
G
2,5
5
5.0000
776.0
1.002076
0.295000
G
1,6
5
6.0000
938.4
1.001911
0.262727
G
2,6
5
6.0000
920.0
1.001949
0.240543
Problem5
G
1,3
6
3.0000
4606.9
1.000238
1.660450
G
2,3
6
3.0000
4412.1
1.000249
1.618180
G
3,3
6
3.0000
5293.3
1.000208
2.425090
G
4,3
6
3.0000
5488.3
1.000200
2.437910
G
5,3
6
3.0000
7802.1
1.000141
3.006550
G
1,5
4
5.0000
6226.7
1.000259
1.515910
G
2,5
4
5.0000
6538.7
1.000246
1.799640
G
1,6
4
6.0000
6614.1
1.000271
1.435450
G
2,6
4
6.0000
6902.7
1.000260
1.510730
Problem6
G1,3
8
3.0000
6658.0
1.000165
15.81950
G
2,3
8
3.0000
6403.0
1.000172
12.79912
G
3,3
8
3.0000
7727.0
1.000142
18.29325
G
4,3
8
3.0000
7982.0
1.000138
18.29750
G
5,3
8
3.0000
11312.0
1.000097
28.52323
G
1,5
6
5.0000
8932.0
1.000180
11.58340
G
2,5
6
5.0000
9509.0
1.000169
13.87700
G
1,6
6
6.0000
9439.0
1.000190
11.37760
G
2,6
6
6.0000
9989.0
1.000179
12.69813
Problem7
G
1,3
7
3.0000
4474.0
1.000246
0.697636
G
2,3
7
3.0000
4014.0
1.000274
0.592818
G
3,3
7
3.0000
7056.0
1.000156
0.882182
G
4,3
7
3.0000
7496.0
1.000147
0.975636
G
5,3
7
3.0000
8076.0
1.000136
0.978727
G
1,5
5
5.0000
5806.0
1.000277
0.584010
G
2,5
5
5.0000
7532.0
1.000214
0.753091
G
1,6
5
6.0000
6682.0
1.000268
0.643909
G
2,6
5
6.0000
8372.0
1.000214
0.779636From numerical results, we can observe that like the existing methods the present methods show consistent convergence behavior. It is also clear that the computational order of convergence overwhelmingly supports the theoretical order of convergence. As far as the verification of the results of Theorem5 is concerned, it is simple to check the theoretical results of statement (i) using the numerical values of the efficiency indices displayed in the second last column of Table 2. However, the results of statement (ii) are not quite obvious to verify. In order to do this, we first find m1, m2, m3, and m4 using the values of l, n, and μ1 obtained for each numerical problem and then we compare the efficiencies as per the rules (a)–(d) of statement (ii). The results are displayed in Table 3. From Table 2 we can see that the numerical values of E1,3, E2,3, E1,5, and E1,6 also confirm the results as shown in Table 3.Table 3
Comparison betweenE1,3, E2,3, E1,5, and E1,6  according to Theorem 5(ii).
Problem
μ
0
m
1
m
2
m
3
m
4
Comparison between efficiencies as per
Theorem5(ii)
Case (a)
Case (b)
Case (c)
Case (d)
1
143
146.50
239.60
460.43
-
120.12
E
2,3
>
E
1,3
E
1,5
>
E
1,3
E
1,6
>
E
1,3
E
1,6
>
E
1,5
2
5
7
-
2.87
-
16.52
19.38
E
2,3
>
E
1,3
E
1,3
>
E
1,5
E
1,3
>
E
1,6
E
1,5
>
E
1,6
3
103.3
141.68
228.72
436.50
-
109.76
E
2,3
>
E
1,3
E
1,5
>
E
1,3
E
1,6
>
E
1,3
E
1,6
>
E
1,5
4
2
14
2.95
-
20.19
40.64
E
2,3
>
E
1,3
E
1,5
>
E
1,3
E
1,3
>
E
1,6
E
1,5
>
E
1,6
5
81
95.99
156.09
268.47
-
26.96
E
2,3
>
E
1,3
E
1,5
>
E
1,3
E
1,6
>
E
1,3
E
1,6
>
E
1,5
6
102
119.00
204.39
358.45
-
46.55
E
2,3
>
E
1,3
E
1,5
>
E
1,3
E
1,6
>
E
1,3
E
1,6
>
E
1,5
7
2
25
71.93
85.29
50.17
E
2,3
>
E
1,3
E
1,5
>
E
1,3
E
1,6
>
E
1,3
E
1,5
>
E
1,6Comparison of numerical results shows thatG2,3 is more efficient among third-order methods. However, the present third-order method G1,3 is efficient among the rest of third-order methods in majority of the problems. The present fifth-order G1,5 and sixth-order G1,6 methods are more efficient than the existing methods of same and inferior order for larger systems. This behavior can be observed in the numerical results of Problems 5–7. From the results of last two columns in Table 2, one can conclude that the more is the efficiency of a method the lesser is the computing time of that method. This shows that the efficiency results are in complete agreement with the CPU time utilized in the execution of program.


## 7. Concluding Remarks

In the foregoing study, we have developed iterative methods of third, fifth and sixth-order of convergence for solving systems of nonlinear equations. The computational efficiency in its general form is discussed. Then a comparison between the efficiencies of proposed methods with existing methods is made. It is proved that the present methods are at least competitive with existing methods of similar nature; in particular, these are especially efficient for larger systems. To illustrate the new techniques, seven numerical examples are presented and completely solved. The performance is compared with some known methods of similar character. The theoretical order of convergence and the analysis of computational efficiency are verified in the considered examples. The numerical results have confirmed the robust and efficient nature of the proposed techniques.



---

*Source: 252798-2013-10-31.xml*