AWeight-Function Generalization of Singh’s Fifth-Order Method for Systems of Nonlinear Equations, with Application to a Discretized Stationary Viscous Burgers Problem

1. Introduction

The numerical solution of systems of nonlinear equations is a central topic in numerical analysis, applied mathematics, and scientific computing, since a broad range of mathematical models arising in physics, engineering, biology, and other applied sciences can be formulated as root-finding problems for vector-valued functions [1,2]. In this context, Newton’s method remains one of the most classical and influential techniques because of its simple structure and local quadratic convergence under suitable regularity assumptions. Nevertheless, its practical use may be severely affected by the need to evaluate and factorize the Jacobian matrix at each iteration, particularly when the dimension of the problem is large or when the nonlinear operator has a complicated structure [1,2]. These difficulties have motivated the construction of multipoint and high-order iterative methods whose purpose is to improve the local convergence behavior while maintaining a reasonable computational cost per iteration.

The design of predictor-corrector schemes provides a natural framework for the development of such methods. In particular, methods based on weighting functions have attracted considerable interest, as they introduce functional degrees of freedom into the correction stage, allowing for the cancellation of dominant terms in the error equation and, consequently, the construction of more accurate and flexible iterative families. This philosophy is closely related to the well-known Kung-Traub conjecture, according to which, for the scalar case, a memoryless multipoint method using n functional evaluations can achieve a maximum order of $2^{n-1}$ [3].

In the scalar case, the use of weight functions has proved to be a particularly fruitful strategy for generating high-order methods from simpler iterative schemes. Instead of prescribing a fixed correction term, one introduces a functional factor depending on suitable ratios or accelerators, thereby obtaining parametric families whose members can be tuned to satisfy desired convergence conditions. In this direction, Jaiswal proposed a class of fourth-order methods obtained by using weight functions, showing that this procedure provides a systematic mechanism for constructing efficient iterative families with improved local behavior [4]. Since then, the weight function technique has become one of the standard tools in the design of scalar high-order methods, especially when the goal is to increase flexibility without losing computational efficiency.

The extension of this strategy to systems of nonlinear equations is significantly more delicate. In the multidimensional setting, the scalar arguments based on ordinary Taylor expansions must be replaced by operator expansions and matrix-based error analysis, and the design of suitable weighted corrections becomes more involved. Even so, important advances have shown that this methodology can be successfully transferred to nonlinear systems. For instance, Artidiello, Cordero, Torregrosa, and Vassileva developed multidimensional generalizations of iterative schemes derived from the weight function procedure, establishing high-order convergence in the vector case and confirming the effectiveness of this approach beyond one-dimensional root-finding problems [5,6]. Later, Sharma et al, introduced efficient weighted-Newton methods for solving systems of nonlinear equations, further illustrating that weighted corrections can be adapted to the vector setting while preserving favorable convergence properties [7]. More recently, Capdevila, Cordero, and Torregrosa constructed a three-step family based on scalar and matrix weight functions, obtaining sixth-order convergence in both scalar and vector formulations [8]. Altogether, these contributions confirm that weight functions constitute a robust and versatile design tool for iterative schemes in nonlinear systems.

Within this line of research, fifth-order methods occupy an especially attractive position, since they provide a good compromise between local accuracy and computational effort. In particular, recent contributions such as those by Singh et al have emphasized the construction of simple and efficient fifth-order solvers for nonlinear systems, underlining the continuing relevance of this order in the search for practically competitive methods [9]. In this framework, Singh’s method represents a highly efficient reference scheme whose structure suggests that the correction stage can be generalized through appropriate weight functions. This observation is mathematically and computationally relevant: instead of working with a single fixed iterative formula, one may construct a whole family of methods that contains the original scheme as a particular case and, at the same time, preserves its essential convergence properties.

Motivated by these considerations, in this paper we introduce a weighted generalization of Singh’s method for solving systems of nonlinear equations. The proposed family is constructed by incorporating suitable weight functions into the iterative process in such a way that the original Singh scheme is recovered as a particular case for specific choices of these functions. The corresponding local convergence analysis makes it possible to determine explicit conditions on the weight functions under which the fifth-order convergence of the original method is preserved. In this way, the new family extends a highly efficient reference scheme without modifying its essential convergence structure, while enlarging the class of admissible iterative formulations and providing additional flexibility for the design of new solvers.

From the viewpoint of applications, nonlinear systems arising from the discretization of differential equations constitute a natural and demanding benchmark for high-order iterative methods. Among them, Burgers-type equations play a prominent role in applied mathematics, since they arise as prototype nonlinear models in fluid mechanics and related transport phenomena, and their stationary forms lead, after spatial discretization, to nonlinear algebraic systems that are well suited for testing iterative solvers [10,11]. In particular, the stationary viscous Burgers equation provides a representative nonlinear boundary-value problem in which convection and diffusion interact in a nontrivial way, giving rise to algebraic systems whose efficient numerical resolution is of clear interest [10]. For this reason, in order to illustrate the practical applicability of the proposed family, we consider a nonlinear system obtained from the discretization of a stationary Burgers equation.

The main contributions of this paper can be summarized as follows. First, we construct a new weighted family of fifth-order iterative methods for systems of nonlinear equations, obtained as a generalization of a highly efficient Singh-type scheme. Second, we derive sufficient conditions on the weight functions that guarantee the preservation of fifth-order local convergence. Third, we provide several admissible choices of weight functions, which generate different particular members of the family and illustrate its flexibility from both the theoretical and computational points of view. Fourth, we carry out a comparative numerical study on several large-scale nonlinear systems in order to assess the practical behavior of the proposed schemes. Finally, we illustrate the applicability of the family to a nonlinear system of differential origin obtained from the discretization of a stationary Burgers equation. Therefore, the paper combines the theoretical construction of a weighted iterative family with a representative application that highlights its usefulness in nonlinear numerical models.

The remainder of the paper is organized as follows. In Section 2, we present the preliminary concepts and auxiliary results used throughout the manuscript. Section 3 introduces the proposed weighted family of iterative methods and establishes its local order of convergence. This section also includes complementary theoretical results, such as propositions, remarks, and particular cases derived from the general formulation. In Section 4, we present the different weight functions selected for the numerical tests, together with the comparison methods considered in the study. This section also includes the analysis of the computational efficiency index and the numerical experiments carried out on several large-scale nonlinear systems in order to evaluate the practical performance of the proposed schemes. In Section 5, we address the discretization and numerical solution of a particular stationary Burgers equation, which leads to a second-order nonlinear ordinary differential equation. Finally, Section 6 summarizes the main conclusions of the work.

2. Preliminary Concepts

Throughout this work, we assume that the nonlinear system possesses a solution $\xi \in \mathcal{D}$, such that

$$F\left(\xi \right)=0,F^{\prime }\left(\xi \right)\mathrm{is}\mathrm{nonsingular},$$

and that all derivatives required for the subsequent analysis exist and remain continuous in a neighborhood of $\xi$.

To formalize this assumption, we begin by recalling the notion of a simple root, which will be used throughout the local convergence analysis.

Definition 1.

A point $\xi \in \mathcal{D}$ is called a simple root of the nonlinear system $F\left(x\right)=0$ if

$$F\left(\xi \right)=0\mathrm{and}{F}^{\prime }\left(\xi \right)\mathrm{is}\mathrm{nonsingular}.$$

To establish the notation used in the local study, we now recall the multilinear structure of the higher-order derivatives of F. For $m\ge 1$, the m-th derivative of the map F at a point x is the m-linear operator

$$F^{\left(m\right)}\left(x\right):\underset{m\mathrm{slots}}{\underbrace{{\mathbb{R}}^n\times \dots \times {\mathbb{R}}^n}}⟶{\mathbb{R}}^n,F^{\left(m\right)}\left(x\right)\in {\mathcal{L}}_m\left({\mathbb{R}}^n;{\mathbb{R}}^n\right).$$

For brevity, we employ the compact notation

$$F^{\left(m\right)}\left(x\right)({\omega }_1,\dots ,{\omega }_m)=F^{\left(m\right)}\left(x\right){\omega }_1\dots {\omega }_m,F^{\left(m\right)}\left(x\right){\omega }^m:=F^{\left(m\right)}\left(x\right)(\omega ,\dots ,\omega ).$$

Moreover, whenever one argument of $F^{\left(m\right)}\left(x\right)$ is the vector $F^{\left(p\right)}\left(x\right){\omega }^p$, we write

$$F^{\left(m\right)}\left(x\right){\omega }^{m-1}\left(F^{\left(p\right)}\left(x\right){\omega }^p\right):=F^{\left(m\right)}\left(x\right)\left(\underset{m-1}{\underbrace{\omega ,\dots ,\omega }},F^{\left(p\right)}\left(x\right){\omega }^p\right).$$

In order to derive the local error equation of the proposed family, we recall the Taylor expansion of F around the solution $\xi$. Let $x=\xi +\eta$, with $\eta$ sufficiently small. Assuming that $F^{\prime }\left(\xi \right)$ is invertible, a Taylor expansion around $\xi$ yields

$$F(\xi +\eta )=F^{\prime }\left(\xi \right)\left[\eta +C_2{\eta }^2+C_3{\eta }^3+C_4{\eta }^4\right]+\mathcal{O}\left({\eta }^5\right),$$

(1)

where the tensors $C_j$ are given by

