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

Preprint

Article

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

Altmetrics

Downloads

Views

Comments

A peer-reviewed article of this preprint also exists.

Alicia Cordero^*

Javier G. Maimó,Antmel Rodríguez-Cabral,

And Juan R. Torregrosa

Alicia Cordero^*

Javier G. Maimó,Antmel Rodríguez-Cabral,

And Juan R. Torregrosa

This version is not peer-reviewed

Submitted:

10 September 2024

Posted:

11 September 2024

You are already at the latest version

Alerts

Abstract

This article introduces a novel two-step fifth-order Jacobian-free iterative method aimed at efficiently solving systems of nonlinear equations. The method leverages the benefits of Jacobian-free approaches, utilizing divided differences to circumvent the computationally intensive calculation of Jacobian matrices. This adaptation significantly reduces computational overhead and simplifies the implementation process while maintaining high convergence rates. We demonstrate that this method achieves fifth-order convergence under specific parameter settings, with a broad applicability across various types of nonlinear systems. The effectiveness of the proposed method is validated through a series of numerical experiments which confirm its superior performance in terms of accuracy and computational efficiency compared to existing methods.

Keywords:

Subject: Computer Science and Mathematics - Computational Mathematics

1. Introduction

Let

F (x) = 0

be a nonlinear system of equations,

F : D \subset R^{n} \to R^{n}

, and the functions

f_{i}

for

i = 1, 2, \dots, n

, are the coordinate components of F, expressed as

F (x) = {(f_{1} (x), f_{2} (x), \dots, f_{n} (x))}^{T}

. Solving nonlinear systems is generally challenging, and solutions

ξ

are typically found by linearizing the problem or employing a fixed-point iteration function

