Two-Step Fifth-Order Efficient Jacobian-Free Iterative Method for Solving Nonlinear Systems

1. Introduction

Let $F\left(x\right)=0$ be a nonlinear system of equations, $F:D\subset {\mathbb{R}}^n\to {\mathbb{R}}^n$, and the functions $f_i$ for $i=1,2,\dots ,n$, are the coordinate components of F, expressed as $F\left(x\right)={\left(f_1\left(x\right),f_2\left(x\right),\dots ,f_n\left(x\right)\right)}^T$. Solving nonlinear systems is generally challenging, and solutions $\xi$ are typically found by linearizing the problem or employing a fixed-point iteration function $G:D\subset {\mathbb{R}}^n\to {\mathbb{R}}^n$, leading to an iterative fixed-point method. Among the various root-finding techniques for nonlinear systems, Newton’s method is the most well-known, which follows the second-order iterative procedure:

$$x^{(k+1)}=x^{\left(k\right)}-{\left[F^{\prime }\left(x^{\left(k\right)}\right)\right]}^{-1}F\left(x^{\left(k\right)}\right),k=0,1,\dots ,$$

being $F^{\prime }\left(x^{\left(k\right)}\right)$ the Jacobian matrix of F at the k-th iterate.

Recenty, many researchers have focused on developing iterative methods that outperform Newton’s method in terms of both efficiency and order of convergence. Numerous approaches need the computation of $F^{\prime }$ at different points along each iteration. Nevertheless, calculating the Jacobian poses significant challenges, particularly in high-dimensional problems, where its computation can be costly or even impractical. In some instances, the Jacobian may not exist at all.

To address this issue, alternative approaches have been proposed, such as replacing the Jacobian matrix with a divided difference operator. One of the simplest alternatives is the multidimensional version of Steffensen’s method, attributed to Samanskii [1,2], which substitutes the Jacobian in Newton’s procedure with a first-order operator of divided differences:

$$x^{(k+1)}=x^{\left(k\right)}-{[x^{\left(k\right)},z^{\left(k\right)};F]}^{-1}F\left(x^{\left(k\right)}\right),k=0,1,\dots ,$$

being $z^{\left(k\right)}=x^{\left(k\right)}+F\left(x^{\left(k\right)}\right)$, and $[·,·;F]:\Omega \times \Omega \subset {\mathbb{R}}^n\times {\mathbb{R}}^n\to \mathcal{L}\left({\mathbb{R}}^n\right)$ is the operator of divided differences related to F [3],

$$[y,x;F](y-x)=F\left(y\right)-F\left(x\right) \mathrm{for} \mathrm{any} x,y\in \Omega .$$

This substitution retains the 2-nd order of convergence while bypassing the calculation of the $F^{\prime }$.

Although both Steffensen and Newton methods exhibit quadratic convergence, it has been shown that Steffensen’s scheme is less stable than Newton’s method, with stability depending more on the initial guess. This has been thoroughly analyzed in [4,5], where it was found that, for scalar cases $f\left(x\right)=0$, derivative-free iterative methods become more stable when selecting $z=x+\alpha f\left(x\right)$ for small real values of $\alpha$.

However, substituting the Jacobian with divided differences can result in lower convergence order for some iterative methods. For example, the multidimensional version of Ostrowski’s fourth-order method (see [6,13]):

$$\begin{array}{cc}\hfill y^{\left(k\right)}& =x^{\left(k\right)}-{\left[F^{\prime }\left(x^{\left(k\right)}\right)\right]}^{-1}F\left(x^{\left(k\right)}\right),\hfill \\ \hfill x^{(k+1)}& =y^{\left(k\right)}-{\left[2[y^{\left(k\right)},x^{\left(k\right)};F]-F^{\prime }\left(x^{\left(k\right)}\right)\right]}^{-1}F\left(y^{\left(k\right)}\right),k=0,1,\dots ,\hfill \end{array}$$

achieves only cubic convergence if $F^{\prime }\left(x\right)$ is replaced by $[y,x;F]$, as follows:

$$\begin{array}{cc}\hfill y^{\left(k\right)}& =x^{\left(k\right)}-{[x^{\left(k\right)}+F\left(x^{\left(k\right)}\right),x^{\left(k\right)};F]}^{-1}F\left(x^{\left(k\right)}\right),k=0,1,\dots ,\hfill \\ \hfill x^{(k+1)}& =y^{\left(k\right)}-{\left[2[y^{\left(k\right)},x^{\left(k\right)};F]-[x^{\left(k\right)}+F\left(x^{\left(k\right)}\right),x^{\left(k\right)};F]\right]}^{-1}F\left(y^{\left(k\right)}\right).\hfill \end{array}$$

Other fourth-order methods also loose their convergence order when the Jacobian is replaced with divided differences, such as Jarratt’s scheme [7], Sharma’s method [8], Montazeri’s method [9], Ostrowski’s vectorial extension ([13,15]), and Sharma-Arora’s fifth-order scheme [10]. In all these cases, Jacobian-free versions of the methods reduce to lower orders of convergence.

Nevertheless, Amiri et al. [11] demonstrated that using a specialized divided difference operator of the form $[x,x+G(x);F]$, where $G\left(x\right)={\left(f_1{\left(x\right)}^m,f_2{\left(x\right)}^m,\dots ,f_n{\left(x\right)}^m\right)}^T$, $m\in \mathbb{N}$, as an approximation of the Jacobian matrix, may preserve the convergence order. By selecting an appropriate parameter m, the original fourth-order convergence of these methods can be maintained.

Despite the reduction in performance observed in some Jacobian-free methods, it is important to highlight that there are iterative methods that are successfully modified in their iterative expressions to preserve the order of convergence, even after fully transitioning to Jacobian-free formulations. This is the case of a combination of the Traub-Steffensen family of methods and a second step with divided differences operators, proposed by Behl et al. in [12]

$$\begin{array}{cc}\hfill y^{\left(k\right)}& =x^{\left(k\right)}-{\left[u^{\left(k\right)},x^{\left(k\right)};F\right]}^{-1}F\left(x^{\left(k\right)}\right),k=0,1,2,\dots \hfill \\ \hfill x^{(k+1)}& =y^{\left(k\right)}-{\left[y^{\left(k\right)},x^{\left(k\right)};F\right]}^{-1}\left[u^{\left(k\right)},x^{\left(k\right)};F\right]{\left[u^{\left(k\right)},y^{\left(k\right)};F\right]}^{-1}F\left(y^{\left(k\right)}\right),\hfill \end{array}$$

(1)

where $u^{\left(k\right)}=x^{\left(k\right)}+\beta F\left(x^{\left(k\right)}\right),\beta \in \mathbb{R}$. The iterative schemes have a fourth order of convergence for every $\beta$, $\beta \ne 0$, for our purposes we choose $\beta =1$ and this be called $\mathbf{Traub}-\mathbf{Ste}$.

Now, we consider several efficient vectorial iterative schemes existing in the literature to transform them in their Jacobian-free versions following the idea of Amiri et al. [11]. In the following sections, we compare these later schemes with our proposed procedures, in terms of efficiency and numerical performance. The first one is the vectorial extension of Ostrowski’s scheme (see [13–15], for instance),

$$\begin{array}{cc}\hfill y^{\left(k\right)}& =x^{\left(k\right)}-{\left[F^{\prime }\left(x^{\left(k\right)}\right)\right]}^{-1}F\left(x^{\left(k\right)}\right),\hfill \\ \hfill x^{(k+1)}& =y^{\left(k\right)}-{\left[2\left[y^{\left(k\right)},x^{\left(k\right)};F\right]-F^{\prime }\left(x^{\left(k\right)}\right)\right]}^{-1}F\left(y^{\left(k\right)}\right),\hfill \end{array}$$

(2)

whose Jacobian-free version obtained by substituting the Jacobian matrix by the divided difference operator (with Amiri et al. approach [11], $m=2$) is

$$\begin{array}{cc}\hfill y^{\left(k\right)}& =x^{\left(k\right)}-{\left[u^{\left(k\right)},x^{\left(k\right)};F\right]}^{-1}F\left(x^{\left(k\right)}\right),\hfill \\ \hfill x^{(k+1)}& =y^{\left(k\right)}-{\left[2\left[y^{\left(k\right)},x^{\left(k\right)};F\right]-\left[u^{\left(k\right)},x^{\left(k\right)};F\right]\right]}^{-1}F\left(y^{\left(k\right)}\right),\hfill \end{array}$$