$$C_j={\displaystyle \frac{1}{j!}}{\left[F^{\prime }\left(\xi \right)\right]}^{-1}{F}^{\left(j\right)}\left(\xi \right),j\ge 2.$$

A similar expansion for the Jacobian operator of F at $\xi +\eta$ around $\xi$ reads

$$F^{\prime }(\xi +\eta )=F^{\prime }\left(\xi \right)\left[I+2C_2\eta +3C_3{\eta }^2+4C_4{\eta }^3\right]+\mathcal{O}\left({\eta }^4\right).$$

(2)

Moreover, the inverse matrix ${\left[F^{\prime }(\xi +\eta )\right]}^{-1}$ admits the Neumann-type expansion

$${\left[F^{\prime }(\xi +\eta )\right]}^{-1}=\left[I-2C_2\eta +(4C_2^2-3C_3){\eta }^2+(-8C_2^3+6C_3{C}_2+6C_2{C}_3-4C_4){\eta }^3\right]{\left[F^{\prime }\left(\xi \right)\right]}^{-1}+\mathcal{O}\left({\eta }^4\right).$$

(3)

Next, we introduce the local error notation associated with an iterative sequence.

Definition 2.

Given a sequence of approximations $\left\{x^{\left(k\right)}\right\}$ to ξ, we define the local error at iteration k by

$$e^{\left(k\right)}=x^{\left(k\right)}-\xi .$$

This notation allows us to express the local behavior of the method through its error equation.

Definition 3.

We say that the sequence $\left\{x^{\left(k\right)}\right\}$ converges locally to ξ with order $p\ge 1$ if

$$e^{(k+1)}=\mathcal{T}{e}^{{\left(k\right)}^p}+\mathcal{O}\left(e^{{\left(k\right)}^{p+1}}\right),$$

where $\mathcal{T}$ denotes a suitable p-linear operator characterizing the leading term of the local behavior.

For the numerical experiments, we also use the following computational approximation of the convergence order.

Definition 4

([12]). The approximated computational order of convergence (ACOC) is given by

$$p\approx ACOC={\displaystyle \frac{ln\left(\parallel x^{(k+1)}-x^{\left(k\right)}\parallel /\parallel x^{\left(k\right)}-x^{(k-1)}\parallel \right)}{ln\left(\parallel x^{\left(k\right)}-x^{(k-1)}\parallel /\parallel x^{(k-1)}-x^{(k-2)}\parallel \right)}}.$$

(4)

Since the proposed family is built by means of weighted corrections, it is convenient to make explicit the notion of weight function used in this work.

Definition 5.

A weight function is an auxiliary scalar-, vector-, or matrix-valued function introduced into an iterative scheme in order to modify the correction step and control the cancellation of dominant terms in the local error equation, thus improving the convergence order without significantly increasing the computational cost.

3. Local Convergence of the Fifth-Order Weighted Family

In this section, we analyze the local convergence of the proposed weighted family. Starting from its general iterative formulation, we derive sufficient conditions on the weight functions to ensure that the methods preserve fifth-order convergence in a neighborhood of a simple root. The corresponding error equation is also obtained, allowing us to identify the influence of the weight functions on the leading term and to justify several admissible particular cases.

Theorem 1.

Let $F:\Omega \subset {\mathbb{R}}^n\to {\mathbb{R}}^n$ be sufficiently Fréchet differentiable in a neighborhood of a simple root $\xi \in \Omega$, that is,

$$F\left(\xi \right)=0,F^{\prime }\left(\xi \right)\mathrm{is}\mathrm{invertible}.$$

Consider the iterative family

$$\left\{\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[F^{\prime }\left(y^{\left(k\right)}\right)\right]}^{-1}\left(H\left({\nu }_k\right)F\left(y^{\left(k\right)}\right)+G\left({\nu }_k\right)F\left(x^{\left(k\right)}\right)\right),\hfill \end{array}\right.$$

where

$${\nu }_k={\displaystyle \frac{{∥F\left(y^{\left(k\right)}\right)∥}^2}{{∥F\left(x^{\left(k\right)}\right)∥}^2}},$$

and H and G are analytic scalar weight functions in a neighborhood of 0.

Let

$$e^{\left(k\right)}=x^{\left(k\right)}-\xi ,C_q={\displaystyle \frac{1}{q!}}{\left[F^{\prime }\left(\xi \right)\right]}^{-1}{F}^{\left(q\right)}\left(\xi \right),q\ge 2,$$

and define

$$P={\left(F^{\prime }\left(\xi \right)\right)}^T{F}^{\prime }\left(\xi \right),Q={\displaystyle \sum _{i=1}^n}\left(R_i^T{S}_i+S_i^T{R}_i\right),$$

where

$$R_i=f_i^{\prime }\left(\xi \right),S_i={\displaystyle \frac{1}{2}}{f}_i^{\prime \prime }\left(\xi \right),i=1,\dots ,n.$$

If the weight functions satisfy

$$H\left(0\right)=1,H^{\prime }\left(0\right)=1,G\left(0\right)=0,G^{\prime }\left(0\right)=0,$$

then the iterative family converges locally to ξ with order five. Moreover, its error equation is

$$e^{(k+1)}=\left[\left(2-{\displaystyle \frac{1}{2}}{G}^{\prime \prime }\left(0\right)\right)C_2^4-2C_3{C}_2^2+P^{-1}Q{C}_2^3\right]e^{{\left(k\right)}^5}+\mathcal{O}\left(e^{{\left(k\right)}^6}\right).$$

In particular, the coefficients of the terms of orders $e^{\left(k\right)}$, $e^{{\left(k\right)}^2}$, $e^{{\left(k\right)}^3}$, and $e^{{\left(k\right)}^4}$ vanish under the previous conditions. In the multidimensional setting, the operators $C_q$ do not commute in general, so the order of the products must be preserved throughout the analysis.

Proof. 

We preserve the compact multilinear notation introduced in Section 2, and we keep the order of the products throughout the proof, since in the multidimensional setting the operators $C_q$ do not commute in general.

The Taylor expansions of $F\left(x^{\left(k\right)}\right)$ and $F^{\prime }\left(x^{\left(k\right)}\right)$ around $\xi$ are

$$F\left(x^{\left(k\right)}\right)=F^{\prime }\left(\xi \right)\left[e^{\left(k\right)}+C_2{e}^{{\left(k\right)}^2}+C_3{e}^{{\left(k\right)}^3}+C_4{e}^{{\left(k\right)}^4}+C_5{e}^{{\left(k\right)}^5}\right]+\mathcal{O}\left(e^{{\left(k\right)}^6}\right),$$

(5)

and

$$F^{\prime }\left(x^{\left(k\right)}\right)=F^{\prime }\left(\xi \right)+F^{\prime \prime }\left(\xi \right)e^{\left(k\right)}+{\displaystyle \frac{1}{2!}}{F}^{\left(3\right)}\left(\xi \right)e^{{\left(k\right)}^2}+{\displaystyle \frac{1}{3!}}{F}^{\left(4\right)}\left(\xi \right)e^{{\left(k\right)}^3}+{\displaystyle \frac{1}{4!}}{F}^{\left(5\right)}\left(\xi \right)e^{{\left(k\right)}^4}+\mathcal{O}\left(e^{{\left(k\right)}^5}\right).$$

Factoring $F^{\prime }\left(\xi \right)$, and using

$$C_q={\displaystyle \frac{1}{q!}}{\left[F^{\prime }\left(\xi \right)\right]}^{-1}{F}^{\left(q\right)}\left(\xi \right),q\ge 2,$$

we obtain

$$F^{\prime }\left(x^{\left(k\right)}\right)=F^{\prime }\left(\xi \right)\left[I+2C_2{e}^{\left(k\right)}+3C_3{e}^{{\left(k\right)}^2}+4C_4{e}^{{\left(k\right)}^3}+5C_5{e}^{{\left(k\right)}^4}\right]+\mathcal{O}\left(e^{{\left(k\right)}^5}\right).$$

(6)

Now, let us calculate the Taylor expansion of ${\left[F^{\prime }\left(x^{\left(k\right)}\right)\right]}^{-1}$. For this purpose, we use the identity

$${\left[F^{\prime }\left(x^{\left(k\right)}\right)\right]}^{-1}{F}^{\prime }\left(x^{\left(k\right)}\right)=I.$$

(7)

Thus, we assume that

$${\left[F^{\prime }\left(x^{\left(k\right)}\right)\right]}^{-1}=\left[I+X_2{e}^{\left(k\right)}+X_3{e}^{{\left(k\right)}^2}+X_4{e}^{{\left(k\right)}^3}+X_5{e}^{{\left(k\right)}^4}\right]{\left[F^{\prime }\left(\xi \right)\right]}^{-1}+\mathcal{O}\left(e^{{\left(k\right)}^5}\right).$$

Substituting into (7), and comparing coefficients of like powers of $e^{\left(k\right)}$, we obtain

$$X_2=-2C_2,$$

$$X_3=4C_2^2-3C_3,$$

$$X_4=-8C_2^3+6C_3{C}_2+6C_2{C}_3-4C_4,$$

$$X_5=16C_2^4-12C_3{C}_2^2-12C_2{C}_3{C}_2-12C_2^2{C}_3+9C_3^2+8C_4{C}_2+8C_2{C}_4-5C_5.$$

Therefore,

