Two Efficient Derivative-Free Iterative Methods for Solving Nonlinear Systems

Abstract

In this work, two multi-step derivative-free iterative methods are presented for solving system of nonlinear equations. The new methods have high computational efficiency and low computational cost. The order of convergence of the new methods is proved by a development of an inverse first-order divided difference operator. The computational efficiency is compared with the existing methods. Numerical experiments support the theoretical results. Experimental results show that the new methods remarkably reduce the computing time in the process of high-precision computing.

1. Introduction

Finding the solutions of system of nonlinear equations $F(x)=0$ is a hot problem with wide applications in sciences and engineering, wherein $F:D\subset R^m\to R^m$ and D is an open convex domain in $R^m$. Many efficient methods have been proposed for solving system of nonlinear equations, see for example [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18] and the references therein. The best known method is the Steffensen method [1,2], which is given by

$$y^{(k)}={\psi }_1(x^{(k)},w^{(k)})=x^{(k)}-{[w^{(k)},x^{(k)};F]}^{-1}F(x^{(k)})$$

(1)

where $w^{(k)}=x^{(k)}+F(x^{(k)}),{[w^{(k)},x^{(k)};F]}^{-1}$ is the inverse of $[w^{(k)},x^{(k)};F]$ and $[w^{(k)},x^{(k)};F]:D\subset R^m\to R^m$ is the first order divided difference on D. Equation (1) does not require the derivative of the system F in per iteration.

To reduce the computational time and improve the efficiency index of the Steffensen method, many modified high-order methods have been proposed in open literatures, see [3,4,5,6,7,8,9,10,11,12,13,14] and the references therein. Liu et al. [3] obtained a fourth-order derivative-free method for solving system of nonlinear equations, which can be written as