(3)

where $\left[u^{\left(k\right)},x^{\left(k\right)},F\right]\approx F^{\prime }\left(x^{\left(k\right)}\right)$, and $u^{\left(k\right)}=x^{\left(k\right)}+\alpha G\left(x^{\left(k\right)}\right)$ (with $m=2$), $\alpha \in \mathbb{R}$. We denote this method as ${\mathbf{Ostro}}_{\mathbf{01}}$.

Another fourth-order method proposed by Sharma in [16] using Jacobian matrices is

$$\begin{array}{cc}\hfill y^{\left(k\right)}& =x^{\left(k\right)}-{\left[F^{\prime }\left(x^{\left(k\right)}\right)\right]}^{-1}F\left(x^{\left(k\right)}\right),\hfill \\ \hfill x^{(k+1)}& =y^{\left(k\right)}-\left[3I-2{\left[F^{\prime }\left(x^{\left(k\right)}\right)\right]}^{-1}[y^{\left(k\right)},x^{\left(k\right)};F]\right]{\left[F^{\prime }\left(x^{\left(k\right)}\right)\right]}^{-1}F\left(y^{\left(k\right)}\right),\hfill \end{array}$$

(4)

to which we apply the same Amiri’s procedure performed for (2), getting its Jacobian-free partner

$$\begin{array}{cc}\hfill y^{\left(k\right)}& =x^{\left(k\right)}-{\left[u^{\left(k\right)},x^{\left(k\right)};F\right]}^{-1}F\left(x^{\left(k\right)}\right),\hfill \\ \hfill x^{(k+1)}& =y^{\left(k\right)}-\left[3I-2{\left[u^{\left(k\right)},x^{\left(k\right)};F\right]}^{-1}[y^{\left(k\right)},x^{\left(k\right)};F]\right]{\left[u^{\left(k\right)},x^{\left(k\right)};F\right]}^{-1}F\left(y^{\left(k\right)}\right),\hfill \end{array}$$

(5)

that we denote by ${\mathbf{M}}_{\mathbf{4},\mathbf{3}}$, where $\left[u^{\left(k\right)},x^{\left(k\right)},F\right]\approx F^{\prime }\left(x^{\left(k\right)}\right)$ and $u^{\left(k\right)}=x^{\left(k\right)}+\alpha G\left(x^{\left(k\right)}\right)$, $m=2$, firstly appeared in [11].

We finish with a sixth-order scheme [16], which is obtained by adding a step to the previous method (4),

$$\begin{array}{cc}\hfill y^{\left(k\right)}& =x^{\left(k\right)}-{\left[F^{\prime }\left(x^{\left(k\right)}\right)\right]}^{-1}F\left(x^{\left(k\right)}\right),k=0,1,\dots ,\hfill \\ \hfill z^{\left(k\right)}& =y^{\left(k\right)}-\left[3I-2{\left[F^{\prime }\left(x^{\left(k\right)}\right)\right]}^{-1}[y^{\left(k\right)},x^{\left(k\right)};F]\right]{\left[F^{\prime }\left(x^{\left(k\right)}\right)\right]}^{-1}F\left(y^{\left(k\right)}\right),\hfill \\ \hfill x^{(k+1)}& =z^{\left(k\right)}-\left[3I-2{\left[F^{\prime }\left(x^{\left(k\right)}\right)\right]}^{-1}[y^{\left(k\right)},x^{\left(k\right)};F]\right]{\left[F^{\prime }\left(x^{\left(k\right)}\right)\right]}^{-1}F\left(z^{\left(k\right)}\right).\hfill \end{array}$$

(6)

Similarly, its Jacobian-free version was constructed in [11] and denoted by ${\mathbf{M}}_{\mathbf{6},\mathbf{3}}$,

$$\begin{array}{cc}\hfill y^{\left(k\right)}& =x^{\left(k\right)}-{\left[u^{\left(k\right)},x^{\left(k\right)};F\right]}^{-1}F\left(x^{\left(k\right)}\right),k=0,1,\dots ,\hfill \\ \hfill z^{\left(k\right)}& =y^{\left(k\right)}-\left[3I-2{\left[u^{\left(k\right)},x^{\left(k\right)};F\right]}^{-1}[y^{\left(k\right)},x^{\left(k\right)};F]\right]{\left[u^{\left(k\right)},x^{\left(k\right)};F\right]}^{-1}F\left(y^{\left(k\right)}\right),\hfill \\ \hfill x^{(k+1)}& =z^{\left(k\right)}-\left[3I-2{\left[u^{\left(k\right)},x^{\left(k\right)};F\right]}^{-1}[y^{\left(k\right)},x^{\left(k\right)};F]\right]{\left[u^{\left(k\right)},x^{\left(k\right)};F\right]}^{-1}F\left(z^{\left(k\right)}\right).\hfill \end{array}$$

(7)

where again $\left[u^{\left(k\right)},x^{\left(k\right)},F\right]\approx F^{\prime }\left(x^{\left(k\right)}\right)$ and $u^{\left(k\right)}=x^{\left(k\right)}+\alpha G\left(x^{\left(k\right)}\right)$, $m=2$. It should be noticed that in schemes (3), (5) and (7), we employed a quadratic element-by-element power of $F\left(x^{\left(k\right)}\right)$ in the divided differences. This adjustment was essential for preserving the convergence order of the original method (see [11]). However, in our proposal, the order of convergence of the original schemes is held avoiding the computational cost of this element-by-element power.

Therefore, to avoid the calculation of Jacobian matrices, which can be a bottleneck in terms of computational efficiency especially for large systems, this article presents a two-step fifth-order efficient Jacobian-free iterative method that addresses these challenges by eliminating the need for direct Jacobian computation. Our approach is grounded in the use of divided differences and scalar accelerators recently developed in some very efficient schemes (using Jacobian matrices), [17,18]. This not only reduces the computational costs, but also accelerates the convergence with simpler iterative expressions. The proposed method’s design and theoretical underpinnings are discussed, emphasizing its ability to achieve high-order convergence without the Jacobian calculations typically required.

In the Section 2, we develop a new parametric class of Jacobian-free iterative methods using scalar accelerators and demonstrate its theoretical order of convergence, depending on the values of the parameters involved. Subsequently, in Section 3 we carry out an efficiency analysis in which we compare our proposed method with the Jacobian-free versions of others previously cited in the literature. Finally, Section 4 presents practical results of these iterative methods applied to different nonlinear systems of equations.

2. Construction and Convergence of New Jacobian-Free Iterative Method

In 2023, Singh, Sharma and Kumar [18] proposed a family of iterative methods,

$$\begin{array}{cc}\hfill w^{\left(k\right)}& =x^{\left(k\right)}-F^{\prime }{\left(x^{\left(k\right)}\right)}^{-1}F\left(x^{\left(k\right)}\right),k=0,1,2,\dots \hfill \\ \hfill x^{(k+1)}& =w^{\left(k\right)}-\left(p_1+p_2{\displaystyle \frac{F{\left(w^{\left(k\right)}\right)}^{T}F\left(w^{\left(k\right)}\right)}{F{\left(x^{\left(k\right)}\right)}^{T}F\left(x^{\left(k\right)}\right)}}\right)F^{\prime }{\left(w^{\left(k\right)}\right)}^{-1}F\left(w^{\left(k\right)}\right),\hfill \end{array}$$

(8)

where ${\displaystyle \frac{F{\left(y^{\left(k\right)}\right)}^{T}F\left(y^{\left(k\right)}\right)}{F{\left(x^{\left(k\right)}\right)}^{T}F\left(x^{\left(k\right)}\right)}}$ is a scalar accelerator that can be interpreted as $F{\left(y^{\left(k\right)}\right)}^{T}F\left(y^{\left(k\right)}\right)={\parallel F\left(y^{\left(k\right)}\right)\parallel }^2$, and $F{\left(x^{\left(k\right)}\right)}^{T}F\left(x^{\left(k\right)}\right)={\parallel F\left(x^{\left(k\right)}\right)\parallel }^2$, respectively. The real parameters $p_1$ and $p_2$ make the order of convergence of the method five if $p_1$=$p_2=1$, order four if $p_1=1$ and $p_2$ arbitrary and order two if $p_1\ne 1$ and $p_2$ arbitrary. It is known that in many practical applications, computing the Jacobian matrix can be very resource-intensive and time-consuming, therefore, Jacobian-free methods are often preferred.