$${\displaystyle \begin{array}{cc}\hfill {\left[F^{\prime }\left(x^{\left(k\right)}\right)\right]}^{-1}=& [I-2C_2{e}^{\left(k\right)}+(4C_2^2-3C_3)e^{{\left(k\right)}^2}+\left(-8C_2^3+6C_3{C}_2+6C_2{C}_3-4C_4\right)e^{{\left(k\right)}^3}\hfill \\ & +\left(16C_2^4-12C_3{C}_2^2-12C_2{C}_3{C}_2-12C_2^2{C}_3+9C_3^2+8C_4{C}_2+8C_2{C}_4-5C_5\right)e^{{\left(k\right)}^4}\hfill \\ & ]{\left[F^{\prime }\left(\xi \right)\right]}^{-1}+\mathcal{O}\left(e^{{\left(k\right)}^5}\right).\hfill \end{array}}$$

Multiplying by (5), we get

$$\begin{array}{cc}\hfill {\left[F^{\prime }\left(x^{\left(k\right)}\right)\right]}^{-1}F\left(x^{\left(k\right)}\right)=& e^{\left(k\right)}-C_2{e}^{{\left(k\right)}^2}+\left(2C_2^2-2C_3\right)e^{{\left(k\right)}^3}\hfill \\ & +\left(-4C_2^3+3C_3{C}_2+4C_2{C}_3-3C_4\right)e^{{\left(k\right)}^4}\hfill \\ & +\left(8C_2^4-6C_3{C}_2^2-6C_2{C}_3{C}_2-8C_2^2{C}_3+6C_3^2+4C_4{C}_2+6C_2{C}_4-4C_5\right)e^{{\left(k\right)}^5}\hfill \\ & +\mathcal{O}\left(e^{{\left(k\right)}^6}\right).\hfill \end{array}$$

Hence,

$$\begin{array}{cc}\hfill e_y^{\left(k\right)}=y^{\left(k\right)}-\xi =& C_2{e}^{{\left(k\right)}^2}+\left(2C_3-2C_2^2\right)e^{{\left(k\right)}^3}\hfill \\ \hfill & +\left(4C_2^3-3C_3{C}_2-4C_2{C}_3+3C_4\right)e^{{\left(k\right)}^4}\hfill \\ \hfill & +\left(-8C_2^4+6C_3{C}_2^2+6C_2{C}_3{C}_2+8C_2^2{C}_3-6C_3^2-4C_4{C}_2-6C_2{C}_4+4C_5\right)e^{{\left(k\right)}^5}\hfill \\ \hfill & +\mathcal{O}\left(e^{{\left(k\right)}^6}\right).\hfill \end{array}$$

(8)

Now, let us compute $F\left(y^{\left(k\right)}\right)$. Since $e_y^{\left(k\right)}=\mathcal{O}\left(e^{{\left(k\right)}^2}\right)$, up to order $e^{{\left(k\right)}^5}$ we only need the terms

$$F\left(y^{\left(k\right)}\right)=F^{\prime }\left(\xi \right)\left[e_y^{\left(k\right)}+C_2{\left(e_y^{\left(k\right)}\right)}^2\right]+\mathcal{O}\left(e^{{\left(k\right)}^6}\right).$$

From (8),

$${\left(e_y^{\left(k\right)}\right)}^2=C_2^2{e}^{{\left(k\right)}^4}+\left(2C_2{C}_3+2C_3{C}_2-4C_2^3\right)e^{{\left(k\right)}^5}+\mathcal{O}\left(e^{{\left(k\right)}^6}\right).$$

Therefore,

$$\begin{array}{cc}\hfill F\left(y^{\left(k\right)}\right)=F^{\prime }\left(\xi \right)[& C_2{e}^{{\left(k\right)}^2}+\left(2C_3-2C_2^2\right)e^{{\left(k\right)}^3}\hfill \\ & +\left(5C_2^3-3C_3{C}_2-4C_2{C}_3+3C_4\right)e^{{\left(k\right)}^4}\hfill \\ & +\left(-12C_2^4+6C_3{C}_2^2+8C_2{C}_3{C}_2+10C_2^2{C}_3-6C_3^2-4C_4{C}_2-6C_2{C}_4+4C_5\right)e^{{\left(k\right)}^5}]\hfill \\ & +\mathcal{O}\left(e^{{\left(k\right)}^6}\right).\hfill \end{array}$$

Next, let us compute ${\left[F^{\prime }\left(y^{\left(k\right)}\right)\right]}^{-1}$. First,

$$F^{\prime }\left(y^{\left(k\right)}\right)=F^{\prime }\left(\xi \right)+F^{\prime \prime }\left(\xi \right)e_y^{\left(k\right)}+{\displaystyle \frac{1}{2!}}{F}^{\left(3\right)}\left(\xi \right){\left(e_y^{\left(k\right)}\right)}^2+\mathcal{O}\left(e^{{\left(k\right)}^5}\right).$$

Factoring $F^{\prime }\left(\xi \right)$, and using (8), we get

$$\begin{array}{cc}\hfill F^{\prime }\left(y^{\left(k\right)}\right)=& F^{\prime }\left(\xi \right)[I+2C_2^2{e}^{{\left(k\right)}^2}+\left(-4C_2^3+4C_2{C}_3\right)e^{{\left(k\right)}^3}\hfill \\ & +\left(8C_2^4-8C_2^2{C}_3-6C_2{C}_3{C}_2+3C_3{C}_2^2+6C_2{C}_4\right)e^{{\left(k\right)}^4}]+\mathcal{O}\left(e^{{\left(k\right)}^5}\right).\hfill \end{array}$$

Now, we assume that

$${\left[F^{\prime }\left(y^{\left(k\right)}\right)\right]}^{-1}=\left[I+Y_2{e}^{{\left(k\right)}^2}+Y_3{e}^{{\left(k\right)}^3}+Y_4{e}^{{\left(k\right)}^4}\right]{\left[F^{\prime }\left(\xi \right)\right]}^{-1}+\mathcal{O}\left(e^{{\left(k\right)}^5}\right).$$

Using again the identity

$${\left[F^{\prime }\left(y^{\left(k\right)}\right)\right]}^{-1}{F}^{\prime }\left(y^{\left(k\right)}\right)=I,$$

and comparing coefficients, we obtain

$$Y_2=-2C_2^2,Y_3=4C_2^3-4C_2{C}_3,$$

$$Y_4=-4C_2^4+8C_2^2{C}_3+6C_2{C}_3{C}_2-3C_3{C}_2^2-6C_2{C}_4.$$

Hence,

$$\begin{array}{cc}\hfill {\left[F^{\prime }\left(y^{\left(k\right)}\right)\right]}^{-1}=& [I-2C_2^2{e}^{{\left(k\right)}^2}+\left(4C_2^3-4C_2{C}_3\right)e^{{\left(k\right)}^3}\hfill \\ \hfill & +\left(-4C_2^4+8C_2^2{C}_3+6C_2{C}_3{C}_2-3C_3{C}_2^2-6C_2{C}_4\right)e^{{\left(k\right)}^4}]{\left[F^{\prime }\left(\xi \right)\right]}^{-1}+\mathcal{O}\left(e^{{\left(k\right)}^5}\right).\hfill \end{array}$$

(9)

We now analyze the quotient

$${\nu }_k={\displaystyle \frac{{∥F\left(y^{\left(k\right)}\right)∥}^2}{{∥F\left(x^{\left(k\right)}\right)∥}^2}}.$$

Using the notation

$$R_i=f_i^{\prime }\left(\xi \right),S_i={\displaystyle \frac{1}{2}}{f}_i^{\prime \prime }\left(\xi \right),i=1,\dots ,m,$$

we have

$$f_i\left(x^{\left(k\right)}\right)=R_i{e}^{\left(k\right)}+S_i{e}^{{\left(k\right)}^2}+\mathcal{O}\left(e^{{\left(k\right)}^3}\right),$$

$$f_i\left(y^{\left(k\right)}\right)=R_i{e}_y^{\left(k\right)}+S_i{\left(e_y^{\left(k\right)}\right)}^2+\mathcal{O}\left(e^{{\left(k\right)}^6}\right).$$

Since each $f_i$ is scalar-valued,

$$f_i^2\left(x^{\left(k\right)}\right)=R_i^T{R}_i{e}^{{\left(k\right)}^2}+\left(R_i^T{S}_i+S_i^T{R}_i\right)e^{{\left(k\right)}^3}+\mathcal{O}\left(e^{{\left(k\right)}^4}\right),$$

and

$$f_i^2\left(y^{\left(k\right)}\right)=R_i^T{R}_i{\left(e_y^{\left(k\right)}\right)}^2+\left(R_i^T{S}_i+S_i^T{R}_i\right){\left(e_y^{\left(k\right)}\right)}^3+\mathcal{O}\left(e^{{\left(k\right)}^8}\right).$$

Defining

$$P={\displaystyle \sum _{i=1}^m}{R}_i^T{R}_i,Q={\displaystyle \sum _{i=1}^m}\left(R_i^T{S}_i+S_i^T{R}_i\right),$$

it follows that

$${\nu }_k={\displaystyle \frac{P{\left(e_y^{\left(k\right)}\right)}^2+Q{\left(e_y^{\left(k\right)}\right)}^3+\mathcal{O}\left(e^{{\left(k\right)}^8}\right)}{P{e}^{{\left(k\right)}^2}+Q{e}^{{\left(k\right)}^3}+\mathcal{O}\left(e^{{\left(k\right)}^4}\right)}}.$$