G : D \subset R^{n} \to 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^{(k)} - {[F^{'} (x^{(k)})]}^{- 1} F (x^{(k)}), k = 0, 1, \dots,

being

F^{'} (x^{(k)})

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^{'}

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^{(k)} - {[x^{(k)}, z^{(k)}; F]}^{- 1} F (x^{(k)}), k = 0, 1, \dots,

being

z^{(k)} = x^{(k)} + F (x^{(k)})

, and

[\cdot, \cdot; F] : Ω \times Ω \subset R^{n} \times R^{n} \to L (R^{n})

is the operator of divided differences related to F [3],

[y, x; F] (y - x) = F (y) - F (x) for any x, y \in Ω .

This substitution retains the 2-nd order of convergence while bypassing the calculation of the

F^{'}

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 (x) = 0

, derivative-free iterative methods become more stable when selecting

z = x + α f (x)

for small real values of

α

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{matrix} y^{(k)} & = x^{(k)} - {[F^{'} (x^{(k)})]}^{- 1} F (x^{(k)}), \\ x^{(k + 1)} & = y^{(k)} - {[2 [y^{(k)}, x^{(k)}; F] - F^{'} (x^{(k)})]}^{- 1} F (y^{(k)}), k = 0, 1, \dots, \end{matrix}

achieves only cubic convergence if

F^{'} (x)

is replaced by

[y, x; F]

, as follows:

\begin{matrix} y^{(k)} & = x^{(k)} - {[x^{(k)} + F (x^{(k)}), x^{(k)}; F]}^{- 1} F (x^{(k)}), k = 0, 1, \dots, \\ x^{(k + 1)} & = y^{(k)} - {[2 [y^{(k)}, x^{(k)}; F] - [x^{(k)} + F (x^{(k)}), x^{(k)}; F]]}^{- 1} F (y^{(k)}) . \end{matrix}

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 (x) = {(f_{1} {(x)}^{m}, f_{2} {(x)}^{m}, \dots, f_{n} {(x)}^{m})}^{T}

m \in 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{matrix} y^{(k)} & = x^{(k)} - {[u^{(k)}, x^{(k)}; F]}^{- 1} F (x^{(k)}), k = 0, 1, 2, \dots \\ x^{(k + 1)} & = y^{(k)} - {[y^{(k)}, x^{(k)}; F]}^{- 1} [u^{(k)}, x^{(k)}; F] {[u^{(k)}, y^{(k)}; F]}^{- 1} F (y^{(k)}), \end{matrix}

(1)

where

u^{(k)} = x^{(k)} + β F (x^{(k)}), β \in R

. The iterative schemes have a fourth order of convergence for every

β

β \neq 0

, for our purposes we choose

β = 1

and this be called

Traub - 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{matrix} y^{(k)} & = x^{(k)} - {[F^{'} (x^{(k)})]}^{- 1} F (x^{(k)}), \\ x^{(k + 1)} & = y^{(k)} - {[2 [y^{(k)}, x^{(k)}; F] - F^{'} (x^{(k)})]}^{- 1} F (y^{(k)}), \end{matrix}

(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{matrix} y^{(k)} & = x^{(k)} - {[u^{(k)}, x^{(k)}; F]}^{- 1} F (x^{(k)}), \\ x^{(k + 1)} & = y^{(k)} - {[2 [y^{(k)}, x^{(k)}; F] - [u^{(k)}, x^{(k)}; F]]}^{- 1} F (y^{(k)}), \end{matrix}

(3)

where

[u^{(k)}, x^{(k)}, F] \approx F^{'} (x^{(k)})

, and

u^{(k)} = x^{(k)} + α G (x^{(k)})

(with

m = 2

α \in R

. We denote this method as

{Ostro}_{01}

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

\begin{matrix} y^{(k)} & = x^{(k)} - {[F^{'} (x^{(k)})]}^{- 1} F (x^{(k)}), \\ x^{(k + 1)} & = y^{(k)} - [3 I - 2 {[F^{'} (x^{(k)})]}^{- 1} [y^{(k)}, x^{(k)}; F]] {[F^{'} (x^{(k)})]}^{- 1} F (y^{(k)}), \end{matrix}

(4)

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

\begin{matrix} y^{(k)} & = x^{(k)} - {[u^{(k)}, x^{(k)}; F]}^{- 1} F (x^{(k)}), \\ x^{(k + 1)} & = y^{(k)} - [3 I - 2 {[u^{(k)}, x^{(k)}; F]}^{- 1} [y^{(k)}, x^{(k)}; F]] {[u^{(k)}, x^{(k)}; F]}^{- 1} F (y^{(k)}), \end{matrix}

(5)

that we denote by

M_{4, 3}

, where

[u^{(k)}, x^{(k)}, F] \approx F^{'} (x^{(k)})

and

u^{(k)} = x^{(k)} + α G (x^{(k)})

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{matrix} y^{(k)} & = x^{(k)} - {[F^{'} (x^{(k)})]}^{- 1} F (x^{(k)}), k = 0, 1, \dots, \\ z^{(k)} & = y^{(k)} - [3 I - 2 {[F^{'} (x^{(k)})]}^{- 1} [y^{(k)}, x^{(k)}; F]] {[F^{'} (x^{(k)})]}^{- 1} F (y^{(k)}), \\ x^{(k + 1)} & = z^{(k)} - [3 I - 2 {[F^{'} (x^{(k)})]}^{- 1} [y^{(k)}, x^{(k)}; F]] {[F^{'} (x^{(k)})]}^{- 1} F (z^{(k)}) . \end{matrix}

(6)

Similarly, its Jacobian-free version was constructed in [11] and denoted by

M_{6, 3}

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

(7)

where again

[u^{(k)}, x^{(k)}, F] \approx F^{'} (x^{(k)})

and

u^{(k)} = x^{(k)} + α G (x^{(k)})

m = 2

. It should be noticed that in schemes (3), (5) and (7), we employed a quadratic element-by-element power of

F (x^{(k)})

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{matrix} w^{(k)} & = x^{(k)} - F^{'} {(x^{(k)})}^{- 1} F (x^{(k)}), k = 0, 1, 2, \dots \\ x^{(k + 1)} & = w^{(k)} - (p_{1} + p_{2} \frac{F {(w^{(k)})}^{T} F (w^{(k)})}{F {(x^{(k)})}^{T} F (x^{(k)})}) F^{'} {(w^{(k)})}^{- 1} F (w^{(k)}), \end{matrix}

(8)

where

\frac{F {(y^{(k)})}^{T} F (y^{(k)})}{F {(x^{(k)})}^{T} F (x^{(k)})}

is a scalar accelerator that can be interpreted as

F {(y^{(k)})}^{T} F (y^{(k)}) = {∥ F (y^{(k)}) ∥}^{2}

, and

F {(x^{(k)})}^{T} F (x^{(k)}) = {∥ F (x^{(k)}) ∥}^{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} \neq 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{matrix} y^{(k)} & = x^{(k)} - {[u_{x}^{(k)}, x^{(k)}, F]}^{- 1} F (x^{(k)}), k = 0, 1, 2, \dots \\ x^{(k + 1)} & = y^{(k)} - (p_{1} + p_{2} \frac{F {(y^{(k)})}^{T} F (y^{(k)})}{F {(x^{(k)})}^{T} F (x^{(k)})}) {[u_{y}^{(k)}, y^{(k)}, F]}^{- 1} F (y^{(k)}), \end{matrix}

(9)

where

[u_{x}^{(k)}, x^{(k)}; F] \approx F^{'} (x^{(k)})

[u_{y}^{(k)}, y^{(k)}; F] \approx F^{'} (y^{(k)})

u_{x}^{(k)} = x^{(k)} + α F (x^{(k)})

and

u_{y}^{(k)} = y^{(k)} + α F (y^{(k)})

. From now on, we will refer to our modified scheme as

M S (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 : Ω \subseteq R^{n} \to R^{n}

defined in the open convex neighbourhood Ω of ξ, solution of

F (x) = 0

. Let us also consider an initial seed

x^{(0)}

near enough to ξ and let

F^{'} (x)

be continuous and invertible at ξ. Then, the parametric class of iterative schemes presented in (9) locally converges for all

α \in R

, with the order of convergence given by:

(a): Fifth-order convergence if $p_{1} = p_{2} = 1$ , being the corresponding error equation

$\begin{matrix} x^{(k + 1)} - ξ & = (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} \\ - P^{- 1} M_{1}^{2} Q) M_{1}) e^{{(k)}^{5}} + O (e^{{(k)}^{6}}) . \end{matrix}$

(10)
(b): Fourth-order convergence if $p_{1} = 1, p_{2} \neq 1$ , being the corresponding error equation

$\begin{matrix} x^{(k + 1)} - ξ & = (- A_{2} M_{1}^{2} + C_{1} M_{1} - p_{2} M_{1}^{3}) e^{{(k)}^{4}} \\ + (M_{4} - p_{1} (A_{2} (M_{2} M_{1} + M_{1} M_{2}) - (C_{1} + C_{2}) M_{2}) - (p_{2} ((M_{1}^{2} + M_{1}) M_{2} \\ + M_{2} M_{1} - P^{- 1} M_{1}^{2} Q) M_{1}) e^{{(k)}^{5}} + O (e^{{(k)}^{6}}) . \end{matrix}$

(11)
(c): Second-order convergence if $p_{1} \neq 1, p_{2} a r b i t r a r y$ , being the corresponding error equation

$\begin{matrix} x^{(k + 1)} - ξ & = M_{1} (1 - p_{1}) e^{{(k)}^{2}} + M_{2} (1 - p_{1}) e^{{(k)}^{3}} \\ + (M_{3} (1 - p_{1}) - p_{1} A_{2} M_{1}^{2} + p_{1} C_{1} M_{1} - p_{2} M_{1}^{3}) e^{{(k)}^{4}} \\ + (M_{4} - p_{1} (A_{2} (M_{2} M_{1} + M_{1} M_{2}) - (C_{1} + C_{2}) M_{2}) - (p_{2} ((M_{1}^{2} + M_{1}) M_{2} \\ + M_{2} M_{1} - P^{- 1} M_{1}^{2} Q) M_{1}) e^{{(k)}^{5}} + O (e^{{(k)}^{6}}), \end{matrix}$

(12)

being

A_{j} = \frac{1}{j!} {[F^{'} (ξ)]}^{- 1} F^{(j)} (ξ)

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^{(k)} = x^{(k)} - ξ

Proof.

Let

e^{(k)} = x^{(k)} - ξ

be the error at the k-th iteration and let

ξ \in R^{n}

be a solution of

F (x) = 0

. Then, expanding

F (x^{(k)})

in the neighborhood of

ξ

, we have

\begin{matrix} F (x^{(k)}) & = F^{'} (ξ) [e^{(k)} + A_{2} e^{{(k)}^{2}} + A_{3} e^{{(k)}^{3}} + A_{4} e^{{(k)}^{4}} + A_{5} e^{{(k)}^{5}} + A_{6} e^{{(k)}^{6}} + O (e^{{(k)}^{7}})], \\ F^{'} (x^{(k)}) & = F^{'} (ξ) [I + 2 A_{2} e^{(k)} + 3 A_{3} e^{{(k)}^{2}} + 4 A_{4} e^{{(k)}^{3}} + 5 A_{5} e^{{(k)}^{4}} + 6 A_{6} e^{{(k)}^{5}} + O (e^{{(k)}^{6}})], \\ F^{″} (x^{(k)}) & = F^{'} (ξ) [2 A_{2} + 6 A_{3} e^{(k)} + 12 A_{4} e^{{(k)}^{2}} + 20 A_{5} e^{{(k)}^{3}} + 30 A_{6} e^{{(k)}^{4}} + O (e^{{(k)}^{5}})], \\ F^{‴} (x^{(k)}) & = F^{'} (ξ) [6 A_{3} + 24 A_{4} e^{(k)} + 60 A_{5} e^{{(k)}^{2}} + 120 A_{6} e^{{(k)}^{3}} + O (e^{{(k)}^{4}})], \\ F^{(i v)} (x^{(k)}) & = F^{'} (ξ) [24 A_{4} + 120 A_{5} e^{(k)} + 360 A_{6} e^{{(k)}^{2}} + O (e^{{(k)}^{3}})], \\ F^{(v)} (x^{(k)}) & = F^{'} (ξ) [120 A_{5} + 720 A_{6} e^{(k)} + O (e^{{(k)}^{2}})], \\ F^{(v i)} (x^{(k)}) & = F^{'} (ξ) [720 A_{6} + O (e^{(k)})] . \end{matrix}

(13)

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

\begin{matrix} [u_{x}^{(k)}, x^{(k)}; F] & = F^{'} (x^{(k)}) + \frac{1}{2!} F^{″} (x^{(k)}) (u_{x}^{(k)} - x^{(k)}) + \frac{1}{3!} F^{‴} (x^{(k)}) {(u_{x}^{(k)} - x^{(k)})}^{2} + \frac{1}{4!} F^{(i v)} (x^{(k)}) {(u_{x}^{(k)} - x^{(k)})}^{3} \\ + \frac{1}{5!} F^{(v)} (x^{(k)}) {(u_{x}^{(k)} - x^{(k)})}^{4} + \frac{1}{6!} F^{(v i)} (x^{(k)}) {(u_{x}^{(k)} - x^{(k)})}^{5} + O ({(u_{x^{(k)}} - x^{(k)})}^{6}) . \end{matrix}

(14)

Taking into account that

u_{x}^{(k)} - x^{(k)} = α F (x^{(k)})

, and performing a series expansion up to fifth-order, we get

\begin{matrix} α^{2} {(F (x^{(k)}))}^{2} & = α^{2} {(F^{'} (ξ))}^{2} e^{{(k)}^{2}} + α [{(F^{'} (ξ))}^{2} A_{2} + F^{'} (ξ) A_{2} F^{'} (ξ)] e^{{(k)}^{3}} \\ + α^{2} [{(F^{'} (ξ))}^{2} A_{3} + F^{'} (ξ) A_{2} F^{'} (ξ) A_{2} + F^{'} (ξ) A_{3} F^{'} (ξ)] e^{{(k)}^{4}} \\ + α^{2} [{(F^{'} (ξ))}^{2} A_{4} + F^{'} (ξ) A_{2} F^{'} (ξ) A_{3} + F^{'} (ξ) A_{3} F^{'} (ξ) A_{2} + F^{'} (ξ) A_{4} F^{'} (ξ)] e^{{(k)}^{5}} + O (e^{{(k)}^{6}}), \\ α^{3} {(F (x^{(k)}))}^{3} & = α^{3} {(F^{'} (ξ))}^{3} e^{{(k)}^{3}} + α^{3} [{(F^{'} (ξ))}^{3} A_{2} + {(F^{'} (ξ))}^{2} A_{2} F^{'} (ξ) + F^{'} (ξ) A_{2} {(F^{'} (ξ))}^{2}] e^{{(k)}^{4}} \\ + α^{3} [{(F^{'} (ξ))}^{3} A_{3} + {(F^{'} (ξ))}^{2} A_{2} F^{'} (ξ) A_{2} + F^{'} (ξ) A_{2} {(F^{'} (ξ))}^{2} A_{2} + {(F^{'} (ξ))}^{2} A_{3} F^{'} (ξ) \\ + F^{'} (ξ) A_{2} F^{'} (ξ) A_{2} F^{'} (ξ) + F^{'} (ξ) A_{3} {(F^{'} (ξ))}^{2}] e^{{(k)}^{5}} + O (e^{{(k)}^{6}}), \\ α^{4} {(F (x^{(k)}))}^{4} & = α^{4} {(F^{'} (ξ))}^{4} e^{{(k)}^{4}} + α^{4} [{(F^{'} (ξ))}^{4} A_{2} + {(F^{'} (ξ))}^{3} A_{2} F^{'} (ξ) \\ + {(F^{'} (ξ))}^{2} A_{2} {(F^{'} (ξ))}^{2} + F^{'} (ξ) A_{2} ({(F^{'} (ξ))}^{3}] e^{{(k)}^{5}} + O (e^{k^{6}}), \\ α^{5} {(F (x^{(k)}))}^{5} & = α {(F^{'} (ξ))}^{5} e^{{(k)}^{5}} + O (e^{{(k)}^{6}}) . \end{matrix}

(15)

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

\begin{matrix} [u_{x}^{(k)}, x^{(k)}; F] & = F^{'} (ξ) [I + B_{1} e^{k} + B_{2} e^{{(k)}^{2}} + B_{3} e^{{(k)}^{3}} + B_{4} e^{{(k)}^{4}} + B_{5} e^{{(k)}^{5}} + O (e^{{(k)}^{6}})], \end{matrix}

(16)

where

\begin{matrix} B_{1} & = 2 A_{2} + α A_{2} F^{'} (ξ), \\ B_{2} & = 3 A_{3} + α A_{2} F^{'} (ξ) A_{2} + 3 α A_{3} F^{'} (ξ) + α^{2} A_{3} {(F^{'} (ξ))}^{2}, \\ B_{3} & = 4 A_{4} + α A_{2} F^{'} (ξ) A_{3} + 3 α A_{3} F^{'} (ξ) A_{2} + 6 α A_{4} F^{'} (ξ) + α^{2} A_{3} {(F^{'} (ξ))}^{2} A_{2} + α^{2} A_{3} F^{'} (ξ) A_{2} F^{'} (ξ) \\ + α^{2} 4 A_{4} {(F^{'} (ξ))}^{2} + α^{3} A_{4} {(F^{'} (ξ))}^{3}, \\ B_{4} & = 5 A_{5} + α A_{2} F^{'} (ξ) A_{4} + 3 α A_{3} F^{'} (ξ) A_{3} + 6 α A_{4} F^{'} (ξ) A_{2} + 10 α A_{5} F^{'} (ξ) + α^{2} A_{3} {(F^{'} (ξ))}^{2} \\ + α^{2} A_{3} F^{'} (ξ) A_{2} F^{'} (ξ) A_{2} + α^{2} A_{3} F^{'} (ξ) A_{3} F^{'} (ξ) + 4 α^{2} A_{4} {(F^{'} (ξ))}^{2} A_{2} + 4 α^{2} A_{4} F^{'} (ξ) A_{2} F^{'} (ξ) \\ + 10 α^{2} A_{5} {(F^{'} (ξ))}^{2} + α^{3} A_{4} {(F^{'} (ξ))}^{3} A_{2} + α^{3} A_{4} {(F^{'} (ξ))}^{2} A_{2} F^{'} (ξ) + α^{3} A_{4} F^{'} (ξ) A_{2} {(F^{'} (ξ))}^{2} \\ + 5 α^{3} A_{5} {(F^{'} (ξ))}^{3} + α^{4} A_{5} {(F^{'} (ξ))}^{4}, \\ a n d \\ B_{5} & = 6 A_{6} + α A_{2} F^{'} (ξ) A_{5} + 3 α A_{3} F^{'} (ξ) A_{4} + 6 α A_{4} F^{'} (ξ) A_{3} + 10 α A_{5} F^{'} (ξ) A_{2} + 15 α A_{6} F^{'} (ξ) \\ + α^{2} A_{3} {(F^{'} (ξ))}^{2} A_{4} + α^{2} A_{3} F^{'} (ξ) A_{2} F^{'} (ξ) A_{3} + α^{2} A_{3} F^{'} (ξ) A_{3} F^{'} (ξ) A_{2} + α^{2} A_{3} F^{'} (ξ) A_{4} F^{'} (ξ) \\ + 4 α^{2} A_{4} {(F^{'} (ξ))}^{2} A_{3} + 4 α^{2} A_{4} F^{'} (ξ) A_{2} F^{'} (ξ) A_{2} + 4 α^{2} A_{4} F^{'} (ξ) A_{3} F^{'} (ξ) + 10 α^{2} A_{5} {(F^{'} (ξ))}^{2} A_{2} \\ + 10 α^{2} A_{5} F^{'} (ξ) A_{2} F^{'} (ξ) + 20 α^{2} A_{6} {(F^{'} (ξ))}^{2} + α^{3} A_{4} {(F^{'} (ξ))}^{3} A_{3} + α^{3} A_{4} {(F^{'} (ξ))}^{2} A_{2} F^{'} (ξ)) A_{2} \\ + α^{3} A_{4} F^{'} (ξ) A_{2} {(F^{'} (ξ))}^{2} A_{2} + α^{3} A_{4} {(F^{'} (ξ))}^{2} A_{3} F^{'} (ξ) + α^{3} A_{4} F^{'} (ξ) A_{2} F^{'} (ξ) A_{2} F^{'} (ξ) \\ + α^{3} A_{4} F^{'} (ξ) A_{3} {(F^{'} (ξ))}^{2} + 5 α^{3} A_{5} {(F^{'} (ξ))}^{3} A_{2} + 5 α^{3} A_{5} {(F^{'} (ξ))}^{2} A_{2} F^{'} (ξ) + 5 α^{3} A_{5} F^{'} (ξ) A_{2} {(F^{'} (ξ))}^{2} \\ + 15 α^{3} A_{6} {(F^{'} (ξ))}^{3} + α^{4} A_{5} {(F^{'} (ξ))}^{4} A_{2} + α^{4} A_{5} {(F^{'} (ξ))}^{3} A_{2} F^{'} (ξ) + α^{4} A_{5} {(F^{'} (ξ))}^{2} A_{2} {(F^{'} (ξ))}^{2} \\ + α^{4} A_{5} F^{'} (ξ) A_{2} {(F^{'} (ξ))}^{3} + 6 α^{4} A_{6} {(F^{'} (ξ))}^{4} + α^{5} A_{6} {(F^{'} (ξ))}^{5} . \end{matrix}

Next, we expand the inverse of the divided difference operator

[u_{x}^{(k)}, x^{(k)}; F]

, forcing it to satisfy

{[u_{x}^{(k)}, x^{(k)}; F]}^{- 1} [u_{x}^{(k)}, x^{(k)}; F] = I

\begin{matrix} {[u_{x}^{(k)}, x^{(k)}; F]}^{- 1} & = [I + X_{2} e^{k} + X_{3} e^{{(k)}^{2}} + X_{4} e^{{(k)}^{3}} + X_{5} e^{{(k)}^{4}} + O (e^{{(k)}^{5}})] {(F^{'} (ξ))}^{- 1}, \end{matrix}

(17)

where

\begin{matrix} X_{2} & = - B_{1}, \\ X_{3} & = B_{1}^{2} - B_{2}, \\ X_{4} & = B_{1} B_{2} + B_{2} B_{1} - B_{1}^{3} - B_{3}, \\ 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} . \end{matrix}

(18)

By using the Taylor expansions of

F (x^{(k)})

defined in (13) and

{[u_{x}^{(k)}, x^{(k)}, F]}^{- 1}

obtained in (17), we get the error equation for the first step:

\begin{matrix} y^{(k)} - ξ & = M_{1} e^{{(k)}^{2}} + M_{2} e^{{(k)}^{3}} + M_{3} e^{{(k)}^{4}} + M_{4} e^{{(k)}^{5}} + O (e^{{(k)}^{6}}), \end{matrix}

(19)

where

\begin{matrix} M_{1} & = - (X_{2} + A_{2}), \\ M_{2} & = - (A_{3} + X_{2} A_{2} + X_{3}), \\ M_{3} & = - (A_{4} + X_{2} A_{3} + X_{3} A_{2} + X_{4}), \\ M_{4} & = - (A_{5} + X_{2} A_{4} + X_{3} A_{3} + X_{4} A_{2} + X_{5}) . \end{matrix}

(20)

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

\begin{matrix} F (y^{(k)}) & = F^{'} (ξ) [e_{y}^{(k)} + A_{2} e_{y}^{{(k)}^{2}} + O (e_{y}^{{(k)}^{3}})], \\ F^{'} (y^{(k)}) & = F^{'} (ξ) [I + 2 A_{2} e_{y}^{(k)} + 3 A_{3} e_{y}^{{(k)}^{2}} + O (e_{y}^{{(k)}^{3}})], \\ F^{″} (y^{(k)}) & = F^{'} (ξ) [2 A_{2} + 6 A_{3} e_{y}^{(k)} + 12 A_{4} e_{y}^{{(k)}^{2}} + O (e_{y}^{{(k)}^{3}})], \\ F^{‴} (y^{(k)}) & = F^{'} (ξ) [6 A_{3} + O (e_{y}^{(k)})], \end{matrix}

(21)

from which arises

\begin{matrix} F (y^{(k)}) & = F^{'} (ξ) [M_{1} e^{{(k)}^{2}} + M_{2} e^{{(k)}^{3}} + (M_{3} + A_{2} M_{1}^{2}) e^{{(k)}^{4}} + A_{2} (M_{2} M_{1} + M_{1} M_{2}) e^{{(k)}^{5}} + O (e^{{(k)}^{6}})], \\ F^{'} (y^{(k)}) & = F^{'} (ξ) [I + 2 A_{2} M_{1} e^{{(k)}^{2}} + 2 A_{2} M_{2} e^{{(k)}^{3}} + (2 A_{2} M_{3} + 3 A_{3} M_{1}^{2}) e^{{(k)}^{4}} \\ + (2 A_{2} M_{4} + 3 A_{3} (M_{1} M_{2} + M_{2} M_{1})) e^{{(k)}^{5}} + O (e^{{(k)}^{6}})], \\ F^{″} (y^{(k)}) & = F^{'} (ξ) [2 A_{2} + 6 A_{3} M_{1} e^{{(k)}^{2}} + 6 A_{3} M_{3} e^{{(k)}^{3}} + O (e^{{(k)}^{4}})], \\ F^{‴} (y^{(k)}) & = F^{'} (ξ) [6 A_{3} + O (e^{(k)})] . \end{matrix}

(22)

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

\begin{matrix} \begin{matrix} [u_{y}^{(k)}, y^{(k)}; F] & = F^{'} (ξ) [I + (2 A_{2} M_{1} + α A_{2} F^{'} (ξ) M_{1}) e^{{(k)}^{2}} + (2 A_{2} M_{2} + α A_{2} F^{'} (ξ) M_{2}) e^{{(k)}^{3}} \\ + (2 A_{3} M_{3} + 3 α A_{3} M_{1}^{2} + α A_{2} F^{'} (ξ) (M_{3} + A_{2} M_{2}^{2}) + 3 α A_{3} M_{1} F^{'} (ξ) M_{1} \\ + α A_{3} F^{'} (ξ) M_{1} F^{'} (ξ) M_{1}) e^{{(k)}^{4}} + (2 A_{2} M_{4} + 3 A_{3} (M_{1} M_{2} + M_{2} M_{1}) \\ + α A_{2} F^{'} (ξ) A_{2} (M_{2} M_{1} + M_{1} M_{2}) + 3 α A_{3} M_{1} F^{'} (ξ) (M_{2} + M_{1}) \\ + α A_{3} F^{'} (ξ) (M_{1} F^{'} (ξ) M_{2} + M_{2} F^{'} (ξ) M_{1})) e^{{(k)}^{5}} + O (e^{{(k)}^{6}})], \end{matrix} \end{matrix}

(23)

that can be expressed as

\begin{matrix} [u_{y}^{(k)}, y^{(k)}; F] & = F^{'} (ξ) [I + C_{1} e^{{(k)}^{2}} + C_{2} e^{{(k)}^{3}} + C_{3} e^{{(k)}^{4}} + C_{4} e^{{(k)}^{5}} + O (e^{{(k)}^{6}})], \end{matrix}

(24)

being

\begin{matrix} C_{1} & = 2 A_{2} M_{1} + α A_{2} F^{'} (ξ) M_{1}, \\ C_{2} & = 2 A_{2} M_{2} + α A_{2} F^{'} (ξ) M_{2}, \\ C_{3} & = 2 A_{3} M_{3} + 3 α A_{3} M_{1}^{2} + α A_{2} F^{'} (ξ) (M_{3} + A_{2} M_{2}^{2}) + 3 α A_{3} M_{1} F^{'} (ξ) M_{1} + α A_{3} F^{'} (ξ) M_{1} F^{'} (ξ) M_{1}, \\ C_{4} & = 2 A_{2} M_{4} + 3 A_{3} (M_{1} M_{2} + M_{2} M_{1}) + α A_{2} F^{'} (ξ) A_{2} (M_{2} M_{1} + M_{1} M_{2}) + 3 α A_{3} M_{1} F^{'} (ξ) (M_{2} + M_{1}) \\ + α A_{3} F^{'} (ξ) (M_{1} F^{'} (ξ) M_{2} + M_{2} F^{'} (ξ) M_{1}) . \end{matrix}

(25)

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

\begin{matrix} {[u_{y}^{(k)}, y^{(k)}; F]}^{- 1} & = - C_{1} e^{{(k)}^{2}} - C_{2} e^{{(k)}^{3}} + O (e^{{(k)}^{4}}) . \end{matrix}

(26)

Now, we proceed to calculate the expansion of

{[u_{y}^{(k)}, y^{(k)}; F]}^{- 1} F (y^{(k)})

, obtaining

\begin{matrix} {[u_{y}^{(k)}, y^{(k)}; F]}^{- 1} F (y^{k}) = & M_{1} e^{{(k)}^{2}} + M_{2} e^{{(k)}^{3}} + (M_{3} - A_{2} M_{1}^{2} - C_{1} M_{1}) e^{{(k)}^{4}} \\ + (A_{2} (M_{2} M_{1} + M_{1} M_{2}) - (C_{1} + C_{2}) M_{2}) e^{{(k)}^{5}} + O (e^{{(k)}^{6}}) . \end{matrix}

(27)

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

\frac{F {(y^{(k)})}^{T} F (y^{(k)})}{F {(x^{(k)})}^{T} F (x^{(k)})} = \frac{P e_{y}^{{(k)}^{2}} + O (e_{y}^{{(k)}^{3}})}{P e^{{(k)}^{2}} + Q e^{{(k)}^{3}} + O (e^{{(k)}^{4}})},

(28)

where

P = \sum_{i = 1}^{n} m_{i} P_{i}, with P_{i} = R_{i}^{T} \times R_{i}, Q = \sum_{i = 1}^{n} m_{i} Q_{i}, with Q_{i} = R_{i}^{T} H_{i} + H_{i}^{T} R_{i},

and

R_{i} = (\frac{\partial f_{i}}{\partial x_{1}}, \frac{\partial f_{i}}{\partial x_{2}}, \dots, \frac{\partial f_{i}}{\partial x_{n}}), H_{i} = \frac{1}{2} {(\frac{\partial^{2} f_{i}}{\partial x_{j} \partial x_{r}})}_{n \times n} .

Using (19), we get

e_{y}^{{(k)}^{2}} = M_{1}^{2} e^{{(k)}^{4}} + (M_{1} M_{2} + M_{2} M_{1}) e^{{(k)}^{5}} + O (e^{{(k)}^{6}}) .

(29)

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

\frac{F {(y^{(k)})}^{T} F (y^{(k)})}{F {(x^{(k)})}^{T} F (x^{(k)})} = M_{1}^{2} e^{{(k)}^{2}} + (M_{1} M_{2} + M_{2} M_{1} - P^{- 1} M_{1}^{2} Q) e^{{(k)}^{3}} + O (e^{{(k)}^{4}}) .

(30)

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

\begin{matrix} x^{(k + 1)} - ξ & = M_{1} (1 - p_{1}) e^{{(k)}^{2}} + M_{2} (1 - p_{1}) e^{{(k)}^{3}} \\ + (M_{3} (1 - p_{1}) - p_{1} A_{2} M_{1}^{2} + p_{1} C_{1} M_{1} - p_{2} M_{1}^{3}) e^{{(k)}^{4}} \\ + (M_{4} - p_{1} (A_{2} (M_{2} M_{1} + M_{1} M_{2}) - (C_{1} + C_{2}) M_{2}) - (p_{2} ((M_{1}^{2} + M_{1}) M_{2} \\ + M_{2} M_{1} - P^{- 1} M_{1}^{2} Q) M_{1}) e^{{(k)}^{5}} + O (e^{{(k)}^{6}}) . \end{matrix}

(31)

In this last result, it is easy to observe that when

p_{1}

is different from one (

p_{1} \neq 1

) the iterative method has order of convergence equal to two, since the term

e^{{(k)}^{2}}

would not cancel out. However, when

p_{1} = p_{2} = 1

both terms with

e^{{(k)}^{2}}

and

e^{{(k)}^{3}}

cancel out while the term

e^{{(k)}^{4}}

is as follows

- A_{2} M_{1}^{2} + C_{1} M_{1} - M_{1}^{3} .

Let us remember that

C_{1} = 2 A_{2} M_{1} + α A_{2} F^{'} (ξ) M_{1}

, so that

A_{2} M_{1}^{2} + α A_{2} F^{'} (ξ) 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} = 2 A_{2} + α A_{2} F^{'} (ξ)

. Given that

M_{1} = A_{2} + α A_{2} F^{'} (ξ)

, replacing

M_{1}

in equation (32), we get

(A_{2} + α A_{2} F^{'} (ξ)) M_{1}^{2} - M_{1}^{3} \to M_{1}^{3} - M_{1}^{3} = 0,

resulting in the error equation

\begin{matrix} e^{(k + 1)} - ξ & = (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} \\ - P^{- 1} M_{1}^{2} Q) M_{1}) e^{{(k)}^{5}} + O (e^{{(k)}^{6}}) . \end{matrix}

(33)

From the above it is clear to see that if

p_{1} = 1

but

p_{2} \neq 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

\frac{1}{3} n^{3} + n^{2} - \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

M S (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

T r a u b - S t e

, while those with intermediate costs are

O s t r o_{01}

and

M S (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^{\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} (x) = {(f_{1}^{1} (x), f_{2}^{1} (x), \dots, f_{25}^{1} (x))}^{T}$ , where

$\begin{matrix} f_{i}^{1} (x) & = x_{i}^{2} x_{i + 1} - 1, i = 1, 2, \dots, 24, \\ f_{25}^{1} (x) & = x_{25}^{2} x_{1} - 1 . \end{matrix}$

whose solution is $\bar{ξ} = (1, 1, 1, \dots, 1)$ .
We also have $F^{2} (x) = {(f_{1}^{2} (x), f_{2}^{2} (x), \dots, f_{8}^{2} (x))}^{T}$ , where

$f_{i}^{2} (x) = x_{i} - cos (2 x_{i} - \sum_{k = 1}^{8} x_{k}), i = 1, 2, \dots, 8,$

being $\bar{ξ} \approx (0.5149, 0.5149, \dots, 0.5149)$ its solution.
We also consider $F^{3} (x) = {(f_{1}^{3} (x), f_{2}^{3} (x), \dots, f_{5}^{3} (x))}^{T}$ , where

$f_{i}^{3} (x) = \sum_{k = 1}^{5} (x_{k}) - x_{i} - e^{(- x_{i})}, i = 1, 2, 3, 4, 5,$

where the solution is $\bar{ξ} \approx (0.20389, 0.20389, 0.20389, \dots, 0.20389)$ .
We test also $F^{4} (x) = {(f_{1}^{4} (x), f_{2}^{4} (x), \dots, f_{5}^{4} (x))}^{T}$ , where

$f_{i}^{4} (x) = \sum_{k = 1}^{5} (x_{k}) - x_{i} - e^{(- x_{i})} x_{i}, i = 1, 2, 3, 4, 5,$

whose solution is $\bar{ξ} = (0, 0, 0, \dots, 0)$ .
$F^{5} (x) = {(f_{1}^{5} (x), f_{2}^{5} (x), \dots, f_{10}^{5} (x))}^{T}$ , where

$f_{i}^{5} (x) = x_{i} + 1 - 2 log (1 - x_{i} + \sum_{k = 1}^{10} (x_{k})), i = 1, 2, 3, \dots, 10,$

being it solution $\bar{ξ} \approx (7.4370, 7.4370, \dots, 7.4370)$ .
$F^{6} (x) = {(f_{1}^{6} (x), f_{2}^{6} (x), \dots, f_{5}^{6} (x))}^{T}$ , where

$f_{i}^{6} (x) = x_{i} + 1.5 sin (\sum_{k = 1}^{5} (x_{k}) - x_{i}), i = 1, 2, 3, 4, 5,$

and there exist two solutions, ${\bar{ξ}}_{1} \approx (- 0.3004, - 0.3004, \dots, - 0.3004)$ , and ${\bar{ξ}}_{2} \approx (0.4579, 0.4579, \dots, 0.4579)$ .
$F^{7} (x) = {(f_{1}^{7} (x), f_{2}^{7} (x))}^{T}$ , where

$F_{i}^{7} (x) = arctan (x_{i}) + 1 - 2 [(\sum_{k = 1}^{2} x_{k}^{2}) - x_{i}^{2}] = 0, i = 1, 2,$

being $\bar{ξ} \approx (0.936, 0.936)$ .
We also test $F^{8} (x) = {(f_{1}^{8} (x), f_{2}^{8} (x))}^{T}$ , where

$\begin{matrix} f_{1}^{8} (x) & = l o g (| x_{1} |) + | x_{2} |, \\ f_{2}^{8} (x) & = e^{(x_{1})} + x_{2} - 1, \end{matrix}$

with solutions ${\bar{ξ}}_{1} \approx (- 0.6275, 0.4661)$ and ${\bar{ξ}}_{2} \approx (0.5122, - 0.669)$ .

Numerical results have been obtained with Matlab2022b version, using 8000 digits in variable precision arithmetics, a processor AMD

A 12 - 9720 P

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).
$\bar{ξ}$ : obtained solution.
Cpu-time: average time in seconds required by the iterative method to reach the solution of the problem when executed ten times.
$ρ$ : approximated computational order of convergence, ACOC, firstly appeaaring in [19]

$ρ = \frac{ln \frac{∥x^{(k + 1)} - x^{(k)}∥}{∥x^{(k)} - x^{(k - 1)}∥}}{ln \frac{∥x^{(k)} - x^{(k - 1)}∥}{∥x^{(k - 1)} - x^{(k - 2)}∥}}, k = 2, 3, \dots,$

(if $ρ$ is not stable, then "-" appears in the table).
$ϵ_{a p r o x} = ∥x^{(k + 1)} - x^{(k)}∥$ .
$ϵ_{f} = ∥F (x^{(k + 1)})∥$ . If $ϵ_{f}$ or $ϵ_{a p r o x}$ 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 (x^{(k + 1)})∥ < 10^{- 100}$ ,
(ii): $∥x^{(k + 1)} - x^{(k)}∥ < 10^{- 100}$ ,
(iii): 50 iterations are reached.

In Table 2,

O s t r 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

M S (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

M S (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

M S (p 1, p 2)

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

M S (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

M S (p 1, p 2)

. On the other hand, the best errors are obtained by

O s t r o_{01}

scheme, which has a competitive computational time considering it requires four iterations, similar to

T r a u b - S t e

. The latter has shown to be the one that requires the most time to converge.

In Table 7, the best errors were obtained by

M S (1, - 1)

and

T r a u b - S t e

, which converged to different solutions, while in terms of performance,

M S (1, 1)

and

O s t r o_{01}

appeared to be better.

In Table 8,

O s t r o_{01}

method could not converge to the solution because it reached the maximum number of iterations, whereas the

M S (p 1, p 2)

and

M_{4, 3}

methods showed better overall behavior.

In the Table 9,

M S (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

T r a u b - S t e

method is less efficient, with a computational time of

3.1954

but with a larger errors.

O s t r 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

M S

and

O s t r o 01

methods.

In Table 10, method

M S (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

T r a u b - S t e

method, despite having the same convergence order as

M S (1, - 1)

, shows a slightly higher computational time and larger errors in the difference between iterates.

O s t r 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

M S (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

T r a u b - S t e

scheme, although it has the same convergence order as

M S (1, - 1)

, shows a significantly higher computational time but better results in terms of errors.

M S (1, - 1)

, while requiring more iterations than

M S (1, 1)

, still maintains a competitive time compared to

T r a u b - S t e

. 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.

M S (1, 1)

stands out as the most efficient in terms of time and iterations, while

M S (1, - 1)

yields the best errors. The method

O s t r 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

O s t r 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.

Author Contributions

Conceptualization, J.R.T. and A.C.; software, A.R.-C.; methodology, J.G.M.; validation, J.G.M., A.C. and J.R.T.; investigation, A.C.; formal analysis, J.R.; resources, A.R.-C.; writing—original draft preparation, A.R.-C.; visualization, J.G.M.; writing—review and editing, J.R.T. and A.C.; supervision, J.R.T. and A.C. All the authors have read and agreed to the published version of this paper.

Conflicts of Interest

The authors declare no conflict of interest.

References

Samanskii, V.: On a modification of the Newton method. Ukr. Math. J. 19, 133–138 (1967).
Steffensen, J.F.: Remarks on iteration. Skand. Aktuarietidskr. 1, 64–72 (1933). [CrossRef]
Ortega, J.M.; Rheinboldt, W.C. Iterative Solution of Nonlinear Equations in Several Variables; Academic Press: Cambridge, MA, USA, 1970. [CrossRef]
Chicharro, F.I., Cordero, A., Gutiérrez, J.M., Torregrosa, J.R.: Complex dynamics of derivative-free methods for nonlinear equations. Appl. Math. Comput. 219, 7023–7035 (2013). [CrossRef]
Cordero, A., Soleymani, F., Torregrosa, J.R., Shateyi, S.: Basins of attraction for various Steffensen-type methods. Appl. Math. 2014, Article ID 539707 1–17 (2014). [CrossRef]
Cordero, A., García-Maimó, J., Torregrosa, J.R., Vassileva, M.P.: Solving nonlinear problems by Ostrowski–Chun type parametric families. Math. Chem. 53, 430–449 (2015). [CrossRef]
Jarratt, P.: Some fourth order multipoint iterative methods for solving equations. Math. Comp. 20, 434–437 (1966). [CrossRef]
Sharma, J.R., Arora, H.: On efficient weighted-Newton methods for solving systems of nonlinear equations. Appl. Math. Comput. 222, 497–506 (2013). [CrossRef]
Montazeri, H., Soleymani, F., Shateyi, S., Motsa, S.S.: On a new method for computing the numerical solution of systems of nonlinear equations. Appl. Math. 2012 ID. 751975, 1–15 (2012). [CrossRef]
Sharma, J.R., Arora, H.: Efficient derivative-free numerical methods for solving systems of nonlinear equations. Comp. Appl. Math. 35, 269–284 (2016). [CrossRef]
Amiri, A.R., Cordero, A., Darvishi, M.T., Torregrosa, J.R.: Preserving the order of convergence: Low-complexity Jacobian-free iterative schemes for solving nonlinear systems. Comput. Appl. Math. 337, 87–97 (2018). [CrossRef]
Behl, R., Cordero, A., Torregrosa, J. R., Bhalla, S., A New High-Order Jacobian-Free Iterative Method with Memory for Solving Nonlinear Systems. Mathematics 9(17), 2122 (2021). [CrossRef]
Ostrowski, A.M., Solutions of Equations and Systems of Equations, Academic Press, New York, London, 1966.
Abad, M.F., Cordero, A., Torregrosa, J.R., A family of seventh-order schemes for solving nonlinear systems, Bull. Math. Soc. Sci. Math. Roumanie 57(105)(2), 133–145 (2014).
Cordero, A. , Maimó, J.G. ,Torregrosa, J.R., Vassileva, M.P., Solving nonlinear problems by Ostrowski-Chun type parametric families, Math. Chem. 53, 430–449 (2015). [CrossRef]
Sharma, J.R., Arora, H., On efficient weighted-Newton methods for solving systems of nonlinear equations, Applied Mathematics and Computation 222, 497–506 (2013). [CrossRef]
Cordero, A, Rojas–Hiciano, R. V., Torregrosa, J. R., Vassileva, M. P: A highly efficient class of optimal fourth-order methods for solving nonlinear systems. Numerical Algorithms 95, 1879–1904 (2024). [CrossRef]
Singh, H., Sharma, J. R., Kumar, S.: A simple yet efficient two-step fifth-order weighted-Newton method for nonlinear models. Numerical Algorithms 93, 203–225 (2023). [CrossRef]
Cordero, A., Torregrosa, J.R., Variants of Newton’s method using fifth order quadrature formulas. Appl. Math. Comput. 190, 686–698 (2007). [CrossRef]
Traub, I.F. Iterative Methods for the Solution of Equations; Prentice-Hall: Englewood Cliffs, NJ, USA, 1964. [CrossRef]

Figure 1. I indices for

M S (p_{1}, p_{2})

and comparison methods

Figure 1. I indices for

M S (p_{1}, p_{2})

and comparison methods

Figure 2. I indices for

M S (p_{1}, p_{2})

and comparison schemes

Figure 2. I indices for

M S (p_{1}, p_{2})

and comparison schemes

Table 1. Computational effort of new and comparison schemes

Method	Complexity C
$M S (p_{1}, p_{2})$	$\frac{2}{3} n^{3} + 6 n^{2} - \frac{2}{3} n$
$T r a u b - S t e$	$n^{3} + 10 n^{2} - 2 n$
$O s t r o_{01}$	$\frac{2}{3} n^{3} + 6 n^{2} + \frac{1}{3} n$
$M_{4, 3}$	$\frac{1}{3} n^{3} + 8 n^{2} + \frac{2}{3} n$
$M_{6, 3}$	$\frac{1}{3} n^{3} + 11 n^{2} + \frac{5}{3} n$

Table 2. Results for function

F^{1}

, using as seed

x^{(0)} = (1.5, 1.5, 1.5, \dots, 1.5)

Table 2. Results for function

F^{1}

, using as seed

x^{(0)} = (1.5, 1.5, 1.5, \dots, 1.5)

Iterative method	k	$ρ$	$ϵ_{a p r o x}$	$ϵ_{f}$	$\bar{ξ}$	Cpu-time
$M S (1, - 1)$	7	4.00	1.108e-51	4.624e-204	${\bar{ξ}}_{1}$	191.7242
$M S (1, 1)$	5	4.97	4.904e-24	3.995e-117	${\bar{ξ}}_{1}$	143.1027
$T r a u b - S t e$	5	4.00	1.241e-35	1.517e-140	${\bar{ξ}}_{1}$	227.8747
$O s t r o_{01}$	-	-	-	-	-	-
$M_{4, 3}$	6	4.00	7.338e-40	3.015e-158	${\bar{ξ}}_{1}$	163.9808
$M_{6, 3}$	5	5.96	1.336e-47	7.876e-284	${\bar{ξ}}_{1}$	142.6446

Table 3. Numerical results for

F^{2}

x^{(0)} = (1, 1, 1, \dots, 1)

Table 3. Numerical results for

F^{2}

x^{(0)} = (1, 1, 1, \dots, 1)

Iterative method	k	$ρ$	$ϵ_{a p r o x}$	$ϵ_{f}$	$\bar{ξ}$	Cpu-time
$M S (1, - 1)$	8	3.99	2.794e-29	4.092e-115	${\bar{ξ}}_{1}$	73.3126
$M S (1, 1)$	8	5.00	7.534e-74	2.162e-366	${\bar{ξ}}_{1}$	72.6143
$T r a u b - S t e$	16	4.00	9.389e-68	3.459e-269	${\bar{ξ}}_{1}$	317.9304
$O s t r o_{01}$	15	4.00	6.341e-70	1.641e-278	${\bar{ξ}}_{1}$	139.3808
$M_{4, 3}$	7	4.00	8.510e-87	3.949e-346	${\bar{ξ}}_{1}$	76.9114
$M_{6, 3}$	6	6.02	4.223e-47	8.063e-281	${\bar{ξ}}_{1}$	79.7898

Table 4. Numerical results for

F^{3}

x^{(0)} = (0.5, 0.5, 0.5, 0.5, 0.5)

Table 4. Numerical results for

F^{3}

x^{(0)} = (0.5, 0.5, 0.5, 0.5, 0.5)

Iterative method	k	$ρ$	$ϵ_{a p r o x}$	$ϵ_{f}$	$\bar{ξ}$	Cpu-time
$M S (1, - 1)$	4	4.00	1.305e-55	2.980e-221	${\bar{ξ}}_{1}$	12.8166
$M S (1, 1)$	4	5.00	4.997e-101	1.143e-503	${\bar{ξ}}_{1}$	12.6105
$T r a u b - S t e$	4	4.00	9.461e-68	1.372e-270	${\bar{ξ}}_{1}$	18.9457
$O s t r o_{01}$	4	4.00	1.420e-47	2.882e-189	${\bar{ξ}}_{1}$	12.4768
$M_{4, 3}$	4	4.00	3.440e-44	1.007e-175	${\bar{ξ}}_{1}$	12.8707
$M_{6, 3}$	3	6.07	1.482e-21	3.012e-127	${\bar{ξ}}_{1}$	10.5536

Table 5. Numerical results for

F^{4}

x^{(0)} = (0.5, 0.5, 0.5, 0.5, 0.5)

Table 5. Numerical results for

F^{4}

x^{(0)} = (0.5, 0.5, 0.5, 0.5, 0.5)

Iterative method	k	$ρ$	$ϵ_{a p r o x}$	$ϵ_{f}$	$\bar{ξ}$	Cpu-time
$M S (1, - 1)$	4	3.99	7.838e-31	4.801e-121	${\bar{ξ}}_{1}$	14.6365
$M S (1, 1)$	4	5.00	2.864e-52	2.514e-258	${\bar{ξ}}_{1}$	14.6532
$T r a u b - S t e$	4	4.00	8.656e-34	3.124e-133	${\bar{ξ}}_{1}$	21.7138
$O s t r o_{01}$	4	4.00	2.268e-35	6.443e-140	${\bar{ξ}}_{1}$	14.0140
$M_{4, 3}$	4	4.00	2.594e-47	9.227e-188	${\bar{ξ}}_{1}$	14.8093
$M_{6, 3}$	3	5.32	6.593e-26	7.254e-153	${\bar{ξ}}_{1}$	12.0667

Table 6. Numerical results for

F^{5}

x^{(0)} = (7, 7, 7, \dots, 7)

Table 6. Numerical results for

F^{5}

x^{(0)} = (7, 7, 7, \dots, 7)

Iterative method	k	$ρ$	$ϵ_{a p r o x}$	$ϵ_{f}$	$\bar{ξ}$	Cpu-time
$M S (1, - 1)$	3	4.03	1.343e-24	1.076e-101	${\bar{ξ}}_{1}$	58.6649
$M S (1, 1)$	3	5.03	5.715e-36	1.098e-183	${\bar{ξ}}_{1}$	58.6010
$T r a u b - S t e$	4	4.00	2.252e-97	1.397e-392	${\bar{ξ}}_{1}$	138.3238
$O s t r o_{01}$	4	4.00	1.317e-99	1.705e-401	${\bar{ξ}}_{1}$	76.6871
$M_{4, 3}$	3	4.00	1.490e-24	2.175e-101	${\bar{ξ}}_{1}$	64.3101
$M_{6, 3}$	3	6.01	1.565e-53	4.782e-326	${\bar{ξ}}_{1}$	69.4223

Table 7. Numerical results for

F^{6}

x^{(0)} = (0.75, 0.75, 0.75, 0.75, 0.75)

Table 7. Numerical results for

F^{6}

x^{(0)} = (0.75, 0.75, 0.75, 0.75, 0.75)

Iterative method	k	$ρ$	$ϵ_{a p r o x}$	$ϵ_{f}$	$\bar{ξ}$	Cpu-time
$M S (1, - 1)$	6	4.00	5.425e-88	5.270e-349	${\bar{ξ}}_{1}$	20.2517
$M S (1, 1)$	5	4.97	5.922e-39	1.965e-190	${\bar{ξ}}_{1}$	16.9191
$T r a u b - S t e$	20	4.00	9.536e-99	4.618e-391	${\bar{ξ}}_{2}$	154.5691
$O s t r o_{01}$	5	4.00	1.050e-33	2.164e-132	${\bar{ξ}}_{1}$	16.5869
$M_{4, 3}$	10	4.01	3.127e-35	1.500e-138	${\bar{ξ}}_{1}$	51.2143
$M_{6, 3}$	17	5.90	6.480e-23	2.076e-132	${\bar{ξ}}_{2}$	105.4212

Table 8. Numerical results for

F^{7}

x^{(0)} = (0.25, 0.25)

Table 8. Numerical results for

F^{7}

x^{(0)} = (0.25, 0.25)

Iterative method	k	$ρ$	$ϵ_{a p r o x}$	$ϵ_{f}$	$\bar{ξ}$	Cpu-time
$M S (1, - 1)$	5	4.00	5.203e-41	6.322e-161	${\bar{ξ}}_{1}$	3.6588
$M S (1, 1)$	4	5.00	5.282e-73	6.419e-362	${\bar{ξ}}_{1}$	3.0030
$T r a u b - S t e$	6	4.00	6.615e-81	7.658e-321	${\bar{ξ}}_{1}$	4.6891
$O s t r o_{01}$	-	-	-	-	-	-
$M_{4, 3}$	6	4.00	5.613e-98	3.787e-389	${\bar{ξ}}_{1}$	3.7011
$M_{6, 3}$	7	5.92	2.032e-32	5.723e-190	${\bar{ξ}}_{1}$	4.4699

Table 9. Numerical results for

F^{8}

x^{(0)} = (0.25, 0.25)

Table 9. Numerical results for

F^{8}

x^{(0)} = (0.25, 0.25)

Iterative method	k	$ρ$	$ϵ_{a p r o x}$	$ϵ_{f}$	$\bar{ξ}$	Cpu-time
$M S (1, - 1)$	4	4.00	1.139e-54	9.420e-217	${\bar{ξ}}_{1}$	2.9146
$M S (1, 1)$	4	4.20	1.415e-72	2.272e-301	${\bar{ξ}}_{1}$	2.8907
$T r a u b - S t e$	4	4.11	4.417e-48	1.391e-191	${\bar{ξ}}_{1}$	3.1954
$O s t r o_{01}$	6	4.00	2.520e-78	1.241e-310	${\bar{ξ}}_{1}$	3.4855
$M_{4, 3}$	7	4.02	1.401e-72	3.093e-288	${\bar{ξ}}_{1}$	4.0501
$M_{6, 3}$	6	6.27	4.634e-58	7.902e-346	${\bar{ξ}}_{1}$	3.6475

Table 10. Numerical results for

F^{8}

x^{(0)} = (1.25, 1.25)

Table 10. Numerical results for

F^{8}

x^{(0)} = (1.25, 1.25)

Iterative method	k	$ρ$	$ϵ_{a p r o x}$	$ϵ_{f}$	$\bar{ξ}$	Cpu-time
$M S (1, - 1)$	6	4.00	1.019e-65	6.029e-261	${\bar{ξ}}_{1}$	4.1403
$M S (1, 1)$	5	4.05	2.474e-42	4.929e-177	${\bar{ξ}}_{1}$	3.3999
$T r a u b - S t e$	6	4.00	8.946e-86	2.341e-342	${\bar{ξ}}_{1}$	4.4433
$O s t r o_{01}$	8	4.00	2.050e-55	5.429e-219	${\bar{ξ}}_{1}$	4.4246
$M_{4, 3}$	-	-	-	-	-	-
$M_{6, 3}$	-	-	-	-	-	-

Table 11. Numerical results for

F^{8}

x^{(0)} = (2, 2)

Table 11. Numerical results for

F^{8}

x^{(0)} = (2, 2)

Iterative method	k	$ρ$	$ϵ_{a p r o x}$	$ϵ_{f}$	$\bar{ξ}$	Cpu-time
$M S (1, - 1)$	10	4.00	6.679e-71	1.113e-281	${\bar{ξ}}_{1}$	6.9568
$M S (1, 1)$	8	4.28	4.387e-77	3.113e-321	${\bar{ξ}}_{1}$	5.4126
$T r a u b - S t e$	15	4.00	4.978e-100	2.442e-397	${\bar{ξ}}_{2}$	10.6624
$O s t r o_{01}$	-	-	-	-	-	-
$M_{4, 3}$	-	-	-	-	-	-
$M_{6, 3}$	-	-	-	-	-	-

Table 12. Numerical results for

F^{8}

x^{(0)} = (2.25, 2.25)

Table 12. Numerical results for

F^{8}

x^{(0)} = (2.25, 2.25)

Iterative method	k	$ρ$	$ϵ_{a p r o x}$	$ϵ_{f}$	$\bar{ξ}$	Cpu-time
$M S (1, - 1)$	16	4.00	2.821e-94	3.543e-375	${\bar{ξ}}_{1}$	10.5084
$M S (1, 1)$	6	4.19	8.959e-45	1.550e-184	${\bar{ξ}}_{2}$	4.2948
$T r a u b - S t e$	-	-	-	-	-	-
$O s t r o_{01}$	$k > 50$	-	-	-	-	-
$M_{4, 3}$	-	-	-	-	-	-
$M_{6, 3}$	-	-	-	-	-	-

Table 13. Numerical results for

F^{8}

x^{(0)} = (2.5, 2.5)

Table 13. Numerical results for

F^{8}

x^{(0)} = (2.5, 2.5)

Iterative method	k	$ρ$	$ϵ_{a p r o x}$	$ϵ_{f}$	$\bar{ξ}$	Cpu-time
$M S (1, - 1)$	10	4.00	2.063e-86	6.071e-344	${\bar{ξ}}_{2}$	6.8361
$M S (1, 1)$	7	4.17	3.281e-52	6.327e-218	${\bar{ξ}}_{1}$	4.7535
$T r a u b - S t e$	-	-	-	-	-	-
$O s t r o_{01}$	-	-	-	-	-	-
$M_{4, 3}$	-	-	-	-	-	-
$M_{6, 3}$	-	-	-	-	-	-

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

© 2024 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

Copyright: This open access article is published under a Creative Commons CC BY 4.0 license, which permit the free download, distribution, and reuse, provided that the author and preprint are cited in any reuse.

MDPI Initiatives

Important Links

Choose an area of interest and we will send you notifications of new preprints at your preferred frequency.

Disclaimer

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

Abstract

1. Introduction

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

3. Efficiency Analysis

4. Numerical Results

5. Conclusions

Author Contributions

Conflicts of Interest

References

MDPI Initiatives

Important Links

Subscribe