Making a modification to the scheme (8) by replacing the Jacobian matrices by specific divided differences, we obtain the following family:

$$\begin{array}{cc}\hfill y^{\left(k\right)}& =x^{\left(k\right)}-{\left[u_x^{\left(k\right)},x^{\left(k\right)},F\right]}^{-1}F\left(x^{\left(k\right)}\right),k=0,1,2,\dots \hfill \\ \hfill x^{(k+1)}& =y^{\left(k\right)}-\left(p_1+p_2{\displaystyle \frac{F{\left(y^{\left(k\right)}\right)}^{T}F\left(y^{\left(k\right)}\right)}{F{\left(x^{\left(k\right)}\right)}^{T}F\left(x^{\left(k\right)}\right)}}\right){\left[u_y^{\left(k\right)},y^{\left(k\right)},F\right]}^{-1}F\left(y^{\left(k\right)}\right),\hfill \end{array}$$

(9)

where $\left[u_x^{\left(k\right)},x^{\left(k\right)};F\right]\approx F^{\prime }\left(x^{\left(k\right)}\right)$, $\left[u_y^{\left(k\right)},y^{\left(k\right)};F\right]\approx F^{\prime }\left(y^{\left(k\right)}\right)$, $u_x^{\left(k\right)}=x^{\left(k\right)}+\alpha F\left(x^{\left(k\right)}\right)$ and $u_y^{\left(k\right)}=y^{\left(k\right)}+\alpha F\left(y^{\left(k\right)}\right)$. From now on, we will refer to our modified scheme as $MS(p_1,p_2)$.

The following result shows the error equations arising from method (9) for the possible parameter values, thereby demonstrating that the convergence results of the family (8) hold.

Theorem 1.

Let F be a differentiable enough function $F:\Omega \subseteq {\mathbb{R}}^n\to {\mathbb{R}}^n$ defined in the open convex neighbourhood Ω of ξ, solution of $F\left(x\right)=0$. Let us also consider an initial seed $x^{\left(0\right)}$ near enough to ξ and let $F^{\prime }\left(x\right)$ be continuous and invertible at ξ. Then, the parametric class of iterative schemes presented in (9) locally converges for all $\alpha \in \mathbb{R}$, with the order of convergence given by:

(a) 

Fifth-order convergence if $p_1=p_2=1$, being the corresponding error equation

$$\begin{array}{cc}\hfill x^{(k+1)}-\xi & =(M_4-A_2(M_2{M}_1+M_1{M}_2)-(C_1+C_2)M_2-((M_1^2+M_1)M_2+M_2{M}_1\hfill \\ \hfill & -P^{-1}{M}_1^{2}Q\left)M_1\right)e^{{\left(k\right)}^5}+O\left(e^{{\left(k\right)}^6}\right).\hfill \end{array}$$

(10)

(b) 

Fourth-order convergence if $p_1=1,p_2\ne 1$, being the corresponding error equation