Since

$${\left(e_y^{\left(k\right)}\right)}^2=C_2^2{e}^{{\left(k\right)}^4}+\left(2C_2{C}_3+2C_3{C}_2-4C_2^3\right)e^{{\left(k\right)}^5}+\mathcal{O}\left(e^{{\left(k\right)}^6}\right),$$

and ${\left(e_y^{\left(k\right)}\right)}^3=\mathcal{O}\left(e^{{\left(k\right)}^6}\right)$, we obtain

$${\nu }_k=C_2^2{e}^{{\left(k\right)}^2}+\left(2C_2{C}_3+2C_3{C}_2-4C_2^3-P^{-1}Q{C}_2^2\right)e^{{\left(k\right)}^3}+\mathcal{O}\left(e^{{\left(k\right)}^4}\right).$$

Next, we expand the weight functions. Let

$$h_j={\displaystyle \frac{H^{\left(j\right)}\left(0\right)}{j!}},g_j={\displaystyle \frac{G^{\left(j\right)}\left(0\right)}{j!}},j\ge 0.$$

Then

$$H\left({\nu }_k\right)=h_0+h_1{\nu }_k+h_2{\nu }_k^2+\mathcal{O}\left({\nu }_k^3\right),G\left({\nu }_k\right)=g_0+g_1{\nu }_k+g_2{\nu }_k^2+\mathcal{O}\left({\nu }_k^3\right).$$

Since

$${\nu }_k=\mathcal{O}\left(e^{{\left(k\right)}^2}\right),{\nu }_k^2=C_2^4{e}^{{\left(k\right)}^4}+\mathcal{O}\left(e^{{\left(k\right)}^5}\right),$$

we get

$$H\left({\nu }_k\right)=h_0+h_1{C}_2^2{e}^{{\left(k\right)}^2}+h_1\left(2C_2{C}_3+2C_3{C}_2-4C_2^3-P^{-1}Q{C}_2^2\right)e^{{\left(k\right)}^3}+\mathcal{O}\left(e^{{\left(k\right)}^4}\right),$$

because the term $h_2{\nu }_k^2$ contributes only from order $e^{{\left(k\right)}^6}$ onward when multiplied by ${\left[F^{\prime }\left(y^{\left(k\right)}\right)\right]}^{-1}F\left(y^{\left(k\right)}\right)$, and

$$\begin{array}{cc}\hfill G\left({\nu }_k\right)=& g_0+g_1{C}_2^2{e}^{{\left(k\right)}^2}+g_1\left(2C_2{C}_3+2C_3{C}_2-4C_2^3-P^{-1}Q{C}_2^2\right)e^{{\left(k\right)}^3}+g_2{C}_2^4{e}^{{\left(k\right)}^4}+\mathcal{O}\left(e^{{\left(k\right)}^5}\right),\hfill \end{array}$$

because now the quadratic term must be kept, since it is multiplied by ${\left[F^{\prime }\left(y^{\left(k\right)}\right)\right]}^{-1}F\left(x^{\left(k\right)}\right)$, whose leading term is of order $e^{\left(k\right)}$.

Using (9) and the previous expansion of $F\left(y^{\left(k\right)}\right)$, we obtain

$$\begin{array}{cc}\hfill {\left[F^{\prime }\left(y^{\left(k\right)}\right)\right]}^{-1}F\left(y^{\left(k\right)}\right)=& C_2{e}^{{\left(k\right)}^2}+\left(2C_3-2C_2^2\right)e^{{\left(k\right)}^3}\hfill \\ & +\left(3C_2^3-3C_3{C}_2-4C_2{C}_3+3C_4\right)e^{{\left(k\right)}^4}\hfill \\ & +\left(-4C_2^4+6C_3{C}_2^2+4C_2{C}_3{C}_2+6C_2^2{C}_3-6C_3^2-4C_4{C}_2-6C_2{C}_4+4C_5\right)e^{{\left(k\right)}^5}\hfill \\ & +\mathcal{O}\left(e^{{\left(k\right)}^6}\right),\hfill \end{array}$$

and similarly

$$\begin{array}{cc}\hfill {\left[F^{\prime }\left(y^{\left(k\right)}\right)\right]}^{-1}F\left(x^{\left(k\right)}\right)=& e^{\left(k\right)}+C_2{e}^{{\left(k\right)}^2}+\left(C_3-2C_2^2\right)e^{{\left(k\right)}^3}\hfill \\ & +\left(2C_2^3-4C_2{C}_3+C_4\right)e^{{\left(k\right)}^4}\hfill \\ & +\left(C_5+6C_2^2{C}_3+2C_2{C}_3{C}_2-3C_3{C}_2^2-6C_2{C}_4\right)e^{{\left(k\right)}^5}\hfill \\ & +\mathcal{O}\left(e^{{\left(k\right)}^6}\right).\hfill \end{array}$$

Therefore,

$$\begin{array}{cc}\hfill H\left({\nu }_k\right){\left[F^{\prime }\left(y^{\left(k\right)}\right)\right]}^{-1}F\left(y^{\left(k\right)}\right)=& h_0{C}_2{e}^{{\left(k\right)}^2}+h_0\left(2C_3-2C_2^2\right)e^{{\left(k\right)}^3}\hfill \\ & +\left[h_0\left(3C_2^3-3C_3{C}_2-4C_2{C}_3+3C_4\right)+h_1{C}_2^3\right]e^{{\left(k\right)}^4}\hfill \\ & +[h_0\left(-4C_2^4+6C_3{C}_2^2+4C_2{C}_3{C}_2+6C_2^2{C}_3-6C_3^2-4C_4{C}_2-6C_2{C}_4+4C_5\right)\hfill \\ & +h_1\left(2C_2^2{C}_3+2C_2{C}_3{C}_2+2C_3{C}_2^2-6C_2^4-P^{-1}Q{C}_2^3\right)]e^{{\left(k\right)}^5}\hfill \\ & +\mathcal{O}\left(e^{{\left(k\right)}^6}\right),\hfill \end{array}$$

and

$$\begin{array}{cc}\hfill G\left({\nu }_k\right){\left[F^{\prime }\left(y^{\left(k\right)}\right)\right]}^{-1}F\left(x^{\left(k\right)}\right)=& g_0{e}^{\left(k\right)}+g_0{C}_2{e}^{{\left(k\right)}^2}\hfill \\ & +\left[g_0\left(C_3-2C_2^2\right)+g_1{C}_2^2\right]e^{{\left(k\right)}^3}\hfill \\ & +\left[g_0\left(2C_2^3-4C_2{C}_3+C_4\right)+g_1\left(-3C_2^3+2C_2{C}_3+2C_3{C}_2-P^{-1}Q{C}_2^2\right)\right]e^{{\left(k\right)}^4}\hfill \\ & +[g_0\left(C_5+6C_2^2{C}_3+2C_2{C}_3{C}_2-3C_3{C}_2^2-6C_2{C}_4\right)\hfill \\ & +g_1\left(C_2^2{C}_3+2C_2{C}_3{C}_2+2C_3{C}_2^2-6C_2^4-P^{-1}Q{C}_2^3\right)+g_2{C}_2^4]e^{{\left(k\right)}^5}\hfill \\ & +\mathcal{O}\left(e^{{\left(k\right)}^6}\right).\hfill \end{array}$$

Finally, from the second step of the method,

$$e^{(k+1)}=e_y^{\left(k\right)}-H\left({\nu }_k\right){\left[F^{\prime }\left(y^{\left(k\right)}\right)\right]}^{-1}F\left(y^{\left(k\right)}\right)-G\left({\nu }_k\right){\left[F^{\prime }\left(y^{\left(k\right)}\right)\right]}^{-1}F\left(x^{\left(k\right)}\right),$$

and substituting all the previous expansions, we obtain

$$\begin{array}{cc}\hfill e^{(k+1)}=& -g_0{e}^{\left(k\right)}+\left(1-h_0-g_0\right)C_2{e}^{{\left(k\right)}^2}\hfill \\ & +\left[(1-h_0)\left(2C_3-2C_2^2\right)-g_0\left(C_3-2C_2^2\right)-g_1{C}_2^2\right]e^{{\left(k\right)}^3}\hfill \\ & +[4C_2^3-3C_3{C}_2-4C_2{C}_3+3C_4-h_0\left(3C_2^3-3C_3{C}_2-4C_2{C}_3+3C_4\right)-h_1{C}_2^3\hfill \\ & -g_0\left(2C_2^3-4C_2{C}_3+C_4\right)-g_1\left(-3C_2^3+2C_2{C}_3+2C_3{C}_2-P^{-1}Q{C}_2^2\right)]e^{{\left(k\right)}^4}\hfill \\ & +\Lambda e^{{\left(k\right)}^5}+\mathcal{O}\left(e^{{\left(k\right)}^6}\right),\hfill \end{array}$$

where $\Lambda$ denotes the coefficient of order $e^{{\left(k\right)}^5}$.

If we impose

$$h_0=1,h_1=1,g_0=0,g_1=0,$$

that is,

$$H\left(0\right)=1,H^{\prime }\left(0\right)=1,G\left(0\right)=0,G^{\prime }\left(0\right)=0,$$