$$\begin{array}{c}\hfill \left\{\begin{array}{cc}\hfill y^{(k)}=& {\psi }_1(x^{(k)},w^{(k)})\hfill \\ \hfill x^{(k+1)}=& {\psi }_2(x^{(k)},w^{(k)},y^{(k)})\hfill \\ \hfill =& y^{(k)}-{[y^{(k)},x^{(k)};F]}^{-1}([y^{(k)},x^{(k)};F]-[y^{(k)},w^{(k)};F]\hfill \\ \hfill +& [w^{(k)},x^{(k)};F]){[y^{(k)},x^{(k)};F]}^{-1}F(y^{(k)})\hfill \end{array}\right.\end{array}$$

(2)

where $w^{(k)}=x^{(k)}+F(x^{(k)}).$ Grau-Sánchez et al. [4,5] developed some efficient derivative-free methods. One of the methods is the following sixth-order method

$$\begin{array}{c}\hfill \left\{\begin{array}{cc}\hfill y^{(k)}=& x^{(k)}-{[w^{(k)},s^{(k)};F]}^{-1}F(x^{(k)})\hfill \\ \hfill z^{(k)}=& y^{(k)}-{\left\{2[x^{(k)},y^{(k)};F]-[w^{(k)},s^{(k)};F]\right\}}^{-1}F(y^{(k)})\hfill \\ \hfill x^{(k+1)}=& {\psi }_3(w^{(k)},s^{(k)},x^{(k)},y^{(k)},z^{(k)})\hfill \\ \hfill =& z^{(k)}-{\left\{2[x^{(k)},y^{(k)};F]-[w^{(k)},s^{(k)};F]\right\}}^{-1}F(z^{(k)})\hfill \end{array}\right.\end{array}$$

(3)

where $w^{(k)}=x^{(k)}+F(x^{(k)})$ and $s^{(k)}=x^{(k)}-F(x^{(k)}).$ It should be noted that the Equations (2) and (3) need to compute two LU decompositions in per iteration, respectively. Some derivative-free methods are also discussed by Ezquerro et al. in [6] and by Wang et al. in [7,8]. The above multi-step derivative-free iterative methods can save the computing time in the High-precision computing. Therefore, it is meaningful to study the multi-step derivative-free iterative methods.

It is well-known that we can improve the efficiency index of the iterative method and reduce the computational time of the iterative process by reducing the computational cost of the iterative method. There are many ways to reduce the computational cost of the iterative method. In this paper, we reduce the computational cost of the iterative method by reducing the number of LU (lower upper) decompositions in per iteration. Two new derivative-free iterative methods are proposed for solving system of nonlinear equations in Section 2. We prove the local convergence order of the new methods. The feature of the new methods is that the LU decomposition is computed only once in per iteration. Section 3 compares the efficiency of different methods by computational efficiency index [10]. Section 4 illustrates convergence behavior of our methods by numerical examples. Section 5 is a short conclusion.

2. The New Methods and Analysis of Convergence

Using the central difference $[x^{(k)}+F(x^{(k)}),x^{(k)}-F(x^{(k)});F]$, we propose the following iterative scheme

$$\begin{array}{c}\hfill \left\{\begin{array}{cc}\hfill y^{(k)}=& x^{(k)}-{[w^{(k)},s^{(k)};F]}^{-1}F(x^{(k)})\hfill \\ \hfill x^{(k+1)}=& {\psi }_4(x^{(k)},w^{(k)},s^{(k)},y^{(k)})=y^{(k)}-{\mu }_{1}F(y^{(k)})\hfill \end{array}\right.\end{array}$$

(4)

where ${\mu }_1=(3I-2{[w^{(k)},s^{(k)};F]}^{-1}[y^{(k)},x^{(k)};F]){[w^{(k)},s^{(k)};F]}^{-1}$,$w^{(k)}=x^{(k)}+F(x^{(k)}),$$s^{(k)}=x^{(k)}-F(x^{(k)})$ and I is the identity matrix. Furthermore, if we define $z^{(k)}={\psi }_4(x^{(k)},w^{(k)},s^{(k)},y^{(k)})$, then the order of convergence of the following method is six.

$$x^{(k+1)}={\psi }_5(x^{(k)},w^{(k)},s^{(k)},y^{(k)},z^{(k)})=z^{(k)}-{\mu }_{1}F(z^{(k)})$$

(5)

Compared with the Equation (4), the Equation (5) increases one function evaluation $F(z^{(k)})$. In order to simplify calculation, the new Equation (4) can be written as

$$\left\{\begin{array}{c}[w^{(k)},s^{(k)};F]{\gamma }^{(k)}=F(x^{(k)})\hfill \\ y^{(k)}=x^{(k)}-{\gamma }^{(k)}\hfill \\ [w^{(k)},s^{(k)};F]{\delta }_1^{(k)}=F(y^{(k)})\hfill \\ {\delta }_2^{(k)}=[y^{(k)},x^{(k)};F]{\delta }_1^{(k)}\hfill \\ [w^{(k)},s^{(k)};F]{\delta }_3^{(k)}={\delta }_2^{(k)}\hfill \\ x^{(k+1)}=y^{(k)}-3{\delta }_1^{(k)}+2{\delta }_3^{(k)}\hfill \end{array}\right.$$

(6)

Similar strategy can be used in the Equation (5). For the Equations (4) and (5), we have the following analysis of convergence.

Theorem 1. Let $\alpha \in R^m$ be a solution of the system $F(x)=0$ and $F:D\subset R^m\to R^m$ be sufficiently differentiable in an open neighborhood D of α. Then, for an initial approximation sufficiently close to α, the convergence order of iterative Equation (4) is four with the following error equation

$$\epsilon =(4A_2^2-A_3-A_3{F}^{\prime }{(\alpha )}^2)E{e}^2+A_2{E}^2+O(e^5)$$

(7)

where $e=x^{(k)}-\alpha$ and $E=y^{(k)}-\alpha .$ Iterative Equation (5) is of sixth order convergence and satisfies the following error equation

$$e_{n+1}=2A_{2}E\epsilon -(A_3+A_3{F}^{\prime }{(\alpha )}^2-4A_2^2)e^2\epsilon +O(e^7)$$

(8)

where $e_{n+1}=x^{(k+1)}-\alpha$

Proof. The first order divided difference operator of F as a mapping $[·,·;F]:D\times D\subset R^m\times R^m\to L(R^m)$ (see [5,10,11]) is given by

$$[x+h,x;F]={\int }_0^1{F}^{\prime }(x+th)dt,\forall (x,h)\in R^m\times R^m$$

(9)

Expanding $F^{\prime }(x+th)$ in Taylor series at the point x and integrating, we obtain

$${\int }_0^1{F}^{\prime }(x+th)dt=F^{\prime }(x)+\frac{1}{2}{F}^{″}(x)h+\frac{1}{6}{F}^{‴}(x)h^2+O(h^3)$$

(10)

Developing $F(x^{(k)})$ in a neighborhood of α and assuming that $\Gamma ={[F^{\prime }(\alpha )]}^{-1}$ exists, we have

$$F(x^{(k)})=F^{\prime }(\alpha )[e+A_2{e}^2+A_3{e}^3+A_4{e}^4+A_5{e}^5+O(e^6)]$$

(11)

where $A_i=\frac{1}{i!}\Gamma F^{(i)}(\alpha )\in L_i(R^m,R^m).$ The derivatives of $F(x^{(k)})$ can be given by

$$F^{\prime }(x^{(k)})=F^{\prime }(\alpha )[I+2A_{2}e+3A_3{e}^2+4A_4{e}^3+5A_5{e}^4+O(e^5)]$$

(12)

$$F^{″}(x^{(k)})=F^{\prime }(\alpha )[2A_2+6A_{3}e+12A_4{e}^2+20A_5{e}^3+O(e^4)]$$

(13)

$$F^{‴}(x^{(k)})=F^{\prime }(\alpha )[6A_3+24A_{4}e+60A_5{e}^2+O(e^3)]$$

(14)

Setting $y=x+h$ and $E=y-\alpha$, we have $h=E-e$. Replacing the previous expressions Equations (12)–(14) into Equation (10) we get

$$[x^{(k)},y^{(k)};F]=F^{\prime }(\alpha )(I+A_2(E+e)+A_3(E^2+Ee+e^2)+O(e^5))$$

(15)

Noting that $w^{(k)}-\alpha =e+F(x^k)$ and $s^{(k)}-\alpha =e-F(x^k)$, we replace in Equation (15) E by $e+F(x^k)$, e by $e-F(x^k)$, we obtain

$$[w^{(k)},s^{(k)};F]=F^{\prime }(\alpha )(I+2A_{2}e+(3A_3+A_3{F}^{\prime }{(\alpha )}^2)e^2+O(e^3))=F^{\prime }(\alpha )D(e)+O(e^3)$$

(16)

where $D(e)=I+2A_{2}e+(3A_3+A_3{F}^{\prime }{(\alpha )}^2)e^2$ and I is the identity matrix. Using Equation (16), we find

$${[w^{(k)},s^{(k)};F]}^{-1}=D{(e)}^{-1}\Gamma +O(e^3)$$

(17)

Then, we compel the inverse of $D(e)$ to be (see [12,13])

$$D{(e)}^{-1}=I+X_{2}e+X_3{e}^2+O(e^3)$$

(18)

such that $X_2$ and $X_3$ verify

$$D(e)D{(e)}^{-1}=D{(e)}^{-1}D(e)=I$$

(19)

Solving the system Equation (19), we obtain

$$X_2=-2A_2$$

(20)

$$X_3=(4A_2^2-(3A_3+A_3{F}^{\prime }{(\alpha )}^2))$$

(21)

then,

$${[w^{(k)},s^{(k)};F]}^{-1}=(I-2A_{2}e+(4A_2^2-(3A_3+A_3{F}^{\prime }{(\alpha )}^2))e^2+O(e^3))\Gamma$$

(22)

$$E=y^{(k)}-\alpha =e-{[w^{(k)},s^{(k)};F]}^{-1}F(x^{(k)})=A_2{e}^2+O(e^3)$$

(23)

Similar to Equation (11), we have

$$F(y^{(k)})=F^{\prime }(\alpha )[E+A_2{E}^2+O(E^3)]$$

(24)

From Equations (15) and (22)–(24), we get

$$\begin{array}{c}{\mu }_1=(3I-2{[w^{(k)},s^{(k)};F]}^{-1}[y^{(k)},x^{(k)};F]){[w^{(k)},s^{(k)};F]}^{-1}\hfill \\ =(I-2A_{2}E+(A_3+A_3{F}^{\prime }{(\alpha )}^2-4A_2^2)e^2)\Gamma \hfill \end{array}$$

(25)

Taking into account Equation (4), (24) and (25), we obtain

$$\begin{array}{c}\epsilon ={\psi }_4(x^{(k)},w^{(k)},s^{(k)},y^{(k)})-\alpha =E-{\mu }_{1}F(y^{(k)})\hfill \\ =E-(I-2A_{2}E+(A_3+A_3{F}^{\prime }{(\alpha )}^2-4A_2^2)e^2)(E+A_2{E}^2+O(E^3))\hfill \\ =(4A_2^2-A_3-A_3{F}^{\prime }{(\alpha )}^2)E{e}^2+A_2{E}^2+O(e^5)\hfill \end{array}$$

(26)

This means that the Equation (4) is of fourth-order convergence.

Therefore, from Equations (5) and (24)–(26), we obtain the error equation:

$$\begin{array}{c}{e}_{n+1}=x^{(k+1)}-\alpha =\epsilon -{\mu }_{1}F(z^{(k)})\hfill \\ =\epsilon -(I-2A_{2}E+(A_3+A_3{F}^{\prime }{(\alpha )}^2-4A_2^2)e^2)(\epsilon +O({\epsilon }^2))\hfill \\ =2A_{2}E\epsilon -(A_3+A_3{F}^{\prime }{(\alpha )}^2-4A_2^2)e^2\epsilon +O(e^7)\hfill \end{array}$$

(27)

This means that the Equation (5) is of sixth-order convergence. ☐

3. Computational Efficiency

The classical efficiency index $E={\rho }^{1/c}$ (see [9]) is the most used index, but not the only one. We find that the iterative methods with the same classical efficiency index ($E)$ have the different properties in actual applications. The reason is that the number of functional evaluations of iterative method is not the only influence factor in evaluating the efficiency of the iterative method. The number of matrix products, scalar products, decomposition LU of matrix, and the resolution of the triangular linear systems also play an important role in evaluating the real efficiency of iterative method. In this paper, the computational efficiency index ($CEI)$ [10] is used to compare the efficiency of the iterative methods. Some discussions on the $CEI$ can be found in [4,5,6,7]. The $CEI$ of the iterative methods ${\psi }_i(i=1,2,\cdots ,5)$ is given by

$$CE{I}_i(\mu ,m)={\rho }_i^{\frac{1}{C_i(\mu ,m)}},i=1,2,3,4,5$$

(28)

where ${\rho }_i$ is the order of convergence of the method and $C_i(\mu ,m)$ is the computational cost of method. The $C_i(\mu ,m)$ is given by

$$C_i(\mu ,m)=a_i(m)\mu +p_i(m)$$

(29)

where $a_i(m)$ denotes the number of evaluations of scalar functions used in the evaluations of F and $[x,y;F]$, and $p_i(m)$ represents the operational cost per iteration. To express the value of Equation (29) in terms of products, a ratio $\mu >0$ in Equation (29) between products (and divisions) and evaluations of functions is required, see [5,10]. We must add m products for multiplication of a vector by a scalar and $m^2$ products for matrix-vector multiplication. To compute an inverse linear operator, we need $(m^3-m)/3$ products and divisions in the LU decomposition and $m^2$ products and divisions for solving two triangular linear systems. If we compute the first-order divided difference then we need $m(m-1)$ scalar functional evaluations and $m^2$ quotients. The first-order divided difference $[x,y;F]$ of F is given by

$${[y,x;F]}_{ij}=(F_i(y_1\cdots ,y_{j-1},y_j,x_{j+1},\cdots ,x_m)-F_i(y_1\cdots ,y_{j-1},x_j,x_{j+1},\cdots ,x_m))/(y_j-x_j)$$

where $1\le i,j\le m$, $x=(x_1,\cdots x_{j-1},x_j,x_{j+1},\cdots x_m)$ and $y=(y_1,\cdots y_{j-1},y_j,y_{j+1},\cdots y_m)$ (see [9]). Based on Equations (28) and (29), Table 1 shows the computational cost of different methods.

Table 1. Computational cost of the iterative methods.

Methods	ρ	a(m)	p(m)	C(μ, m)
${\psi }_1$	2	$m(m+1)$	$(m^3-m)/3+2m^2$	$C_1=m(m+1)\mu +(m^3-m)/3+2m^2$
${\psi }_2$	4	$3m^2$	$2(m^3-m)/3+7m^2$	$C_2=3m^2\mu +2(m^3-m)/3+7m^2$
${\psi }_3$	6	$m(2m+3)$	$2(m^3-m)/3+6m^2$	$C_3=m(2m+3)\mu +2(m^3-m)/3+6m^2$
${\psi }_4$	4	$2m(m+1)$	$(m^3-m)/3+6m^2+2m$	$C_4=2m(m+1)\mu +(m^3-m)/3+6m^2+2m$
${\psi }_5$	6	$m(2m+3)$	$(m^3-m)/3+9m^2+4m$	$C_5=m(2m+3)\mu +(m^3-m)/3+9m^2+4m$

Table 1. Computational cost of the iterative methods.

Methods	ρ	a(m)	p(m)	C(μ, m)
${\psi }_1$	2	$m(m+1)$	$(m^3-m)/3+2m^2$	$C_1=m(m+1)\mu +(m^3-m)/3+2m^2$
${\psi }_2$	4	$3m^2$	$2(m^3-m)/3+7m^2$	$C_2=3m^2\mu +2(m^3-m)/3+7m^2$
${\psi }_3$	6	$m(2m+3)$	$2(m^3-m)/3+6m^2$	$C_3=m(2m+3)\mu +2(m^3-m)/3+6m^2$
${\psi }_4$	4	$2m(m+1)$	$(m^3-m)/3+6m^2+2m$	$C_4=2m(m+1)\mu +(m^3-m)/3+6m^2+2m$
${\psi }_5$	6	$m(2m+3)$	$(m^3-m)/3+9m^2+4m$	$C_5=m(2m+3)\mu +(m^3-m)/3+9m^2+4m$

From Table 1, we can see that our methods ${\psi }_i(i=4,5)$ need less number of LU decomposition than methods ${\psi }_2$ and ${\psi }_3$. The computational cost of the fourth-order methods show the following order:

$$C_4<C_2,\mathrm{for} m\ge 2$$

(30)

We use the following expressions [10] to compare the CEI of different methods

$$R_{i,j}=\frac{lnCE{I}_i}{lnCE{I}_j}=\frac{ln({\rho }_i)C_j(\mu ,m)}{ln({\rho }_j)C_i(\mu ,m)},i,j=1,2,3,4,5$$

(31)

For $R_{i,j}>1$ the iterative method $M_i$ is more efficient than $M_j.$

Using the $CEI$ of the iterative methods, we obtain the following theorem:

Theorem 2. 1. For the fourth-order method, we have $CE{I}_4>CE{I}_2$ for all $m\ge 2 and \mu >0.$

2. For the sixth-order method, we have $CE{I}_5>CE{I}_3$ for all $m\ge 11 and \mu >0.$

Proof. 1. From Table 1, we note that the methods ${\psi }_i(i=2,4)$ have the same order ${\rho }_2={\rho }_4=4.$ Based on Equations (29) and (30), we get that $CE{I}_4>CE{I}_2$ for all $m\ge 2$ and $\mu >0.$

2. The methods ${\psi }_i(i=3,5)$ have the same order and the same functional evaluations. The relation between ${\psi }_5$ and ${\psi }_3$ can be given by

$$R_{5,3}=\frac{ln({\rho }_5)C_3(\mu ,m)}{ln({\rho }_3)C_5(\mu ,m)}=\frac{m(2m+3)\mu +2(m^3-m)/3+6m^2}{m(2m+3)\mu +(m^3-m)/3+9m^2+4m}$$

(32)

Subtracting the denominator from the numerator of Equation (32), we have

$$\frac{1}{3}m(m^2-9m-13)$$

(33)

The Equation (33) is positive for $m\ge 10.2662.$ Thus, we obtain that $CE{I}_5>CE{I}_3$ for all $m\ge 11$ and $\mu >0.$ ☐

Then, we compare the $CEI$ of the iterative methods with different convergence order by the following theorem:

Theorem 3. We have 1. $CE{I}_5>CE{I}_4$ for all $m\ge 2$ and $\mu >\frac{m^2{ln}^{2/3}+18m{ln}^{4/3}+(17{ln}^2-5{ln}^3)}{6(m{ln}^{3/2}+{ln}^{3/4})}.$

2. $CE{I}_4>CE{I}_1$ for all $m\ge 8$ and $\mu >0.$

Proof. 1. From the expression Equation (31) and Table 1,We get the following relation between ${\psi }_4$ and ${\psi }_5$

$$R_{5,4}=\frac{ln({\rho }_5)C_4(\mu ,m)}{ln({\rho }_4)C_5(\mu ,m)}=\frac{{ln}^6}{{ln}^4}\frac{2m(m+1)\mu +(m^3-m)/3+6m^2+2m}{m(2m+3)\mu +(m^3-m)/3+9m^2+4m}$$

(34)

We consider the boundary $R_{5,4}=1$. The boundary can is given by the following equation

$$\mu =H_{5,4}(m)=\frac{m^2{ln}^{2/3}+18m{ln}^{4/3}+17{ln}^2-5{ln}^3}{6(m{ln}^{3/2}+{ln}^{3/4})}$$

(35)

where $CE{I}_5>CE{I}_4$ over it (see Figure 1). The boundary Equation (35) cut axes at points $(m,\mu )=(13.888,0)$ and $(2,4.7859)$. Thus, we get that $CE{I}_5>CE{I}_4$ since $R_{5,4}>1$ for all $m\ge 2$ and $\mu >H_{5,4}(m).$

2. The relation between ${\psi }_1$ and ${\psi }_4$ is given by

$$R_{4,1}=\frac{ln({\rho }_4)C_1(\mu ,m)}{ln({\rho }_1)C_4(\mu ,m)}=\frac{{ln}^4}{{ln}^2}\frac{m(m+1)\mu +(m^3-m)/3+2m^2}{2m(m+1)\mu +(m^3-m)/3+6m^2+2m}$$

(36)

Subtracting the denominator from the numerator of Equation (36), we have

$$\frac{1}{3}m(m^2-6m-7)$$

(37)

The Equation (37) is positive for $m>7.$ Thus, we obtain that $CE{I}_4>CE{I}_1$ for all $m\ge 8$ and $\mu >0.$ ☐

Figure 1. The boundary function ${\text{H}}_{5,4}$ in$(m,\mu )-$plain.

Figure 1. The boundary function ${\text{H}}_{5,4}$ in$(m,\mu )-$plain.

4. Numerical Examples

In this section, we compare the performance of related methods by mathematical experiments. The numerical experiments have been carried out using Maple 14 computer algebra system with 2048 digits. The computer specifications are Microsoft Windows 7 Intel(R), Core(TM) i3-2350M CPU, 1.79 GHz with 2 GB of RAM.

According to the Equation (29), the factor μ is claimed by expressing the cost of the evaluation of elementary functions in terms of products [15]. Table 2 gives an estimation of the cost of the elementary functions in amount of equivalent products, where the running time of one product is measured in milliseconds.

Table 2. Estimation of computational cost of elementary functions computed with Maple 14 and using a processor Intel^® Core(TM) i3-2350M CPU, 1.79 GHz (32-bit Machine) Microsoft Windows 7 Professional, where $x=\sqrt{3}-1$ and $y=\sqrt{5}$.

Digits	x·y	x/y	$\sqrt{\mathit{x}}$	$\mathbf{exp}(\mathit{x})$	$\mathbf{ln}(\mathit{x})$	$\mathbf{sin}(\mathit{x})$	$\mathbf{cos}(\mathit{x})$	$\mathbf{arctan}(\mathit{x})$
2048	0.109 ms	1	5	53	12	112	110	95

Table 2. Estimation of computational cost of elementary functions computed with Maple 14 and using a processor Intel^® Core(TM) i3-2350M CPU, 1.79 GHz (32-bit Machine) Microsoft Windows 7 Professional, where $x=\sqrt{3}-1$ and $y=\sqrt{5}$.

Digits	x·y	x/y	$\sqrt{\mathit{x}}$	$\mathbf{exp}(\mathit{x})$	$\mathbf{ln}(\mathit{x})$	$\mathbf{sin}(\mathit{x})$	$\mathbf{cos}(\mathit{x})$	$\mathbf{arctan}(\mathit{x})$
2048	0.109 ms	1	5	53	12	112	110	95

Table 3, Table 4, Table 5, Table 6, Table 7 and Table 8 show the following information of the methods ${\psi }_i(i=1,2,\cdots ,5)$: the number of iterations k needed to converge to the solution, the norm of function $F(x^{(k)})$at the last step, the value of the stopping factors at the last step, the computational cost C, the computational time $Time(s)$, the computational efficiency indices $CEI$ and the computational order of convergence ρ. Using the commond time( ) in Maple 14, we can obtain the computational time of different methods. The computational order of convergence ρ is defined by [16]:

$$\rho \approx \frac{ln(||x^{(k+1)}-x^{(k)}||/||x^{(k)}-x^{(k-1)}||)}{ln(||x^{(k)}-x^{(k-1)}||/||x^{(k-1)}-x^{(k-2)}||)}$$

(38)

The following problems are chosen for numerical tests:

Example 1 Considering the following system

$$\left\{\begin{array}{c}{x}_1+e^{x_1}-cos(x_2)=0\hfill \\ 3x_1-x_2-sin(x_2)=0\hfill \end{array}\right.$$

where $(m,\mu )=(2,\frac{53+110+112+1}{2})=(2,138)$ are the values used in Equation (29). $x^{(0)}={(0.5,0.5)}^T$ is the initial point and $\alpha \approx {(0,0)}^T$ is the solution of the Example 1. $\left|\right|x^{(k)}-x^{(k-1)}\left|\right|<{10}^{-100}$ is the stopping criterion.

The results shown in Table 3 confirm the first assertion of Theorem 2 and the first assertion of Theorem 3 for $m\doteq 2.$ Namely, $CE{I}_5>CE{I}_4$ for $\mu >4.7859$. The new sixth-order method ${\psi }_5$ spends minimum time for finding the numerical solution. The ’nc’ denotes that the method does not converge in the Table 3.

Table 3. Performance of methods for Example 1.

Method	k	$\left|\right|{\mathbf{x}}^{(\mathbf{k})}-{\mathbf{x}}^{(\mathbf{k}\mathbf{-}\mathbf{1})}\left|\right|$	$\left|\right|\mathbf{F}({\mathbf{x}}^{(\mathbf{k})})\left|\right|$	ρ	C	CEI	Time(s)
${\psi }_1$	13	1.792e−161	3.748e−322	2.00000	838	1.0008275	1.127
${\psi }_2$	nc
${\psi }_3$	4	3.558e−743	8.245e−496	6.00314	1960	1.0009148	0.780
${\psi }_4$	5	4.086e−211	5.330e−421	4.00015	1686	1.0008226	0.836
${\psi }_5$	4	6.240e−164	2.389e−489	6.00420	1978	1.0009063	0.546

Table 3. Performance of methods for Example 1.

Method	k	$\left|\right|{\mathbf{x}}^{(\mathbf{k})}-{\mathbf{x}}^{(\mathbf{k}\mathbf{-}\mathbf{1})}\left|\right|$	$\left|\right|\mathbf{F}({\mathbf{x}}^{(\mathbf{k})})\left|\right|$	ρ	C	CEI	Time(s)
${\psi }_1$	13	1.792e−161	3.748e−322	2.00000	838	1.0008275	1.127
${\psi }_2$	nc
${\psi }_3$	4	3.558e−743	8.245e−496	6.00314	1960	1.0009148	0.780
${\psi }_4$	5	4.086e−211	5.330e−421	4.00015	1686	1.0008226	0.836
${\psi }_5$	4	6.240e−164	2.389e−489	6.00420	1978	1.0009063	0.546

Example 2 The second system is defined by [11]

$$\left\{\begin{array}{c}{x}_2+x_3-e^{-x_1}=0\hfill \\ x_1+x_3-e^{-x_3}=0\hfill \\ x_1+x_2-e^{-x_3}=0\hfill \end{array}\right.$$

where $(m,\mu )=(3,\frac{53+53}{3})=(3,35.3)$. The initial point is $x^{(0)}=(0.5,0.5,0.5)$. $\left|\right|x^{(k)}-x^{(k-1)}\left|\right|<{10}^{-200}$ is the stopping criterion.The solution is $\alpha \approx (0.3517337,0.3517337,0.3517337).$

The results shown in Table 4 confirm the first assertion of Theorem 2 and assertion 1 of Theorem 3 for $m=3.$ Namely, $CE{I}_4>CE{I}_2$ and $CE{I}_5>CE{I}_4$ for $\mu >4.7859$. Table 4 shows that sixth-order method ${\psi }_5$ is the most efficient iterative method in both computational time and $CEI$.

Table 4. Performance of methods for Example 2.

Method	k	$\left|\right|{\mathbf{x}}^{(\mathbf{k})}-{\mathbf{x}}^{(\mathbf{k}\mathbf{-}\mathbf{1})}\left|\right|$	$\left|\right|\mathbf{F}({\mathbf{x}}^{(\mathbf{k})})\left|\right|$	ρ	C	CEI	Time(s)
${\psi }_1$	9	2.136e−302	5.945e−604	2	449.6	1.00154289	0.514
${\psi }_2$	5	2.439e−675	2.703e−1350	4	1032.1	1.00134408	0.592
${\psi }_3$	4	1.414e−1080	7.020e−1620	6	1023.1	1.00175284	0.561
${\psi }_4$	5	4.123e−699	9.73957e−1397	4	915.2	1.00151589	0.561
${\psi }_5$	4	9.097e−550	8.57708e−1647	6	1054.1	1.00170125	0.483

Table 4. Performance of methods for Example 2.

Method	k	$\left|\right|{\mathbf{x}}^{(\mathbf{k})}-{\mathbf{x}}^{(\mathbf{k}\mathbf{-}\mathbf{1})}\left|\right|$	$\left|\right|\mathbf{F}({\mathbf{x}}^{(\mathbf{k})})\left|\right|$	ρ	C	CEI	Time(s)
${\psi }_1$	9	2.136e−302	5.945e−604	2	449.6	1.00154289	0.514
${\psi }_2$	5	2.439e−675	2.703e−1350	4	1032.1	1.00134408	0.592
${\psi }_3$	4	1.414e−1080	7.020e−1620	6	1023.1	1.00175284	0.561
${\psi }_4$	5	4.123e−699	9.73957e−1397	4	915.2	1.00151589	0.561
${\psi }_5$	4	9.097e−550	8.57708e−1647	6	1054.1	1.00170125	0.483

Example 3 Now, considering the following large scale nonlinear systems [17]:

$$\left\{\begin{array}{c}{x}_i{x}_{i+1}-1=0,1\le i\le m-1\hfill \\ x_m{x}_1-1=0\hfill \end{array}\right.$$

The initial vector is $x^{(0)}={\{1.5,1.5,\cdots ,1.5\}}^t$ for the solution $\alpha ={\{1,1,\cdots ,1\}}^t.$ The stopping criterion is $\left|\right|x^{(k)}-x^{(k-1)}\left|\right|<{10}^{-100}$.

Table 5. Performance of methods for Example 3, where $(m,\mu )=(199,1)$.

Method	k	$\left|\right|{\mathbf{x}}^{(\mathbf{k})}-{\mathbf{x}}^{(\mathbf{k}\mathbf{-}\mathbf{1})}\left|\right|$	$\left|\right|\mathbf{F}({\mathbf{x}}^{(\mathbf{k})})\left|\right|$	ρ	C	CEI	Time(s)
${\psi }_1$	10	4.993e−150	7.480e−299	2.00000	2,745,802	1.000000252	95.940
${\psi }_2$	5	3.013e−212	1.210e−423	4.00000	5,649,610	1.000000245	126.438
${\psi }_3$	4	3.922e−556	2.197e−833	5.99998	5,571,005	1.000000322	77.111
${\psi }_4$	5	1.404e−269	9.850e−538	4.00000	2,944,404	1.000000471	81.042
${\psi }_5$	4	5.298e−208	2.231e−621	5.99976	3,063,804	1.000000585	64.818

Table 5. Performance of methods for Example 3, where $(m,\mu )=(199,1)$.

Method	k	$\left|\right|{\mathbf{x}}^{(\mathbf{k})}-{\mathbf{x}}^{(\mathbf{k}\mathbf{-}\mathbf{1})}\left|\right|$	$\left|\right|\mathbf{F}({\mathbf{x}}^{(\mathbf{k})})\left|\right|$	ρ	C	CEI	Time(s)
${\psi }_1$	10	4.993e−150	7.480e−299	2.00000	2,745,802	1.000000252	95.940
${\psi }_2$	5	3.013e−212	1.210e−423	4.00000	5,649,610	1.000000245	126.438
${\psi }_3$	4	3.922e−556	2.197e−833	5.99998	5,571,005	1.000000322	77.111
${\psi }_4$	5	1.404e−269	9.850e−538	4.00000	2,944,404	1.000000471	81.042
${\psi }_5$	4	5.298e−208	2.231e−621	5.99976	3,063,804	1.000000585	64.818

Table 6. The computational time (in second) for Example 3 by the methods.

Method	${\mathit{\psi }}_{\mathbf{1}}$	${\mathit{\psi }}_{\mathbf{2}}$	${\mathit{\psi }}_{\mathbf{3}}$	${\mathit{\psi }}_{\mathbf{4}}$	${\mathit{\psi }}_{\mathbf{5}}$
$m=99$	20.982	29.499	16.848	19.219	15.459
$m=199$	95.940	126.438	77.111	81.042	64.818
$m=299$	254.234	328.896	207.340	199.930	156.094

Table 6. The computational time (in second) for Example 3 by the methods.

Method	${\mathit{\psi }}_{\mathbf{1}}$	${\mathit{\psi }}_{\mathbf{2}}$	${\mathit{\psi }}_{\mathbf{3}}$	${\mathit{\psi }}_{\mathbf{4}}$	${\mathit{\psi }}_{\mathbf{5}}$
$m=99$	20.982	29.499	16.848	19.219	15.459
$m=199$	95.940	126.438	77.111	81.042	64.818
$m=299$	254.234	328.896	207.340	199.930	156.094

Application in Integral Equations

The Chandrasekhar integral [18] equation comes from radiative transfer theory, which is given by

$$F(P,c)=0,P:[0,1]\to R$$

with the operator F and parameter c as

$$F(P,c)(u)=P(u)-{\left(1-\frac{c}{2}{\int }_0^1\frac{uP(v)}{u+v}dv\right)}^{-1}$$

We approximate the integrals by the composite midpoint rule:

$${\int }_0^{1}f(t)dt=\frac{1}{m}\sum _{j=1}^{m}f(t_j)$$

where $t_j=(j-1/2)/m$ for $1\le j\le m.$ We obtain the resulting discrete problem is

$$F_i(u,c)=u_i-{\left(1-\frac{c}{2m}\sum _{j=1}^m\frac{t_i{u}_j}{t_i+t_j}\right)}^{-1},1\le i\le m$$

The initial vector is $x^{(0)}={\left\{1.5,1.5,\cdots ,1.5\right\}}^t$, $c=0.9$. Table 7 and Table 8 show the numerical results of this problem. $\left|\right|F(x^{(k)})\left|\right|<{10}^{-200}$ is the stopping criterion of this problem.

Table 7. The computational time (in second) for solving Chandrasekhar Integral equation.

Method	${\mathit{\psi }}_{\mathbf{1}}$	${\mathit{\psi }}_{\mathbf{2}}$	${\mathit{\psi }}_{\mathbf{3}}$	${\mathit{\psi }}_{\mathbf{4}}$	${\mathit{\psi }}_{\mathbf{5}}$
$m=30$	88.468	207.200	87.937	102.055	70.309
$m=60$	422.388	904.602	435.929	488.969	400.345

Table 7. The computational time (in second) for solving Chandrasekhar Integral equation.

Method	${\mathit{\psi }}_{\mathbf{1}}$	${\mathit{\psi }}_{\mathbf{2}}$	${\mathit{\psi }}_{\mathbf{3}}$	${\mathit{\psi }}_{\mathbf{4}}$	${\mathit{\psi }}_{\mathbf{5}}$
$m=30$	88.468	207.200	87.937	102.055	70.309
$m=60$	422.388	904.602	435.929	488.969	400.345

Table 8. The number of iterations for solving Chandrasekhar Integral equation.

Method	${\mathit{\psi }}_{\mathbf{1}}$	${\mathit{\psi }}_{\mathbf{2}}$	${\mathit{\psi }}_{\mathbf{3}}$	${\mathit{\psi }}_{\mathbf{4}}$	${\mathit{\psi }}_{\mathbf{5}}$
$m=30$	8	6	4	5	4
$m=60$	8	6	4	5	4

Table 8. The number of iterations for solving Chandrasekhar Integral equation.

Method	${\mathit{\psi }}_{\mathbf{1}}$	${\mathit{\psi }}_{\mathbf{2}}$	${\mathit{\psi }}_{\mathbf{3}}$	${\mathit{\psi }}_{\mathbf{4}}$	${\mathit{\psi }}_{\mathbf{5}}$
$m=30$	8	6	4	5	4
$m=60$	8	6	4	5	4

The results shown in Table 5 confirm the assertion of Theorem 2 and Theorem 3 for $m=199.$ Namely, $CE{I}_4>CE{I}_2,CE{I}_5>CE{I}_3,$ and $CE{I}_4>CE{I}_1$. From the Table 6, we remark that the computational time of our fourth-order method ${\psi }_4$ is less than that of the sixth order method ${\psi }_3$ for $m=299$. Table 5, Table 6 and Table 7 show that, as the nonlinear system is big-sized, our new methods ${\psi }_i(i=4,5)$ remarkably reduce the computational time.

The numerical results shown in Table 3, Table 4, Table 5, Table 6, Table 7 and Table 8 are in concordance with the theory developed in this paper. The new methods require less number of iterations to obtain higher accuracy in the contrast to the other methods. The most important is that our methods have higher $CEI$ and lower computational time than other methods in this paper. The sixth-order method ${\psi }_5$ is the most efficient iterative methods in both $CEI$ and computational time.

5. Conclusions

In this paper, two high-order iterative methods for solving system of nonlinear equations are obtained. The new methods are derivative free. The order of convergence of the new methods is proved by using a development of an inverse first-order divided difference operator. Moreover, the computational efficiency index for system of nonlinear equations is used to compare the efficiency of different methods. Numerical experiments show that our methods remarkably reduce the computational time for solving big-sized system of nonlinear equations. The main reason is that the LU decomposition of the matrix of our methods is computed only once in per iteration. We concluded that, in order to obtain an efficient iterative method, we should comprehensively consider the number of functional evaluations, the convergence order and the operational cost of the iterative.