$$\begin{array}{cc}\hfill x^{(k+1)}-\xi & =(-A_2{M}_1^2+C_1{M}_1-p_2{M}_1^3)e^{{\left(k\right)}^4}\hfill \\ \hfill & +(M_4-p_1(A_2(M_2{M}_1+M_1{M}_2)-(C_1+C_2)M_2)-\left(p_2\right((M_1^2+M_1)M_2\hfill \\ \hfill & +M_2{M}_1-P^{-1}{M}_1^{2}Q\left)M_1\right)e^{{\left(k\right)}^5}+O\left(e^{{\left(k\right)}^6}\right).\hfill \end{array}$$

(11)

(c) 

Second-order convergence if $p_1\ne 1,p_{2}arbitrary$, being the corresponding error equation

$$\begin{array}{cc}\hfill x^{(k+1)}-\xi & =M_1(1-p_1)e^{{\left(k\right)}^2}+M_2(1-p_1)e^{{\left(k\right)}^3}\hfill \\ \hfill & +(M_3(1-p_1)-p_1{A}_2{M}_1^2+p_1{C}_1{M}_1-p_2{M}_1^3)e^{{\left(k\right)}^4}\hfill \\ \hfill & +(M_4-p_1(A_2(M_2{M}_1+M_1{M}_2)-(C_1+C_2)M_2)-\left(p_2\right((M_1^2+M_1)M_2\hfill \\ \hfill & +M_2{M}_1-P^{-1}{M}_1^{2}Q\left)M_1\right)e^{{\left(k\right)}^5}+O\left(e^{{\left(k\right)}^6}\right),\hfill \end{array}$$

(12)

being $A_j={\displaystyle \frac{1}{j!}}{\left[F^{\prime }\left(\xi \right)\right]}^{-1}{F}^{\left(j\right)}\left(\xi \right)$, $j=2,3,\dots$, and $C_i,M_i$, $i=1,2,\dots$, are combinations of $A_j$, and denoting the error at iteration k by $e^{\left(k\right)}=x^{\left(k\right)}-\xi$.

Proof. 

Let $e^{\left(k\right)}=x^{\left(k\right)}-\xi$ be the error at the k-th iteration and let $\xi \in {\mathbb{R}}^n$ be a solution of $F\left(x\right)=0$. Then, expanding $F\left(x^{\left(k\right)}\right)$ in the neighborhood of $\xi$, we have

$$\begin{array}{cc}\hfill F\left(x^{\left(k\right)}\right)& =F^{\prime }\left(\xi \right)\left[e^{\left(k\right)}+A_2{e}^{{\left(k\right)}^2}+A_3{e}^{{\left(k\right)}^3}+A_4{e}^{{\left(k\right)}^4}+A_5{e}^{{\left(k\right)}^5}+A_6{e}^{{\left(k\right)}^6}+O\left(e^{{\left(k\right)}^7}\right)\right],\hfill \\ \hfill F^{\prime }\left(x^{\left(k\right)}\right)& =F^{\prime }\left(\xi \right)\left[I+2A_2{e}^{\left(k\right)}+3A_3{e}^{{\left(k\right)}^2}+4A_4{e}^{{\left(k\right)}^3}+5A_5{e}^{{\left(k\right)}^4}+6A_6{e}^{{\left(k\right)}^5}+O\left(e^{{\left(k\right)}^6}\right)\right],\hfill \\ \hfill F^{″}\left(x^{\left(k\right)}\right)& =F^{\prime }\left(\xi \right)\left[2A_2+6A_3{e}^{\left(k\right)}+12A_4{e}^{{\left(k\right)}^2}+20A_5{e}^{{\left(k\right)}^3}+30A_6{e}^{{\left(k\right)}^4}+O\left(e^{{\left(k\right)}^5}\right)\right],\hfill \\ \hfill F^{‴}\left(x^{\left(k\right)}\right)& =F^{\prime }\left(\xi \right)\left[6A_3+24A_4{e}^{\left(k\right)}+60A_5{e}^{{\left(k\right)}^2}+120A_6{e}^{{\left(k\right)}^3}+O\left(e^{{\left(k\right)}^4}\right)\right],\hfill \\ \hfill F^{\left(iv\right)}\left(x^{\left(k\right)}\right)& =F^{\prime }\left(\xi \right)\left[24A_4+120A_5{e}^{\left(k\right)}+360A_6{e}^{{\left(k\right)}^2}+O\left(e^{{\left(k\right)}^3}\right)\right],\hfill \\ \hfill F^{\left(v\right)}\left(x^{\left(k\right)}\right)& =F^{\prime }\left(\xi \right)\left[120A_5+720A_6{e}^{\left(k\right)}+O\left(e^{{\left(k\right)}^2}\right)\right],\hfill \\ \hfill F^{\left(vi\right)}\left(x^{\left(k\right)}\right)& =F^{\prime }\left(\xi \right)\left[720A_6+O\left(e^{\left(k\right)}\right)\right].\hfill \end{array}$$

(13)

Then, based on the formula of Genochi and Hermite (see [3]) we have

$$\begin{array}{cc}\hfill \left[u_x^{\left(k\right)},x^{\left(k\right)};F\right]& =F^{\prime }\left(x^{\left(k\right)}\right)+{\displaystyle \frac{1}{2!}}{F}^{″}\left(x^{\left(k\right)}\right)(u_x^{\left(k\right)}-x^{\left(k\right)})+{\displaystyle \frac{1}{3!}}{F}^{‴}\left(x^{\left(k\right)}\right){(u_x^{\left(k\right)}-x^{\left(k\right)})}^2+{\displaystyle \frac{1}{4!}}{F}^{\left(iv\right)}\left(x^{\left(k\right)}\right){(u_x^{\left(k\right)}-x^{\left(k\right)})}^3\hfill \\ \hfill & +{\displaystyle \frac{1}{5!}}{F}^{\left(v\right)}\left(x^{\left(k\right)}\right){(u_x^{\left(k\right)}-x^{\left(k\right)})}^4+{\displaystyle \frac{1}{6!}}{F}^{\left(vi\right)}\left(x^{\left(k\right)}\right){(u_x^{\left(k\right)}-x^{\left(k\right)})}^5+O\left({(u_{x^{\left(k\right)}}-x^{\left(k\right)})}^6\right).\hfill \end{array}$$

(14)

Taking into account that $u_x^{\left(k\right)}-x^{\left(k\right)}=\alpha F\left(x^{\left(k\right)}\right)$, and performing a series expansion up to fifth-order, we get

$$\begin{array}{cc}\hfill {\alpha }^2{\left(F\left(x^{\left(k\right)}\right)\right)}^2& ={\alpha }^2{\left(F^{\prime }\left(\xi \right)\right)}^2{e}^{{\left(k\right)}^2}+\alpha [{\left(F^{\prime }\left(\xi \right)\right)}^2{A}_2+F^{\prime }\left(\xi \right)A_2{F}^{\prime }\left(\xi \right)]e^{{\left(k\right)}^3}\hfill \\ \hfill & +{\alpha }^2[{\left(F^{\prime }\left(\xi \right)\right)}^2{A}_3+F^{\prime }\left(\xi \right)A_2{F}^{\prime }\left(\xi \right)A_2+F^{\prime }\left(\xi \right)A_3{F}^{\prime }\left(\xi \right)]e^{{\left(k\right)}^4}\hfill \\ \hfill & +{\alpha }^2[{\left(F^{\prime }\left(\xi \right)\right)}^2{A}_4+F^{\prime }\left(\xi \right)A_2{F}^{\prime }\left(\xi \right)A_3+F^{\prime }\left(\xi \right)A_3{F}^{\prime }\left(\xi \right)A_2+F^{\prime }\left(\xi \right)A_4{F}^{\prime }\left(\xi \right)]e^{{\left(k\right)}^5}+O\left(e^{{\left(k\right)}^6}\right),\hfill \\ \\ \hfill {\alpha }^3{\left(F\left(x^{\left(k\right)}\right)\right)}^3& ={\alpha }^3{\left(F^{\prime }\left(\xi \right)\right)}^3{e}^{{\left(k\right)}^3}+{\alpha }^3[{\left(F^{\prime }\left(\xi \right)\right)}^3{A}_2+{\left(F^{\prime }\left(\xi \right)\right)}^2{A}_2{F}^{\prime }\left(\xi \right)+F^{\prime }\left(\xi \right)A_2{\left(F^{\prime }\left(\xi \right)\right)}^2]e^{{\left(k\right)}^4}\hfill \\ \hfill & +{\alpha }^3[{\left(F^{\prime }\left(\xi \right)\right)}^3{A}_3+{\left(F^{\prime }\left(\xi \right)\right)}^2{A}_2{F}^{\prime }\left(\xi \right)A_2+F^{\prime }\left(\xi \right)A_2{\left(F^{\prime }\left(\xi \right)\right)}^2{A}_2+{\left(F^{\prime }\left(\xi \right)\right)}^2{A}_3{F}^{\prime }\left(\xi \right)\hfill \\ \hfill & +F^{\prime }\left(\xi \right)A_2{F}^{\prime }\left(\xi \right)A_2{F}^{\prime }\left(\xi \right)+F^{\prime }\left(\xi \right)A_3{\left(F^{\prime }\left(\xi \right)\right)}^2]e^{{\left(k\right)}^5}+O\left(e^{{\left(k\right)}^6}\right),\hfill \\ \\ \hfill {\alpha }^4{\left(F\left(x^{\left(k\right)}\right)\right)}^4& ={\alpha }^4{\left(F^{\prime }\left(\xi \right)\right)}^4{e}^{{\left(k\right)}^4}+{\alpha }^4[{\left(F^{\prime }\left(\xi \right)\right)}^4{A}_2+{\left(F^{\prime }\left(\xi \right)\right)}^3{A}_2{F}^{\prime }\left(\xi \right)\hfill \\ \hfill & +{\left(F^{\prime }\left(\xi \right)\right)}^2{A}_2{\left(F^{\prime }\left(\xi \right)\right)}^2+F^{\prime }\left(\xi \right)A_2\left({\left(F^{\prime }\left(\xi \right)\right)}^3\right]e^{{\left(k\right)}^5}+O\left(e^{k^6}\right),\hfill \\ \\ \hfill {\alpha }^5{\left(F\left(x^{\left(k\right)}\right)\right)}^5& =\alpha {\left(F^{\prime }\left(\xi \right)\right)}^5{e}^{{\left(k\right)}^5}+O\left(e^{{\left(k\right)}^6}\right).\hfill \end{array}$$

(15)

By combining formulas (13), (15) in the Taylor series expansion (14) we obtain:

$$\begin{array}{cc}\hfill \left[u_x^{\left(k\right)},x^{\left(k\right)};F\right]& =F^{\prime }\left(\xi \right)\left[I+B_1{e}^k+B_2{e}^{{\left(k\right)}^2}+B_3{e}^{{\left(k\right)}^3}+B_4{e}^{{\left(k\right)}^4}+B_5{e}^{{\left(k\right)}^5}+O\left(e^{{\left(k\right)}^6}\right)\right],\hfill \end{array}$$

(16)

where

$$\begin{array}{cc}\hfill B_1& =2A_2+\alpha A_2{F}^{\prime }\left(\xi \right),\hfill \\ \hfill B_2& =3A_3+\alpha A_2{F}^{\prime }\left(\xi \right)A_2+3\alpha A_3{F}^{\prime }\left(\xi \right)+{\alpha }^2{A}_3{\left(F^{\prime }\left(\xi \right)\right)}^2,\hfill \\ \hfill B_3& =4A_4+\alpha A_2{F}^{\prime }\left(\xi \right)A_3+3\alpha A_3{F}^{\prime }\left(\xi \right)A_2+6\alpha A_4{F}^{\prime }\left(\xi \right)+{\alpha }^2{A}_3{\left(F^{\prime }\left(\xi \right)\right)}^2{A}_2+{\alpha }^2{A}_3{F}^{\prime }\left(\xi \right)A_2{F}^{\prime }\left(\xi \right)\hfill \\ \hfill & +{\alpha }^{2}4A_4{\left(F^{\prime }\left(\xi \right)\right)}^2+{\alpha }^3{A}_4{\left(F^{\prime }\left(\xi \right)\right)}^3,\hfill \\ \hfill B_4& =5A_5+\alpha A_2{F}^{\prime }\left(\xi \right)A_4+3\alpha A_3{F}^{\prime }\left(\xi \right)A_3+6\alpha A_4{F}^{\prime }\left(\xi \right)A_2+10\alpha A_5{F}^{\prime }\left(\xi \right)+{\alpha }^2{A}_3{\left(F^{\prime }\left(\xi \right)\right)}^2\hfill \\ \hfill & +{\alpha }^2{A}_3{F}^{\prime }\left(\xi \right)A_2{F}^{\prime }\left(\xi \right)A_2+{\alpha }^2{A}_3{F}^{\prime }\left(\xi \right)A_3{F}^{\prime }\left(\xi \right)+4{\alpha }^2{A}_4{\left(F^{\prime }\left(\xi \right)\right)}^2{A}_2+4{\alpha }^2{A}_4{F}^{\prime }\left(\xi \right)A_2{F}^{\prime }\left(\xi \right)\hfill \\ \hfill & +10{\alpha }^2{A}_5{\left(F^{\prime }\left(\xi \right)\right)}^2+{\alpha }^3{A}_4{\left(F^{\prime }\left(\xi \right)\right)}^3{A}_2+{\alpha }^3{A}_4{\left(F^{\prime }\left(\xi \right)\right)}^2{A}_2{F}^{\prime }\left(\xi \right)+{\alpha }^3{A}_4{F}^{\prime }\left(\xi \right)A_2{\left(F^{\prime }\left(\xi \right)\right)}^2\hfill \\ \hfill & +5{\alpha }^3{A}_5{\left(F^{\prime }\left(\xi \right)\right)}^3+{\alpha }^4{A}_5{\left(F^{\prime }\left(\xi \right)\right)}^4,\hfill \\ \hfill and\\ \\ \hfill B_5& =6A_6+\alpha A_2{F}^{\prime }\left(\xi \right)A_5+3\alpha A_3{F}^{\prime }\left(\xi \right)A_4+6\alpha A_4{F}^{\prime }\left(\xi \right)A_3+10\alpha A_5{F}^{\prime }\left(\xi \right)A_2+15\alpha A_6{F}^{\prime }\left(\xi \right)\hfill \\ \hfill & +{\alpha }^2{A}_3{\left(F^{\prime }\left(\xi \right)\right)}^2{A}_4+{\alpha }^2{A}_3{F}^{\prime }\left(\xi \right)A_2{F}^{\prime }\left(\xi \right)A_3+{\alpha }^2{A}_3{F}^{\prime }\left(\xi \right)A_3{F}^{\prime }\left(\xi \right)A_2+{\alpha }^2{A}_3{F}^{\prime }\left(\xi \right)A_4{F}^{\prime }\left(\xi \right)\hfill \\ \hfill & +4{\alpha }^2{A}_4{\left(F^{\prime }\left(\xi \right)\right)}^2{A}_3+4{\alpha }^2{A}_4{F}^{\prime }\left(\xi \right)A_2{F}^{\prime }\left(\xi \right)A_2+4{\alpha }^2{A}_4{F}^{\prime }\left(\xi \right)A_3{F}^{\prime }\left(\xi \right)+10{\alpha }^2{A}_5{\left(F^{\prime }\left(\xi \right)\right)}^2{A}_2\hfill \\ \hfill & +10{\alpha }^2{A}_5{F}^{\prime }\left(\xi \right)A_2{F}^{\prime }\left(\xi \right)+20{\alpha }^2{A}_6{\left(F^{\prime }\left(\xi \right)\right)}^2+{\alpha }^3{A}_4{\left(F^{\prime }\left(\xi \right)\right)}^3{A}_3+{\alpha }^3{A}_4{\left(F^{\prime }\left(\xi \right)\right)}^2{A}_2{F}^{\prime }\left(\xi \right))A_2\hfill \\ \hfill & +{\alpha }^3{A}_4{F}^{\prime }\left(\xi \right)A_2{\left(F^{\prime }\left(\xi \right)\right)}^2{A}_2+{\alpha }^3{A}_4{\left(F^{\prime }\left(\xi \right)\right)}^2{A}_3{F}^{\prime }\left(\xi \right)+{\alpha }^3{A}_4{F}^{\prime }\left(\xi \right)A_2{F}^{\prime }\left(\xi \right)A_2{F}^{\prime }\left(\xi \right)\hfill \\ \hfill & +{\alpha }^3{A}_4{F}^{\prime }\left(\xi \right)A_3{\left(F^{\prime }\left(\xi \right)\right)}^2+5{\alpha }^3{A}_5{\left(F^{\prime }\left(\xi \right)\right)}^3{A}_2+5{\alpha }^3{A}_5{\left(F^{\prime }\left(\xi \right)\right)}^2{A}_2{F}^{\prime }\left(\xi \right)+5{\alpha }^3{A}_5{F}^{\prime }\left(\xi \right)A_2{\left(F^{\prime }\left(\xi \right)\right)}^2\hfill \\ \hfill & +15{\alpha }^3{A}_6{\left(F^{\prime }\left(\xi \right)\right)}^3+{\alpha }^4{A}_5{\left(F^{\prime }\left(\xi \right)\right)}^4{A}_2+{\alpha }^4{A}_5{\left(F^{\prime }\left(\xi \right)\right)}^3{A}_2{F}^{\prime }\left(\xi \right)+{\alpha }^4{A}_5{\left(F^{\prime }\left(\xi \right)\right)}^2{A}_2{\left(F^{\prime }\left(\xi \right)\right)}^2\hfill \\ \hfill & +{\alpha }^4{A}_5{F}^{\prime }\left(\xi \right)A_2{\left(F^{\prime }\left(\xi \right)\right)}^3+6{\alpha }^4{A}_6{\left(F^{\prime }\left(\xi \right)\right)}^4+{\alpha }^5{A}_6{\left(F^{\prime }\left(\xi \right)\right)}^5.\hfill \end{array}$$

Next, we expand the inverse of the divided difference operator $\left[u_x^{\left(k\right)},x^{\left(k\right)};F\right]$, forcing it to satisfy ${\left[u_x^{\left(k\right)},x^{\left(k\right)};F\right]}^{-1}\left[u_x^{\left(k\right)},x^{\left(k\right)};F\right]=I$,

$$\begin{array}{cc}\hfill {\left[u_x^{\left(k\right)},x^{\left(k\right)};F\right]}^{-1}& =\left[I+X_2{e}^k+X_3{e}^{{\left(k\right)}^2}+X_4{e}^{{\left(k\right)}^3}+X_5{e}^{{\left(k\right)}^4}+O\left(e^{{\left(k\right)}^5}\right)\right]{\left(F^{\prime }\left(\xi \right)\right)}^{-1},\hfill \end{array}$$

(17)

where

$$\begin{array}{cc}\hfill X_2& =-B_1,\hfill \\ \hfill X_3& =B_1^2-B_2,\hfill \\ \hfill X_4& =B_1{B}_2+B_2{B}_1-B_1^3-B_3,\hfill \\ \hfill X_5& =B_1{B}_3+B_3{B}_1+B_1^4-B_4-B_1^2+B_2^2-B_1{B}_2{B}_1-B_2{B}_1^2.\hfill \end{array}$$

(18)

By using the Taylor expansions of $F\left(x^{\left(k\right)}\right)$ defined in (13) and ${\left[u_x^{\left(k\right)},x^{\left(k\right)},F\right]}^{-1}$ obtained in (17), we get the error equation for the first step:

$$\begin{array}{cc}\hfill y^{\left(k\right)}-\xi & =M_1{e}^{{\left(k\right)}^2}+M_2{e}^{{\left(k\right)}^3}+M_3{e}^{{\left(k\right)}^4}+M_4{e}^{{\left(k\right)}^5}+O\left(e^{{\left(k\right)}^6}\right),\hfill \end{array}$$

(19)

where

$$\begin{array}{cc}\hfill M_1& =-(X_2+A_2),\hfill \\ \hfill M_2& =-(A_3+X_2{A}_2+X_3),\hfill \\ \hfill M_3& =-(A_4+X_2{A}_3+X_3{A}_2+X_4),\hfill \\ \hfill M_4& =-(A_5+X_2{A}_4+X_3{A}_3+X_4{A}_2+X_5).\hfill \end{array}$$

(20)

Now, we find the error equation for the second step,

$$\begin{array}{cc}\hfill F\left(y^{\left(k\right)}\right)& =F^{\prime }\left(\xi \right)\left[e_y^{\left(k\right)}+A_2{e}_y^{{\left(k\right)}^2}+O\left(e_y^{{\left(k\right)}^3}\right)\right],\hfill \\ \hfill F^{\prime }\left(y^{\left(k\right)}\right)& =F^{\prime }\left(\xi \right)\left[I+2A_2{e}_y^{\left(k\right)}+3A_3{e}_y^{{\left(k\right)}^2}+O\left(e_y^{{\left(k\right)}^3}\right)\right],\hfill \\ \hfill F^{″}\left(y^{\left(k\right)}\right)& =F^{\prime }\left(\xi \right)\left[2A_2+6A_3{e}_y^{\left(k\right)}+12A_4{e}_y^{{\left(k\right)}^2}+O\left(e_y^{{\left(k\right)}^3}\right)\right],\hfill \\ \hfill F^{‴}\left(y^{\left(k\right)}\right)& =F^{\prime }\left(\xi \right)\left[6A_3+O\left(e_y^{\left(k\right)}\right)\right],\hfill \end{array}$$

(21)

from which arises

$$\begin{array}{cc}\hfill F\left(y^{\left(k\right)}\right)& =F^{\prime }\left(\xi \right)\left[M_1{e}^{{\left(k\right)}^2}+M_2{e}^{{\left(k\right)}^3}+(M_3+A_2{M}_1^2)e^{{\left(k\right)}^4}+A_2(M_2{M}_1+M_1{M}_2)e^{{\left(k\right)}^5}+O\left(e^{{\left(k\right)}^6}\right)\right],\hfill \\ \hfill F^{\prime }\left(y^{\left(k\right)}\right)& =F^{\prime }\left(\xi \right)\left[I+2A_2{M}_1{e}^{{\left(k\right)}^2}+2A_2{M}_2{e}^{{\left(k\right)}^3}+(2A_2{M}_3+3A_3{M}_1^2)e^{{\left(k\right)}^4}\right.\hfill \\ \hfill & \left.+(2A_2{M}_4+3A_3(M_1{M}_2+M_2{M}_1))e^{{\left(k\right)}^5}+O\left(e^{{\left(k\right)}^6}\right)\right],\hfill \\ \hfill F^{″}\left(y^{\left(k\right)}\right)& =F^{\prime }\left(\xi \right)\left[2A_2+6A_3{M}_1{e}^{{\left(k\right)}^2}+6A_3{M}_3{e}^{{\left(k\right)}^3}+O\left(e^{{\left(k\right)}^4}\right)\right],\hfill \\ \hfill F^{‴}\left(y^{\left(k\right)}\right)& =F^{\prime }\left(\xi \right)\left[6A_3+O\left(e^{\left(k\right)}\right)\right].\hfill \end{array}$$

(22)

Then, following the process seen in (14), the expansion of the second difference operator is given by

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \left[u_y^{\left(k\right)},y^{\left(k\right)};F\right]& =F^{\prime }\left(\xi \right)\left[I+(2A_2{M}_1+\alpha A_2{F}^{\prime }\left(\xi \right)M_1)e^{{\left(k\right)}^2}+(2A_2{M}_2+\alpha A_2{F}^{\prime }\left(\xi \right)M_2)e^{{\left(k\right)}^3}\right.\hfill \\ \hfill & +(2A_3{M}_3+3\alpha A_3{M}_1^2+\alpha A_2{F}^{\prime }\left(\xi \right)(M_3+A_2{M}_2^2)+3\alpha A_3{M}_1{F}^{\prime }\left(\xi \right)M_1\hfill \\ \hfill & +\alpha A_3{F}^{\prime }\left(\xi \right)M_1{F}^{\prime }\left(\xi \right)M_1)e^{{\left(k\right)}^4}+\left(2A_2{M}_4+3A_3(M_1{M}_2+M_2{M}_1)\right.\hfill \\ \hfill & +\alpha A_2{F}^{\prime }\left(\xi \right)A_2(M_2{M}_1+M_1{M}_2)+3\alpha A_3{M}_1{F}^{\prime }\left(\xi \right)(M_2+M_1)\hfill \\ \hfill & \left.\left.+\alpha A_3{F}^{\prime }\left(\xi \right)(M_1{F}^{\prime }\left(\xi \right)M_2+M_2{F}^{\prime }\left(\xi \right)M_1)\right)e^{{\left(k\right)}^5}+O\left(e^{{\left(k\right)}^6}\right)\right],\hfill \end{array}\end{array}$$

(23)

that can be expressed as

$$\begin{array}{cc}\hfill \left[u_y^{\left(k\right)},y^{\left(k\right)};F\right]& =F^{\prime }\left(\xi \right)\left[I+C_1{e}^{{\left(k\right)}^2}+C_2{e}^{{\left(k\right)}^3}+C_3{e}^{{\left(k\right)}^4}+C_4{e}^{{\left(k\right)}^5}+O\left(e^{{\left(k\right)}^6}\right)\right],\hfill \end{array}$$

(24)

being

$$\begin{array}{cc}\hfill C_1& =2A_2{M}_1+\alpha A_2{F}^{\prime }\left(\xi \right)M_1,\hfill \\ \hfill C_2& =2A_2{M}_2+\alpha A_2{F}^{\prime }\left(\xi \right)M_2,\hfill \\ \hfill C_3& =2A_3{M}_3+3\alpha A_3{M}_1^2+\alpha A_2{F}^{\prime }\left(\xi \right)(M_3+A_2{M}_2^2)+3\alpha A_3{M}_1{F}^{\prime }\left(\xi \right)M_1+\alpha A_3{F}^{\prime }\left(\xi \right)M_1{F}^{\prime }\left(\xi \right)M_1,\hfill \\ \hfill C_4& =2A_2{M}_4+3A_3(M_1{M}_2+M_2{M}_1)+\alpha A_2{F}^{\prime }\left(\xi \right)A_2(M_2{M}_1+M_1{M}_2)+3\alpha A_3{M}_1{F}^{\prime }\left(\xi \right)(M_2+M_1)\hfill \\ \hfill & +\alpha A_3{F}^{\prime }\left(\xi \right)(M_1{F}^{\prime }\left(\xi \right)M_2+M_2{F}^{\prime }\left(\xi \right)M_1).\hfill \end{array}$$

(25)

Again, we get in a similar way as in (17),

$$\begin{array}{cc}\hfill {\left[u_y^{\left(k\right)},y^{\left(k\right)};F\right]}^{-1}& =-C_1{e}^{{\left(k\right)}^2}-C_2{e}^{{\left(k\right)}^3}+O\left(e^{{\left(k\right)}^4}\right).\hfill \end{array}$$

(26)

Now, we proceed to calculate the expansion of ${\left[u_y^{\left(k\right)},y^{\left(k\right)};F\right]}^{-1}F\left(y^{\left(k\right)}\right)$, obtaining

$$\begin{array}{cc}\hfill {\left[u_y^{\left(k\right)},y^{\left(k\right)};F\right]}^{-1}F\left(y^k\right)=& M_1{e}^{{\left(k\right)}^2}+M_2{e}^{{\left(k\right)}^3}+(M_3-A_2{M}_1^2-C_1{M}_1)e^{{\left(k\right)}^4}\hfill \\ \hfill & +(A_2(M_2{M}_1+M_1{M}_2)-(C_1+C_2)M_2)e^{{\left(k\right)}^5}+O\left(e^{{\left(k\right)}^6}\right).\hfill \end{array}$$

(27)

According to Theorem 1 proven by Singh, Sharma and Kurmar in [18], we have

$${\displaystyle \frac{F{\left(y^{\left(k\right)}\right)}^{T}F\left(y^{\left(k\right)}\right)}{F{\left(x^{\left(k\right)}\right)}^{T}F\left(x^{\left(k\right)}\right)}}={\displaystyle \frac{P{e}_y^{{\left(k\right)}^2}+O\left(e_y^{{\left(k\right)}^3}\right)}{P{e}^{{\left(k\right)}^2}+Q{e}^{{\left(k\right)}^3}+O\left(e^{{\left(k\right)}^4}\right)}},$$

(28)

where

$${\displaystyle P=\sum _{i=1}^n{m}_i{P}_i, \mathrm{with} P_i=R_i^T\times R_i,Q=\sum _{i=1}^n{m}_i{Q}_i, \mathrm{with} Q_i=R_i^T{H}_i+H_i^T{R}_i,}$$

and

$$R_i=\left({\displaystyle \frac{\partial f_i}{\partial x_1}},{\displaystyle \frac{\partial f_i}{\partial x_2}},\dots ,{\displaystyle \frac{\partial f_i}{\partial x_n}}\right),H_i={\displaystyle \frac{1}{2}}{\left({\displaystyle \frac{{\partial }^2{f}_i}{\partial x_j\partial x_r}}\right)}_{n\times n}.$$

Using (19), we get

$$e_y^{{\left(k\right)}^2}=M_1^2{e}^{{\left(k\right)}^4}+(M_1{M}_2+M_2{M}_1)e^{{\left(k\right)}^5}+O\left(e^{{\left(k\right)}^6}\right).$$

(29)

After substituting (29) in (28) when performing the quotient, we obtain:

$${\displaystyle \frac{F{\left(y^{\left(k\right)}\right)}^{T}F\left(y^{\left(k\right)}\right)}{F{\left(x^{\left(k\right)}\right)}^{T}F\left(x^{\left(k\right)}\right)}}=M_1^2{e}^{{\left(k\right)}^2}+\left(M_1{M}_2+M_2{M}_1-P^{-1}{M}_1^{2}Q\right)e^{{\left(k\right)}^3}+O\left(e^{{\left(k\right)}^4}\right).$$

(30)

Finally, fitting (19), (27) and the last result obtained in (30), we have

$$\begin{array}{cc}\hfill x^{(k+1)}-\xi & =M_1(1-p_1)e^{{\left(k\right)}^2}+M_2(1-p_1)e^{{\left(k\right)}^3}\hfill \\ \hfill & +(M_3(1-p_1)-p_1{A}_2{M}_1^2+p_1{C}_1{M}_1-p_2{M}_1^3)e^{{\left(k\right)}^4}\hfill \\ \hfill & +(M_4-p_1(A_2(M_2{M}_1+M_1{M}_2)-(C_1+C_2)M_2)-\left(p_2\right((M_1^2+M_1)M_2\hfill \\ \hfill & +M_2{M}_1-P^{-1}{M}_1^{2}Q\left)M_1\right)e^{{\left(k\right)}^5}+O\left(e^{{\left(k\right)}^6}\right).\hfill \end{array}$$

(31)

In this last result, it is easy to observe that when $p_1$ is different from one ($p_1\ne 1$) the iterative method has order of convergence equal to two, since the term $e^{{\left(k\right)}^2}$ would not cancel out. However, when $p_1=p_2=1$ both terms with $e^{{\left(k\right)}^2}$ and $e^{{\left(k\right)}^3}$ cancel out while the term $e^{{\left(k\right)}^4}$ is as follows

$$-A_2{M}_1^2+C_1{M}_1-M_1^3.$$

Let us remember that $C_1=2A_2{M}_1+\alpha A_2{F}^{\prime }\left(\xi \right)M_1$, so that

$$A_2{M}_1^2+\alpha A_2{F}^{\prime }\left(\xi \right)M_1^2-M_1^3,$$

(32)

but since $M_1=-(X_2+A_2)$, then $X_2=-B_1$ , so $M_1=-(-B_1+A_2)$. Also, $B_1=2A_2+\alpha A_2{F}^{\prime }\left(\xi \right)$. Given that $M_1=A_2+\alpha A_2{F}^{\prime }\left(\xi \right)$, replacing $M_1$ in equation (32), we get

$$(A_2+\alpha A_2{F}^{\prime }\left(\xi \right))M_1^2-M_1^3\to M_1^3-M_1^3=0,$$

resulting in the error equation

$$\begin{array}{cc}\hfill e^{(k+1)}-\xi & =(M_4-A_2(M_2{M}_1+M_1{M}_2)-(C_1+C_2)M_2-((M_1^2+M_1)M_2+M_2{M}_1\hfill \\ \hfill & -P^{-1}{M}_1^{2}Q\left)M_1\right)e^{{\left(k\right)}^5}+O\left(e^{{\left(k\right)}^6}\right).\hfill \end{array}$$

(33)

From the above it is clear to see that if $p_1=1$ but $p_2\ne 1$ only order four is reached. □

3. Efficiency Analysis

We have demonstrated the order of convergence of the proposed class of the iterative method for the different values of the parameters $p_1,p_2$. In this section, we perform a computational effort study considering the effort of solving the involved linear systems per iteration and the other computational cost (functional evaluations, amount of product/quotients,...), not only for the proposed class but also for some Jacobian-free schemes presented in the introductory section.

In order to get this aim, it is known that the needed operations (products/quotients) of solving a $n\times n$ linear system is

$${\displaystyle \frac{1}{3}}{n}^3+n^2-{\displaystyle \frac{1}{3}}n.$$

However, if other linear systems with the same coefficient matrix are solved, then the cost upgrades only in $n^2$ operations each; for each divided difference we calculate $n^2$ quotients; for each functional evaluation of F at different points, a cost of n real evaluations; for each evaluation of a divided differences, $n^2-n$ scalar evaluations; Indeed, a matrix–vector product needs $n^2$ product/quotients. Based on the above, the computational cost for each method appears in Table 1. From family (9), which we call $MS(p_1,p_2)$, we consider the fifth-order member $p_1=p_2=1$ and its fourth-order partner $p_1=1,p_2=-1$. They have the same computational cost, which will be reflected, along with the others, in Table 1.

The results presented in Table 1 show that the method with the highest computational cost is $Traub-Ste$, while those with intermediate costs are $Ostr{o}_{01}$ and $MS(p_1,p_2)$, with the latter being slightly better. The ones that offer the lowest cost are $M_{4,3}$ and $M_{6,3}$, although the latter has sixth-order convergence, which is a factor to consider when obtaining the efficiency index.

In order to show more clearly how computational cost influences the efficiency index $I=p^{{\displaystyle \frac{1}{C}}}$ , where p is the convergence order of the corresponding scheme (see [20]), we present Figure 1 and Figure 2 for different sizes of the nonlinear system to be solved.

In Figure 1, we observe that for systems with dimensions $n=2,3,\dots ,10$, the proposed class of vectorial iterative methods (9) shows a better computational efficiency than the other schemes, for both members of the family. On the other hand, as the system grows to sizes $n=10,11,\dots ,50$ it becomes apparent that methods $M_{4,3}$ and $M_{6,3}$ yield better results. Let us notice that for sizes of $n<11$, our method has better efficiency than $M_{6,3}$, and for sizes $n<12$, it has better efficiency compared to $M_{4,3}$ (see Figure 2).

In the next section, we check the theoretical convergence order of our proposed method and assess the efficiency of different schemes on nonlinear systems of various sizes.

4. Numerical Results

In this section, we test numerically that the theoretical order holds for practical purposes in the proposed Jacobian-free class (9). Moreover, we compare it with the methods appearing in the efficiency section to show their accuracy and computational performance.

Below we show some nonlinear problems, including one whose related nonlinear function is not differentiable, in order to test the applicability of the methods on different kind of problems.

• We consider $F^1\left(x\right)={\left(f_1^1\left(x\right),f_2^1\left(x\right),\dots ,f_{25}^1\left(x\right)\right)}^T$, where

  $$\begin{array}{cc}\hfill f_i^1\left(x\right)& =x_i^2{x}_{i+1}-1,i=1,2,\dots ,24,\hfill \\ \hfill f_{25}^1\left(x\right)& =x_{25}^2{x}_1-1.\hfill \end{array}$$

  whose solution is $\overline{\xi }=(1,1,1,\dots ,1)$.
• We also have $F^2\left(x\right)={\left(f_1^2\left(x\right),f_2^2\left(x\right),\dots ,f_8^2\left(x\right)\right)}^T$, where

  $$f_i^2\left(x\right)=x_i-cos\left(2x_i-\sum _{k=1}^8{x}_k\right),i=1,2,\dots ,8,$$

  being $\overline{\xi }\approx (0.5149,0.5149,\dots ,0.5149)$ its solution.
• We also consider $F^3\left(x\right)={\left(f_1^3\left(x\right),f_2^3\left(x\right),\dots ,f_5^3\left(x\right)\right)}^T$, where

  $$f_i^3\left(x\right)=\sum _{k=1}^5\left(x_k\right)-x_i-e^{(-x_i)},i=1,2,3,4,5,$$

  where the solution is $\overline{\xi }\approx (0.20389,0.20389,0.20389,\dots ,0.20389)$.
• We test also $F^4\left(x\right)={\left(f_1^4\left(x\right),f_2^4\left(x\right),\dots ,f_5^4\left(x\right)\right)}^T$, where

  $$f_i^4\left(x\right)=\sum _{k=1}^5\left(x_k\right)-x_i-e^{(-x_i)}{x}_i,i=1,2,3,4,5,$$

  whose solution is $\overline{\xi }=(0,0,0,\dots ,0)$.
• $F^5\left(x\right)={\left(f_1^5\left(x\right),f_2^5\left(x\right),\dots ,f_{10}^5\left(x\right)\right)}^T$, where

  $$f_i^5\left(x\right)=x_i+1-2log\left(1-x_i+\sum _{k=1}^{10}\left(x_k\right)\right),i=1,2,3,\dots ,10,$$

  being it solution $\overline{\xi }\approx (7.4370,7.4370,\dots ,7.4370)$.
• $F^6\left(x\right)={\left(f_1^6\left(x\right),f_2^6\left(x\right),\dots ,f_5^6\left(x\right)\right)}^T$, where

  $$f_i^6\left(x\right)=x_i+1.5sin\left(\sum _{k=1}^5\left(x_k\right)-x_i\right),i=1,2,3,4,5,$$

  and there exist two solutions, ${\overline{\xi }}_1\approx (-0.3004,-0.3004,\dots ,-0.3004)$, and ${\overline{\xi }}_2\approx (0.4579,0.4579,\dots ,0.4579)$.
• $F^7\left(x\right)={\left(f_1^7\left(x\right),f_2^7\left(x\right)\right)}^T$, where

  $$F_i^7\left(x\right)=arctan\left(x_i\right)+1-2\left[\left({\displaystyle \sum _{k=1}^2{x}_k^2}\right)-x_i^2\right]=0,\mathrm{i}=1,2,$$

  being $\overline{\xi }\approx (0.936,0.936)$.
• We also test $F^8\left(x\right)={\left(f_1^8\left(x\right),f_2^8\left(x\right)\right)}^T$, where

  $$\begin{array}{cc}\hfill f_1^8\left(x\right)& =log\left(\right|x_1\left|\right)+|x_2|,\hfill \\ \hfill f_2^8\left(x\right)& =e^{\left(x_1\right)}+x_2-1,\hfill \end{array}$$

  with solutions ${\overline{\xi }}_1\approx (-0.6275,0.4661)$ and ${\overline{\xi }}_2\approx (0.5122,-0.669)$.

Numerical results have been obtained with Matlab2022b version, using 8000 digits in variable precision arithmetics, a processor AMD $A12-9720P$ RADEON R7, 12 COMPUTE CORES 4C, +8G-Ram, 2.70 GHz. These results are shown in Table 2, Table 3, Table 4, Table 5, Table 6, Table 7, Table 8, Table 9, Table 10, Table 11, Table 12 and Table 13, including the following information, where the appearing norms are Euclidean:

• k: amount of iterations needed ("-" appears when the scheme does not converge or it needs more iterations than the maximum allowed).
• $\overline{\xi }$: obtained solution.
• Cpu-time: average time in seconds required by the iterative method to reach the solution of the problem when executed ten times.
• $\rho$: approximated computational order of convergence, ACOC, firstly appeaaring in [19]

  $$\rho ={\displaystyle \frac{ln\frac{∥x^{(k+1)}-x^{\left(k\right)}∥}{∥x^{\left(k\right)}-x^{(k-1)}∥}}{ln\frac{∥x^{\left(k\right)}-x^{(k-1)}∥}{∥x^{(k-1)}-x^{(k-2)}∥}}},k=2,3,\dots ,$$
  (if $\rho$ is not stable, then "-" appears in the table).
• ${ϵ}_{aprox}=∥x^{(k+1)}-x^{\left(k\right)}∥$.
• ${ϵ}_f=∥F\left(x^{(k+1)}\right)∥$. If ${ϵ}_f$ or ${ϵ}_{aprox}$ are very far from zero or we get infinity or NAN, then "-" appears in the table.

Regarding the stopping criterium, the iterative process ends when one of the following conditions is fulfilled:

  (i)

$∥F\left(x^{(k+1)}\right)∥<{10}^{-100}$,

 (ii)

$∥x^{(k+1)}-x^{\left(k\right)}∥<{10}^{-100}$,

(iii)

50 iterations are reached.

In Table 2, $Ostr{o}_{01}$ method reached the maximum number of iterations without converging to the solution, while the most notable scheme is $M_{6,3}$ in almost all aspects such as errors, iterations and computational time. On the other hand, although $MS(1,1)$ exhibits a lower computational time than $M_{4,3}$, the latter has an additional iteration, a relatively similar time, and better errors, proving that it might even be superior in terms of efficiency.

In Table 3, we observe that our proposed schemes yield the best computational times, with $MS(1,1)$ method standing out for having the smallest error norm among all schemes. The $M_{4,3}$ method, while showing good overall performance in terms of errors, is more computationally expensive as it takes one iteration longer to reach the solution to the system and slightly more time compared to $MS(p1,p2)$ methods, indicating that the program takes more time to generate each iteration.

The results in Table 4 and Table 5 are very balanced, with the $MS(1,1)$ method showing the best errors while the shortest computational time is achieved by the $M_{6,3}$ method. However, this last one takes more effort to generate an iteration because, despite producing one less iteration than the other iterative methods, their execution times are very similar.

For Table 6, we note that the best computational times are obtained by the methods $MS(p1,p2)$. On the other hand, the best errors are obtained by $Ostr{o}_{01}$ scheme, which has a competitive computational time considering it requires four iterations, similar to $Traub-Ste$. The latter has shown to be the one that requires the most time to converge.

In Table 7, the best errors were obtained by $MS(1,-1)$ and $Traub-Ste$, which converged to different solutions, while in terms of performance, $MS(1,1)$ and $Ostr{o}_{01}$ appeared to be better.

In Table 8, $Ostr{o}_{01}$ method could not converge to the solution because it reached the maximum number of iterations, whereas the $MS(p1,p2)$ and $M_{4,3}$ methods showed better overall behavior.

In the Table 9, $MS(1,1)$ method stands out with the lowest computational time, being more efficient than the others, despite requiring a similar number of iterations and showing a convergence order of $4.20$, close to the other methods. The $Traub-Ste$ method is less efficient, with a computational time of $3.1954$ but with a larger errors. $Ostr{o}_{01}$, is competitive only falling slightly behind in computational time. The methods $M_{4,3}$ and $M_{6,3}$ are the most expensive in terms of computational time, requiring 7 and 6 iterations respectively, with the highest time recorded for $M_{4,3}$. Despite their higher orders of convergence, both methods show lower overall efficiency compared to the $MS$ and $Ostro01$ methods.

In Table 10, method $MS(1,1)$ stands out as the one with the lowest computational time, requiring fewer iterations compared to the other methods and showing a convergence order of 4.05. The $Traub-Ste$ method, despite having the same convergence order as $MS(1,-1)$, shows a slightly higher computational time and larger errors in the difference between iterates. $Ostr{o}_{01}$, which requires more iterations, falls only slightly behind in computational time but remains competitive overall. Methods $M_{4,3}$ and $M_{6,3}$ do not converge to the solution, as division by zero occurred when calculating divided differences.

Method $MS(1,1)$ stands out as the one with the lowest computational time, requiring fewer iterations compared to the other methods and showing a convergence order of $4.28$. $Traub-Ste$ scheme, although it has the same convergence order as $MS(1,-1)$, shows a significantly higher computational time but better results in terms of errors. $MS(1,-1)$, while requiring more iterations than $MS(1,1)$, still maintains a competitive time compared to $Traub-Ste$. Division by zero is the reason why the other methods could not reach the solution.

In Table 12, it is observed that only the members of our proposed class converge, and they do it to different solutions. $MS(1,1)$ stands out as the most efficient in terms of time and iterations, while $MS(1,-1)$ yields the best errors. The method $Ostr{o}_{01}$ do not converge because it exceeded the maximum number of iterations, while the other methods could not reach the solution due to division by zero during the computation of divided differences.

The observations for Table 13 are the same as those in the previous one, but now method $Ostr{o}_{01}$ does not converge to the solution due to a division by zero, just like the other methods that do not converge.

5. Conclusions

The fifth-order Jacobian-Free iterative method with scalar accelerators developed in this study (along with its fourth-order partners) has proven to be an efficient tool for solving nonlinear systems of equations. It preserves the convergence order of the original scheme, despite the elimination of Jacobian matrix calculations, with low computational cost. The substitution of these matrices with divided differences not only reduces computational complexity but also facilitates implementation, maintaining high precision and efficiency. The numerical results highlight its superiority in terms of performance compared to conventional methods, especially for not-so-large systems.