then all the terms of orders $e^{\left(k\right)}$, $e^{{\left(k\right)}^2}$, $e^{{\left(k\right)}^3}$, and $e^{{\left(k\right)}^4}$ vanish. Under these conditions, the coefficient $\Lambda$ becomes

$$\begin{array}{cc}\hfill \Lambda =& \left(-8C_2^4+6C_3{C}_2^2+6C_2{C}_3{C}_2+8C_2^2{C}_3-6C_3^2-4C_4{C}_2-6C_2{C}_4+4C_5\right)\hfill \\ & -\left(-4C_2^4+6C_3{C}_2^2+4C_2{C}_3{C}_2+6C_2^2{C}_3-6C_3^2-4C_4{C}_2-6C_2{C}_4+4C_5\right)\hfill \\ & -\left(2C_2^2{C}_3+2C_2{C}_3{C}_2+2C_3{C}_2^2-6C_2^4-P^{-1}Q{C}_2^3\right)-g_2{C}_2^4.\hfill \end{array}$$

After simplification,

$$\Lambda =(2-g_2)C_2^4-2C_3{C}_2^2+P^{-1}Q{C}_2^3.$$

Since

$$g_2={\displaystyle \frac{1}{2}}{G}^{\prime \prime }\left(0\right),$$

we finally obtain

$$e^{(k+1)}=\left[\left(2-{\displaystyle \frac{1}{2}}{G}^{\prime \prime }\left(0\right)\right)C_2^4-2C_3{C}_2^2+P^{-1}Q{C}_2^3\right]e^{{\left(k\right)}^5}+\mathcal{O}\left(e^{{\left(k\right)}^6}\right).$$

Therefore, the iterative family has local order of convergence five.    □

Corollary 1.

Under the hypotheses of Theorem 1, assume that

$$G\left(\nu \right)\equiv 0.$$

Then the proposed weighted family reduces to

$$\left\{\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)}-H\left({\nu }_k\right){\left[F^{\prime }\left(y^{\left(k\right)}\right)\right]}^{-1}F\left(y^{\left(k\right)}\right),\hfill \end{array}\right.$$

(10)

where

$${\nu }_k={\displaystyle \frac{{∥F\left(y^{\left(k\right)}\right)∥}^2}{{∥F\left(x^{\left(k\right)}\right)∥}^2}}.$$

If the weight function H satisfies

$$H\left(0\right)=1,H^{\prime }\left(0\right)=1,$$

then the iterative scheme (10) has local order of convergence five.

In particular, for the choice

$$H\left(\nu \right)=1+\nu ,$$

the method (10) becomes

$$\left\{\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(1+{\nu }_k\right){\left[F^{\prime }\left(y^{\left(k\right)}\right)\right]}^{-1}F\left(y^{\left(k\right)}\right).\hfill \end{array}\right.$$

(11)

Therefore, the original fifth-order method of Singh et al. [9] is recovered as a particular member of the proposed family.

Proof. 

If $G\left(\nu \right)\equiv 0$, the general weighted scheme of Theorem 1 immediately reduces to (10). Moreover, the conditions required in Theorem 1 become

$$H\left(0\right)=1,H^{\prime }\left(0\right)=1,$$

which guarantee the cancellation of the lower-order terms in the error equation and, consequently, fifth-order convergence.

For $H\left(\nu \right)=1+\nu$, one has

$$H\left(0\right)=1,H^{\prime }\left(0\right)=1.$$

Substituting this weight function into (10) gives exactly (11), which coincides with the fifth-order method proposed by Singh et al. [9]. Hence, Singh’s method is a particular case of the present weighted family.    □

Remark 1.

Corollary 1 shows that the proposed family is not only inspired by Singh’s method, but actually contains it as a particular case. The introduction of the additional weight function G and the more general admissible choices of H provide a wider class of fifth-order schemes while preserving the same local convergence order.

Proposition 1.

Let $\alpha ,\beta ,\gamma ,\delta ,a,b\in \mathbb{R}$, with $a-b=1$ in case $\left(vii\right)$. Then, the following choices of weight functions satisfy

$$H\left(0\right)=1,H^{\prime }\left(0\right)=1,G\left(0\right)=0,G^{\prime }\left(0\right)=0,$$

and therefore generate fifth-order methods:

$$\begin{array}{cccc}\hfill \left(\mathrm{i}\right)& H\left(\nu \right)=1+\nu ,\hfill & \hfill & G\left(\nu \right)=0,\hfill \\ \hfill \left(\mathrm{ii}\right)& H\left(\nu \right)=1+\nu +\alpha {\nu }^2,\hfill & \hfill & G\left(\nu \right)=0,\hfill \\ \hfill \left(\mathrm{iii}\right)& H\left(\nu \right)=1+\nu ,\hfill & \hfill & G\left(\nu \right)=\gamma {\nu }^2,\hfill \\ \hfill \left(\mathrm{iv}\right)& H\left(\nu \right)=1+\nu +\alpha {\nu }^2,\hfill & \hfill & G\left(\nu \right)=\beta {\nu }^2,\hfill \\ \hfill \left(\mathrm{v}\right)& H\left(\nu \right)=1+\nu +a{\nu }^3,\hfill & \hfill & G\left(\nu \right)=0,\hfill \\ \hfill \left(\mathrm{vi}\right)& H\left(\nu \right)={\displaystyle \frac{1}{1-\nu }},\hfill & \hfill & G\left(\nu \right)=\gamma {\nu }^2,\hfill \\ \hfill \left(\mathrm{vii}\right)& H\left(\nu \right)={\displaystyle \frac{1+a\nu }{1+b\nu }},\hfill & \hfill & G\left(\nu \right)=\delta {\nu }^2,a-b=1.\hfill \end{array}$$

(12)

4. Numerical Experiments

For comparison purposes, the proposed weighted family is tested against several reference schemes from the literature, including methods of comparable convergence order and the classical Newton method. More precisely, the numerical experiments include the modified Newton–Jarratt composition of Cordero et al. [13], the fifth-order method denoted in this work by Vassileva and taken from Arroyo et al. [14], the fifth-order method of Sharma et al., labeled SHM5 [15], and the method denoted by O6, corresponding in our tests to the fifth-order member of the arbitrary odd-order family introduced by Solaiman and Hashim [16]. These benchmark procedures are employed, together with Newton’s method, to compare the accuracy, robustness, and computational performance of the proposed schemes on the selected nonlinear systems.

To facilitate the interpretation and reproducibility of the numerical experiments, Table 1 shows the correspondence between the labels NMA1–NMA5 and the particular choices of the weight functions selected from Proposition 1.

The following experiments are intended to assess practical behavior and do not replace the local convergence analysis established in Section 3.

Newton–Jarratt.

$$\begin{array}{cc}& \left\{\begin{array}{cc}\hfill z^{\left(k\right)}& =x^{\left(k\right)}-{\displaystyle \frac{2}{3}}{\left[F^{\prime }\left(x^{\left(k\right)}\right)\right]}^{-1}F\left(x^{\left(k\right)}\right),\hfill \\ \hfill y^{\left(k\right)}& =x^{\left(k\right)}-{\displaystyle \frac{1}{2}}{\left[3F^{\prime }\left(z^{\left(k\right)}\right)-F^{\prime }\left(x^{\left(k\right)}\right)\right]}^{-1}\left[3F^{\prime }\left(z^{\left(k\right)}\right)+F^{\prime }\left(x^{\left(k\right)}\right)\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[\alpha F^{\prime }\left(x^{\left(k\right)}\right)+\beta F^{\prime }\left(z^{\left(k\right)}\right)\right]}^{-1}F\left(y^{\left(k\right)}\right),\alpha +\beta =1.\hfill \end{array}\right.\hfill \end{array}$$

Vassileva.

$$\begin{array}{cc}& \left\{\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[F^{\prime }\left(x^{\left(k\right)}\right)-5F^{\prime }\left(y^{\left(k\right)}\right)\right]}^{-1}\left[3F^{\prime }\left(x^{\left(k\right)}\right)+F^{\prime }\left(y^{\left(k\right)}\right)\right]{\left[F^{\prime }\left(x^{\left(k\right)}\right)\right]}^{-1}F\left(y^{\left(k\right)}\right).\hfill \end{array}\right.\hfill \end{array}$$

Newton.

$$\begin{array}{cc}& \left\{\begin{array}{cc}\hfill 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).\hfill \end{array}\right.\hfill \end{array}$$

SHM5.

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

O6.

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

4.1. Computational efficiency analysis

In order to complement the local convergence analysis established in Section 3, we now study the computational efficiency of the proposed weighted family. In the vectorial setting, the classical efficiency index

$$E=p^{1/m},$$

where p is the convergence order and m is the number of evaluations per iteration, is not sufficiently discriminating, since it does not distinguish between evaluations of the nonlinear operator, evaluations of the Jacobian matrix, and the cost of solving the associated linear systems. This issue is well known in the literature on iterative methods for nonlinear systems; see, for instance, [17,18]. Therefore, we also consider the computational efficiency index

$$\mathrm{CI}\left(n\right)=p^{1/C\left(n\right)},$$

where $C\left(n\right)$ denotes the total cost per iteration for a nonlinear system of dimension n.

4.1.1. Cost Model

We adopt the following cost model:

$$\mathrm{one}\mathrm{evaluation}\mathrm{of}F\sim n,\mathrm{one}\mathrm{evaluation}\mathrm{of}{F}^{\prime }\sim n^2,$$

$$\mathrm{one}\mathrm{LU}\mathrm{factorization}\sim {\displaystyle \frac{n^3}{3}}-{\displaystyle \frac{n}{3}},\mathrm{one}\mathrm{solve}\mathrm{with}\mathrm{one}\mathrm{vector}\mathrm{right}-\mathrm{hand}\mathrm{side}\sim n^2,$$

$$\mathrm{one}\mathrm{matrix}-\mathrm{vector}\mathrm{product}\sim n^2,\mathrm{one}\mathrm{linear}\mathrm{combination}\mathrm{of}\mathrm{two}n\times n\mathrm{matrices}\sim n^2.$$

Products and quotients are counted, whereas vector additions and subtractions are not included. In the following notation, $n_F$ denotes the number of evaluations of the nonlinear operator F, $n_J$ the number of Jacobian evaluations, $n_{\mathrm{LU}}$ the number of LU factorizations, $n_V$ the number of linear systems solved with a vector right-hand side, $n_{MV}$ the number of matrix–vector products, and $n_A$ the number of linear combinations of matrices. Therefore, for a method with these per-iteration counts, the functional and algebraic costs are given by

$$d\left(n\right)=n_{F}n+n_J{n}^2$$

and

$$op\left(n\right)=n_{\mathrm{LU}}\left({\displaystyle \frac{n^3}{3}}-{\displaystyle \frac{n}{3}}\right)+n_V{n}^2+(n_{MV}+n_A)n^2+\mathrm{additional}\mathrm{scalar}/\mathrm{vector}\mathrm{products}.$$

Thus,

$$C\left(n\right)=d\left(n\right)+op\left(n\right),\mathrm{CI}\left(n\right)=p^{1/C\left(n\right)}.$$

4.1.2. Cost of the Weighted Family

The proposed family is

$$\left\{\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[F^{\prime }\left(y^{\left(k\right)}\right)\right]}^{-1}\left(H\left({\nu }_k\right)F\left(y^{\left(k\right)}\right)+G\left({\nu }_k\right)F\left(x^{\left(k\right)}\right)\right),\hfill \end{array}\right.$$

where

$${\nu }_k={\displaystyle \frac{{∥F\left(y^{\left(k\right)}\right)∥}^2}{{∥F\left(x^{\left(k\right)}\right)∥}^2}}.$$

Each iteration requires two evaluations of F, two evaluations of $F^{\prime }$, two LU factorizations, and two vector solves. Therefore,

$$d_{\mathrm{NMA}}\left(n\right)=2n+2n^2.$$

The two LU factorizations and the two vector solves contribute

$$2\left({\displaystyle \frac{n^3}{3}}-{\displaystyle \frac{n}{3}}\right)+2n^2={\displaystyle \frac{2}{3}}{n}^3+2n^2-{\displaystyle \frac{2}{3}}n.$$

Moreover, the computation of ${\nu }_k$ requires two inner products and one scalar quotient, contributing $2n+1$. The formation of

$$H\left({\nu }_k\right)F\left(y^{\left(k\right)}\right)+G\left({\nu }_k\right)F\left(x^{\left(k\right)}\right)$$

requires two scalar-vector products, contributing $2n$. Hence,

$${op}_{\mathrm{NMA}}\left(n\right)={\displaystyle \frac{2}{3}}{n}^3+2n^2-{\displaystyle \frac{2}{3}}n+(2n+1)+2n={\displaystyle \frac{2}{3}}{n}^3+2n^2+{\displaystyle \frac{10}{3}}n+1.$$

Consequently,

$$C_{\mathrm{NMA}}\left(n\right)={\displaystyle \frac{2}{3}}{n}^3+4n^2+{\displaystyle \frac{16}{3}}n+1,$$

and

$${\mathrm{CI}}_{{\mathrm{NMA}}_i}\left(n\right)=5^{{\displaystyle \frac{1}{\frac{2}{3}{n}^3+4n^2+\frac{16}{3}n+1}}},i=1,\dots ,5.$$

Since the weight functions only modify scalar coefficients, all the members NMA1–NMA5 have the same computational efficiency index.

4.1.3. Per-Iteration Counts

The operation counts of the proposed family and the comparison methods are summarized in Table 2.

4.1.4. Computational efficiency indices

Using the counts in Table 2, we obtain the following expressions.

For the proposed weighted family,

$${\mathrm{CI}}_{{\mathrm{NMA}}_i}\left(n\right)=5^{{\displaystyle \frac{1}{\frac{2}{3}{n}^3+4n^2+\frac{16}{3}n+1}}},i=1,\dots ,5.$$

For Newton–Jarratt,

$$d_{\mathrm{NJ}}\left(n\right)=2n+2n^2,$$

$${op}_{\mathrm{NJ}}\left(n\right)=3\left({\displaystyle \frac{n^3}{3}}-{\displaystyle \frac{n}{3}}\right)+3n^2+4n^2=n^3+7n^2-n,$$

and therefore

$${\mathrm{CI}}_{\mathrm{NJ}}\left(n\right)=5^{{\displaystyle \frac{1}{n^3+9n^2+n}}}.$$

For Vassileva,

$$d_{\mathrm{Vass}}\left(n\right)=2n+2n^2,$$

$${op}_{\mathrm{Vass}}\left(n\right)=2\left({\displaystyle \frac{n^3}{3}}-{\displaystyle \frac{n}{3}}\right)+3n^2+3n^2={\displaystyle \frac{2}{3}}{n}^3+6n^2-{\displaystyle \frac{2}{3}}n,$$

and hence

$${\mathrm{CI}}_{\mathrm{Vass}}\left(n\right)=5^{{\displaystyle \frac{1}{\frac{2}{3}{n}^3+8n^2+\frac{4}{3}n}}}.$$

For Newton,

$$d_{\mathrm{Newton}}\left(n\right)=n+n^2,$$

$${op}_{\mathrm{Newton}}\left(n\right)=\left({\displaystyle \frac{n^3}{3}}-{\displaystyle \frac{n}{3}}\right)+n^2={\displaystyle \frac{n^3}{3}}+n^2-{\displaystyle \frac{n}{3}},$$

and thus

$${\mathrm{CI}}_{\mathrm{Newton}}\left(n\right)=2^{{\displaystyle \frac{1}{\frac{n^3}{3}+2n^2+\frac{2}{3}n}}}.$$

For SHM5,

$$d_{\mathrm{SHM}5}\left(n\right)=2n+2n^2,$$

$${op}_{\mathrm{SHM}5}\left(n\right)=2\left({\displaystyle \frac{n^3}{3}}-{\displaystyle \frac{n}{3}}\right)+4n^2={\displaystyle \frac{2}{3}}{n}^3+4n^2-{\displaystyle \frac{2}{3}}n,$$

and consequently

$${\mathrm{CI}}_{\mathrm{SHM}5}\left(n\right)=5^{{\displaystyle \frac{1}{\frac{2}{3}{n}^3+6n^2+\frac{4}{3}n}}}.$$

For O6,

$$d_{\mathrm{O}6}\left(n\right)=2n+2n^2,$$

$${op}_{\mathrm{O}6}\left(n\right)=3\left({\displaystyle \frac{n^3}{3}}-{\displaystyle \frac{n}{3}}\right)+3n^2+n^2=n^3+4n^2-n,$$

and therefore

$${\mathrm{CI}}_{\mathrm{O}6}\left(n\right)=5^{{\displaystyle \frac{1}{n^3+6n^2+n}}}.$$

These expressions are summarized in Table 3.

From Table 3, it follows that all weighted variants NMA1–NMA5 have the same computational efficiency index. This is expected, since the different weight functions do not change the number of evaluations, factorizations, or linear systems solved. In addition, within the adopted cost model, the proposed family has the smallest denominator among the fifth-order methods considered and therefore exhibits the largest computational efficiency index in the comparison.

For a clearer visualization of the results in Table 3, Figure 1 displays the behavior of $\mathrm{CI}\left(n\right)$ for the proposed weighted family and the comparison methods in both small- and large-dimensional regimes.

The graphs confirm that, although the computational efficiency index decreases with the problem dimension for all schemes, the proposed weighted family maintains one of the most favorable behaviors in both the small- and large-dimensional regimes.

4.1.5. Solution of Some Academic Problems

This section presents a numerical assessment of the iterative schemes within a common computational framework. The aim is to compare their practical behavior in terms of accuracy, efficiency, and convergence. For each method, the tables report the total number of iterations required to satisfy the stopping criterion, the quantity $∥x^{(k+1)}-x^{\left(k\right)}∥$, which measures the distance between the last two iterates, the residual norm $∥F\left(x^{(k+1)}\right)∥$ at the final approximation, the execution time in seconds (e-time), together with the approximated computational order of convergence (ACOC), see [12].

To ensure a fair and consistent comparison, all numerical tests were carried out under the same computational conditions. In every experiment, the stopping criterion was defined by

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

with a maximum of 50 iterations allowed. All computations were performed using variable precision arithmetic through vpa with 1500 significant digits. Such a high-precision setting reduces the influence of round-off errors and makes it possible to observe the asymptotic convergence regime more clearly.

In order to obtain representative timing results, each iterative scheme was executed 10 times for every test problem, and the reported e-time values correspond to the average over these runs. This averaging procedure reduces the influence of incidental fluctuations caused by background processes and operating-system activity. In addition, all methods were applied to each system using the same initial approximation, which guarantees methodological uniformity throughout the comparison.

The experiments were carried out on a computer running macOS Tahoe (version 26.4), equipped with an Apple M2 processor and 8 GB of RAM (2023 model). Therefore, the reported timings reflect the methods’ performance on a modern mid-range computational platform.

The nonlinear systems used in the numerical experiments are described below. For each test problem, we present its analytical formulation, which defines the framework for the comparative study developed in this section.

Example 1.

We begin with the nonlinear system of n equations

$$x_i-1.5sin\left({\displaystyle \sum _{\begin{array}{c}j=1\\ j\ne i\end{array}}^n}{x}_j\right)=0,i=1,2,\dots ,n,$$

with $n=40$, whose solution is

$$\xi \approx {\left(0.2375825578,\dots ,0.2375825578\right)}^T.$$

The initial approximation was taken as

$$x^{\left(0\right)}={\left({\displaystyle \frac{1}{4}},\dots ,{\displaystyle \frac{1}{2}}\right)}^T.$$

Example 2.

We next consider the nonlinear system

$$x_i-cos\left(-{\displaystyle \sum _{j=1}^n}{x}_j+2x_i\right)=0,i=1,2,\dots ,n,$$

with $n=30$, whose solution is

$$\xi \approx {\left(0.4867431909,\dots ,0.4867431909\right)}^T.$$

As initial guess, we used

$$x^{\left(0\right)}={\left({\displaystyle \frac{1}{2}},\dots ,{\displaystyle \frac{1}{2}}\right)}^T.$$

Example 3.

Consider now the nonlinear system

$$x_i+5-2log\left(1+{\displaystyle \sum _{\begin{array}{c}j=1\\ j\ne i\end{array}}^n}{x}_j\right)=0,i=1,2,\dots ,n,$$

with $n=25$, whose solution is

$$\xi \approx {\left(4.2863062097,\dots ,4.2863062097\right)}^T.$$

The seed was chosen as

$$x^{\left(0\right)}={\left(4,\dots ,4\right)}^T.$$

Example 4.

We now analyze the nonlinear system

$$\left\{\begin{array}{c}{x}_1+log(2+x_1+x_2)=0,\hfill \\ x_2+log(2+x_2+x_3)=0,\hfill \\ ⋮\hfill \\ x_{998}+log(2+x_{998}+x_{999})=0,\hfill \\ x_{999}+log(2+x_{999}+x_1)=0,\hfill \end{array}\right.$$

where $n=999$. Its solution is approximately

$$\xi \approx {(-0.3149230578,\dots ,-0.3149230578)}^T.$$

The initial approximation was taken as

$$x^{\left(0\right)}={\left({\displaystyle \frac{1}{10}},\dots ,{\displaystyle \frac{1}{10}}\right)}^T.$$

Example 5.

We now consider the nonlinear system

$$\left\{\begin{array}{c}{x}_{1}sin\left(x_2\right)-1=0,\hfill \\ x_{2}sin\left(x_3\right)-1=0,\hfill \\ ⋮\hfill \\ x_{499}sin\left(x_{500}\right)-1=0,\hfill \\ x_{500}sin\left(x_1\right)-1=0,\hfill \end{array}\right.$$

with $n=500$. The corresponding solution is

$$\xi \approx {(1.1141571409,\dots ,1.1141571409)}^T.$$

The initial approximation was taken as

$$x^{\left(0\right)}={\left({\displaystyle \frac{1}{2}},\dots ,{\displaystyle \frac{1}{2}}\right)}^T.$$

Example 6.

We also study the nonlinear system

$$\left\{\begin{array}{c}{x}_1{x}_2-1=0,\hfill \\ x_2{x}_3-1=0,\hfill \\ ⋮\hfill \\ x_{1498}{x}_{1499}-1=0,\hfill \\ x_{1499}{x}_1-1=0,\hfill \end{array}\right.$$

with $n=1499$. In this case, the solution is

$$\xi \approx {\left(1.0000000000,\dots ,1.0000000000\right)}^T.$$

The initial approximation was chosen as

$$x^{\left(0\right)}={\left({\displaystyle \frac{3}{2}},\dots ,{\displaystyle \frac{3}{2}}\right)}^T.$$

Example 7.

Finally, we consider the nonlinear system

$$\left\{\begin{array}{c}{x}_1{x}_2-e^{-x_1}-e^{-x_2}=0,\hfill \\ x_2{x}_3-e^{-x_2}-e^{-x_3}=0,\hfill \\ ⋮\hfill \\ x_{1498}{x}_{1499}-e^{-x_{1498}}-e^{-x_{1499}}=0,\hfill \\ x_{1499}{x}_1-e^{-x_{1499}}-e^{-x_1}=0,\hfill \end{array}\right.$$

with $n=1499$. Its solution is approximately

$$\xi \approx {(0.9012010317,\dots ,0.9012010317)}^T.$$

The initial estimate was taken as

$$x^{\left(0\right)}={\left({\displaystyle \frac{6}{5}},\dots ,{\displaystyle \frac{6}{5}}\right)}^T.$$

In the seven test problems listed in Table 4, Table 5, Table 6 and Table 7 and Table 8, Table 9 and Table 10, the proposed family $\mathrm{NMA}1,\dots ,\mathrm{NMA}5$ demonstrates consistently competitive performance. In all cases, the NMA variants converge in only three or four iterations, whereas Newton requires seven or eight iterations. Moreover, among the high-order methods, at least one NMA variant attains the lowest execution time in every example. When Newton is also included in the comparison, an NMA variant remains the fastest method in six of the seven tests. The only exception is Example (5), where Newton gives the lowest runtime; even in this case, the NMA variants are clearly faster than Newton–Jarratt, Vassileva, SHM5, and O6. In addition, the observed values of p remain close to the theoretical order five in most cases, confirming that the proposed weighted family provides a robust and efficient alternative for solving large-scale nonlinear systems.

5. Application to a Stationary Viscous Burgers Problem

In order to illustrate the applicability of the proposed iterative family to nonlinear differential models, we consider a stationary one-dimensional viscous Burgers problem. This equation is a classical nonlinear model describing the interaction between convection and diffusion, and it is frequently used as a prototype in fluid mechanics, transport phenomena, and nonlinear wave propagation [19].

We study the boundary value problem

$$-\nu u^{\prime \prime }\left(x\right)+u\left(x\right)u^{\prime }\left(x\right)=f\left(x\right),x\in (0,1),$$

(13)

subject to the homogeneous Dirichlet boundary conditions

$$u\left(0\right)=0,u\left(1\right)=0,$$

(14)

where $\nu >0$ denotes the viscosity parameter.

In order to validate the numerical solution, we use the method of manufactured solutions [20,21]. More precisely, we prescribe the exact solution

$$u\left(x\right)=x(1-x),$$

(15)

which satisfies the boundary conditions (14). Its derivatives are

$$u^{\prime }\left(x\right)=1-2x,u^{\prime \prime }\left(x\right)=-2.$$

Substituting (15) into (13), the forcing term is obtained as

$$f\left(x\right)=-\nu u^{\prime \prime }\left(x\right)+u\left(x\right)u^{\prime }\left(x\right)=2\nu +x(1-x)(1-2x).$$

(16)

Therefore, the model problem under consideration is

$$-\nu u^{\prime \prime }\left(x\right)+u\left(x\right)u^{\prime }\left(x\right)=2\nu +x(1-x)(1-2x),x\in (0,1),$$

(17)

with

$$u\left(0\right)=0,u\left(1\right)=0.$$

The choice of this model is convenient for two reasons. First, it preserves the nonlinear convection–diffusion structure characteristic of Burgers-type equations [19]. Second, the exact solution is explicitly known, which allows us to verify the quality of the discrete approximation and to assess the performance of the iterative methods on the nonlinear algebraic system arising from the discretization.

Finite difference discretization

Let N be the number of interior nodes and define the uniform mesh

$$x_i=ih,i=0,1,\dots ,N+1,$$

with mesh size

$$h={\displaystyle \frac{1}{N+1}}.$$

At the interior nodes $x_i$, $i=1,\dots ,N$, we denote by

$$u_i\approx u\left(x_i\right)$$

the numerical approximation of the exact solution.

To discretize the derivatives, we employ the standard centered finite difference formulas [22]. For the second derivative, we use

$$u^{\prime \prime }\left(x_i\right)\approx {\displaystyle \frac{u_{i-1}-2u_i+u_{i+1}}{h^2}},$$

(18)

while for the first derivative we take

$$u^{\prime }\left(x_i\right)\approx {\displaystyle \frac{u_{i+1}-u_{i-1}}{2h}}.$$

(19)

Substituting (18) and (19) into (17), we obtain, for each interior node $x_i$,

$$-\nu {\displaystyle \frac{u_{i-1}-2u_i+u_{i+1}}{h^2}}+u_i{\displaystyle \frac{u_{i+1}-u_{i-1}}{2h}}-f\left(x_i\right)=0,i=1,\dots ,N.$$

Using the boundary conditions,

$$u_0=0,u_{N+1}=0,$$

the discrete problem becomes a nonlinear system of N equations with N unknowns.

Nonlinear algebraic system

Let $U={(u_1,u_2,\dots ,u_N)}^T\in {\mathbb{R}}^N$. The nonlinear system can be written in compact form as

$$F\left(U\right)=0,$$

(20)

where $F\left(U\right)={(F_1\left(U\right),F_2\left(U\right),\dots ,F_N\left(U\right))}^T$, with components given by

$$F_i\left(U\right)=-\nu {\displaystyle \frac{u_{i-1}-2u_i+u_{i+1}}{h^2}}+u_i{\displaystyle \frac{u_{i+1}-u_{i-1}}{2h}}-f\left(x_i\right),i=1,\dots ,N,$$

(21)

together with the conventions

$$u_0=0,u_{N+1}=0.$$

More explicitly, the first and last equations read

$$F_1\left(U\right)=-\nu {\displaystyle \frac{-2u_1+u_2}{h^2}}+u_1{\displaystyle \frac{u_2}{2h}}-f\left(x_1\right),$$

$$F_N\left(U\right)=-\nu {\displaystyle \frac{u_{N-1}-2u_N}{h^2}}-u_N{\displaystyle \frac{u_{N-1}}{2h}}-f\left(x_N\right),$$

while for $i=2,\dots ,N-1$,

$$F_i\left(U\right)=-\nu {\displaystyle \frac{u_{i-1}-2u_i+u_{i+1}}{h^2}}+u_i{\displaystyle \frac{u_{i+1}-u_{i-1}}{2h}}-f\left(x_i\right).$$

This system is nonlinear because of the convective term

$$u_i{\displaystyle \frac{u_{i+1}-u_{i-1}}{2h}}.$$

Consequently, an iterative method is required in order to compute the discrete solution.

Jacobian matrix

Since the methods considered in this work require the Jacobian matrix, we now derive its explicit form. Differentiating (21) with respect to the neighboring unknowns, we obtain

$${\displaystyle \frac{\partial F_i}{\partial u_{i-1}}}=-{\displaystyle \frac{\nu }{h^2}}-{\displaystyle \frac{u_i}{2h}},$$

$${\displaystyle \frac{\partial F_i}{\partial u_i}}={\displaystyle \frac{2\nu }{h^2}}+{\displaystyle \frac{u_{i+1}-u_{i-1}}{2h}},$$

$${\displaystyle \frac{\partial F_i}{\partial u_{i+1}}}=-{\displaystyle \frac{\nu }{h^2}}+{\displaystyle \frac{u_i}{2h}}.$$

All other partial derivatives are zero. Therefore, the Jacobian matrix

$$J\left(U\right)=F^{\prime }\left(U\right)$$

is tridiagonal and can be written as

$$J\left(U\right)=\left(\begin{array}{ccccc}{\displaystyle \frac{2\nu }{h^2}}+{\displaystyle \frac{u_2}{2h}}& -{\displaystyle \frac{\nu }{h^2}}+{\displaystyle \frac{u_1}{2h}}& 0& \dots & 0\\ -{\displaystyle \frac{\nu }{h^2}}-{\displaystyle \frac{u_2}{2h}}& {\displaystyle \frac{2\nu }{h^2}}+{\displaystyle \frac{u_3-u_1}{2h}}& -{\displaystyle \frac{\nu }{h^2}}+{\displaystyle \frac{u_2}{2h}}& \ddots & ⋮\\ 0& -{\displaystyle \frac{\nu }{h^2}}-{\displaystyle \frac{u_3}{2h}}& {\displaystyle \frac{2\nu }{h^2}}+{\displaystyle \frac{u_4-u_2}{2h}}& \ddots & 0\\ ⋮& \ddots & \ddots & \ddots & -{\displaystyle \frac{\nu }{h^2}}+{\displaystyle \frac{u_{N-1}}{2h}}\\ 0& \dots & 0& -{\displaystyle \frac{\nu }{h^2}}-{\displaystyle \frac{u_N}{2h}}& {\displaystyle \frac{2\nu }{h^2}}-{\displaystyle \frac{u_{N-1}}{2h}}\end{array}\right).$$

(22)

Hence, the discrete problem (20) leads to a sparse nonlinear algebraic system whose Jacobian has a structured tridiagonal form. This makes it especially suitable for assessing the practical performance of the iterative methods studied in this work.

Numerical solution procedure

To solve the nonlinear system (20), one may choose as initial approximation, for instance,

$$U^{\left(0\right)}={(0,0,\dots ,0)}^T,$$

or, if desired, the exact solution sampled at the grid points plus a small perturbation. Once the mesh size h, the viscosity parameter $\nu$, and the forcing term $f\left(x_i\right)$ are fixed, the iterative schemes are applied directly to the nonlinear system (20).

For each computed approximation $U^{\left(k\right)}$, the quality of the numerical solution may be assessed through the residual norm

$$\parallel F\left(U^{\left(k\right)}\right)\parallel ,$$

and, since the exact solution is known, also through the discrete error

$$E^{\left(k\right)}=\underset{1\le i\le N}{max}|u_i-u\left(x_i\right)|.$$

This allows us to compare not only the convergence speed of the iterative methods, but also the accuracy of the final discrete approximation.

The previous construction provides a complete test framework: a nonlinear differential model of physical interest, a manufactured exact solution, a finite difference discretization, and an explicit nonlinear algebraic system. Therefore, the stationary Burgers problem constitutes a suitable benchmark for validating the theoretical and practical performance of the proposed high-order iterative family [19].

For the numerical solution of the nonlinear algebraic system arising from the finite difference discretization of the stationary Burgers problem, we considered $N=300$ interior nodes, viscosity parameter $\nu =0.09$, and mesh size $h=0.00332226$. The stopping criterion was fixed as

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

with a maximum of 50 iterations. Table 11 reports, for each method, the total number of iterations, the increment norm $\parallel x^{(k+1)}-x^{\left(k\right)}\parallel$, the residual norm $\parallel F\left(x^{(k+1)}\right)\parallel$, the execution time in seconds, the estimated computational order of convergence p, and the infinity norm of the error with respect to the exact solution.

From Table 11, it can be observed that the proposed variants $\mathrm{NMA}1,\dots ,\mathrm{NMA}5$ converge in four iterations and attain very small residual and error norms. They are faster than Newton–Jarratt, Vassileva, SHM5, and O6, while Newton is the only method with a lower execution time, although it requires seven iterations and remains quadratically convergent. The ACOC values of the NMA variants are close to four rather than five; however, this empirical behavior is not unusual in nonlinear systems obtained from discretized differential equations. In such problems, the observed order may be influenced by mesh effects, conditioning, stopping criteria, and numerical saturation before the asymptotic regime is fully reflected by the ACOC estimator. Thus, the table still confirms the efficiency and accuracy of the proposed weighted family for the stationary Burgers problem.

5.1. Graphical Analysis for the Stationary Burgers Problem

To complement the numerical results reported above, this subsection presents a graphical analysis of the performance of NMA2 when applied to the stationary Burgers problem. The visual comparison includes the agreement between the exact and numerical solutions, the distribution of the absolute error, the convergence history of the selected methods, and the evolution of the iterates generated during the nonlinear resolution process.

Figure 2a and Figure 2b show that the numerical solution obtained with NMA2 is practically indistinguishable from the exact solu tion and that the absolute error remains negligible over the computational mesh. Figure 3a confirms the fast decrease of the residual norm for NMA2 compared with the reference methods, while Figure 3b illustrates the stable evolution of the iterates toward the final discrete solution. Overall, these plots provide graphical evidence of the accuracy, stability, and fast convergence of NMA2 for this nonlinear differential problem.

6. Conclusions

In this work, we introduced and analyzed a weighted family of fifth-order iterative schemes for solving nonlinear systems of equations. The proposed methods, denoted by $\mathrm{NMA}1,\dots ,\mathrm{NMA}5$, were constructed by incorporating suitable weight functions into a two-step framework, leading to a class of schemes that preserve the expected local order of convergence while maintaining a competitive computational cost.

The theoretical analysis established the conditions that the weight functions must satisfy in order to guarantee fifth-order convergence. Several admissible choices were proposed, giving rise to different variants of the method. This flexibility makes the family adaptable from both the analytical and computational points of view, since the particular form of the weight functions can influence the practical efficiency of the resulting scheme.

The numerical experiments performed on a broad collection of nonlinear systems confirmed the reliability and effectiveness of the proposed methods. In all test problems, the compared schemes converged to the same solution, so the reported differences reflect their numerical performance. The proposed variants were consistently competitive and often outperformed classical reference methods in execution time, while maintaining a convergence order close to the theoretical one. Among them, $\mathrm{NMA}2$ showed a particularly balanced behavior, combining fast convergence, small residuals, and numerical stability.

In addition, the application to the stationary Burgers problem illustrated the suitability of the proposed family for nonlinear systems arising from the discretization of differential models. The graphical analysis showed an excellent agreement between the numerical and exact solutions, a uniformly small absolute error, and a rapid decay of the residual through the iterations. The mesh representation of the iterates also provided a clear visual confirmation of the stability and fast convergence of the method.

Overall, the proposed weighted family constitutes a robust, accurate, and efficient alternative for solving nonlinear systems. The theoretical results, together with the numerical evidence, indicate that these schemes are promising candidates for the treatment of large-scale nonlinear problems and for applications involving discretized differential equations. Future work may include the extension of this family to Jacobian-free formulations, the incorporation of memory techniques, and its application to more general nonlinear partial differential equations.