A Class of Efficient Sixth-Order Iterative Methods for Solving the Nonlinear Shear Model of a Reinforced Concrete Beam

Abstract

In this paper, we present a three-step sixth-order class of iterative schemes to estimate the solutions of a nonlinear system of equations. This procedure is designed by means of a weight function technique. We apply this procedure for predicting the shear strength of a reinforced concrete beam. The values for the parameters of the nonlinear system describing this problem were randomly selected inside the prescribed ranges by technical standards for structural concrete. Moreover, some of these parameters were fixed taking into consideration the solvability region of the adopted steel constitutive model. The effectiveness of the new class is also compared with other current schemes in terms of the computational efficiency and numerical performance, with very good results. The advantages of this new class come from the low computational cost, due to the existence of an only inverse operator.

1. Introduction

Reinforced and prestressed concrete beams represent a structural type that resists internal stresses in a relatively complex manner due to their constitutive nature. Prior to the cracking of the concrete, the shear loads are carried by a set of diagonal compressive stresses complemented by another set of diagonal tensile stresses acting perpendicular to the first ones. Once the concrete tensile strength is reached, cracks form in the direction normal to the diagonal tensile stresses while pre-existing cracks spread and change inclination. Then, the ability of concrete to transmit diagonal tensile stresses is significantly reduced and the appropriate reinforcement is necessary to create a new system of internal stresses that carry the shear acting on the beam after cracking.

Between 1899 and 1902, Ritter [1] and Mörsch [2] proposed a truss model for explaining the field of forces in a cracked reinforced concrete beam, with the principal compressive stresses acting as diagonal members at 45º and the stirrups acting as vertical tension members. This model neglected the tensile stresses in the cracked concrete. In 1910, the first ACI Code modified Mörsch’s 45º truss model through the addition of concrete in order to compensate for the conservatism of the model and to account for the fact that the crack angle is usually less than 45º [3]. Between 1904 and 1922, Talbot and Withey demonstrated that the stirrup stresses were lower than those predicted by the 45º truss model [4].

Thus, before using the equilibrium equations, the inclination of the diagonal compressive struts should be known. In 1929, Wagner [5] treated a similar problem by studying the post-buckling shear response of thin metal girders, and he assumed that the angle of inclination of the diagonal tensile stresses coincided with the corresponding value of the diagonal tensile strains. Between 1974 and 1978, Collins and Mitchell developed a shear design model for reinforced concrete that, based on Wagner’s assumption, predicted the inclination of the compressive struts by considering the strains in the transverse and longitudinal reinforcement, and in the diagonally stressed concrete [6,7]; this last approach became known as the Compression Field Theory (CFT).

In 1982, Vecchio and Collins found that the principal compressive stress in the concrete is a function not only of the principal compressive strain, but also of the coexisting principal tensile strain [8]. In 1986, they published, as a further development of the CFT, the so-called Modified Compression Field Theory (MCFT) [9], which accounts for the influence of the tensile stresses in the cracked concrete. This theory assumes a bilinear constitutive model for the steel and requires the checking of the local bar stresses at the cracks, ensuring thus that the smeared steel stresses between adjacent cracks are lower than the yield value.

Since the MCFT, other alternative approaches have been also developed due to the consideration given to other strength mechanisms, such as the dowel action of the reinforcement intersecting the cracks or the friction between the crack faces [10]. Likewise, different procedures have been proposed to treat the shear response of a reinforced concrete member in a continuum mechanics context (i.e., to account for tensile stresses in the diagonally cracked concrete), such as the developments of Hsu and his co-workers at the University of Houston [11,12], in the framework of the so-called Rotating-Angle Softened Truss Model (RA-STM). The RA-STM proposes a constitutive relationship for the reinforcement as stiffened by the concrete (i.e., the embedded model bar); due to this alternative concept, the steel stress does not exceed the yielding point and the local checking of cracks is no longer necessary. The last contribution in this line is the so-called Refined Compression Field Theory (RCFT) [13,14], where the embedded bar stress–strain relationship is obtained from the concrete tension stiffening model considered in the MCFT, so the crack check is avoided and a new formulation with respect to the traditional MCFT one is no longer needed. Moreover, the numerical results obtained from the RCFT lead to a better fitting of the experimental results, particularly in the region near the peak point in the shear–strain response, where MCFT significantly deviates from the experimental data. This last approach to the constitutive modeling of the steel reinforcement is the one considered in this work.

As it is justified in Section 2, CFT mechanical models involve several types of nonlinearities, among other reasons, due to the constitutive relationships of the reinforced concrete. These models usually require in their implementation the application of iterative fixed-point methods to solve the nonlinear systems that appear. In fact, the most efficient way to solve a nonlinear problem is usually to choose between accuracy and computational cost [15]. Moreover, in this work, the previous determination of a solvability region using algebraic procedures is also necessary in order to improve the efficiency of the numerical solver, as indicated in Section 2.

Solving systems of nonlinear equations is an important problem in science and engineering, as has been previously described. The objective is to find the roots of the nonlinear system $F\left(x\right)=0$, F being a multidimensional function, $F:D\subseteq {\mathbb{R}}^n\to {\mathbb{R}}^n$, on D convex set, of size $n\times n$, $F\left(x\right)={(f_1\left(x\right),f_2\left(x\right),\dots ,f_n\left(x\right))}^T$ being $f_i$, $i=1,2,\dots ,n$, the functional coordinates of F.

One of the most commonly used methods is the classical Newton’s method, which has a quadratic order of convergence and iterative expression:

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

(1)

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

Several Newton-type procedures, by using different techniques, have been published in the last few years, some introductory texts to this area are [16,17,18,19,20]. Their main aim is accelerating their convergence or increasing their efficiency with differently designed techniques used by numerous authors in the literature (see, for example, [21,22,23,24,25,26]). In what follows, we are going to recall some of them for comparison purposes.

All the schemes we are going to mention use, in their iterative expression, the Jacobian matrix of function F and have, under the usual conditions, a convergence order 6. We will compare these methods, from the point of view of the results of the convergence order and computational efficiency, with the methods proposed in this paper that also have an order of 6 and only use ${\left[F^{\prime }\left(x\right)\right]}^{-1}$ in their expressions.

In [27], by using the weight function procedure, the authors designed a Jarratt-type method for solving nonlinear systems, denoted by $M{2}_6$, the iterative expression of which is

$$\left\{\begin{array}{ccc}{y}^{\left(k\right)}\hfill & =& x^{\left(k\right)}-\frac{2}{3}{\left[F^{\prime }\left(x^{\left(k\right)}\right)\right]}^{-1}F\left(x^{\left(k\right)}\right),k=0,1,\dots ,\hfill \\ z^{\left(k\right)}\hfill & =& x^{\left(k\right)}-\left(\frac{5}{8}I+\frac{3}{8}{\left({\left[F^{\prime }\left(y^{\left(k\right)}\right)\right]}^{-1}{F}^{\prime }\left(x^{\left(k\right)}\right)\right)}^2\right){\left[F^{\prime }\left(x^{\left(k\right)}\right)\right]}^{-1}F\left(x^{\left(k\right)}\right),\hfill \\ x^{(k+1)}\hfill & =& z^{\left(k\right)}-\left(\frac{-9}{4}I-\frac{15}{8}{\left[F^{\prime }\left(x^{\left(k\right)}\right)\right]}^{-1}{F}^{\prime }\left(y^{\left(k\right)}\right)+\frac{11}{8}{\left[F^{\prime }\left(y^{\left(k\right)}\right)\right]}^{-1}{F}^{\prime }\left(x^{\left(k\right)}\right)\right){\left[F^{\prime }\left(y^{\left(k\right)}\right)\right]}^{-1}F\left(z^{\left(k\right)}\right),\hfill \end{array}\right.$$

(2)

where I denotes the identity matrix of size $n\times n$. This method aims to evaluate the Jacobian matrix in two points and uses two inverse operators. These elements increase the number of operations per iteration.

In order to reduce the number of inverse operators, Narang et al. in [28], from a Chebyshev–Halley-type family, constructed a class of iterative schemes of the sixth order. One of its members, denoted by $M_{6,2}(1/2,0)$, has the following iterative expression:

$$\left\{\begin{array}{ccc}{y}^{\left(k\right)}\hfill & =& x^{\left(k\right)}-\frac{2}{3}{\left[F^{\prime }\left(x^{\left(k\right)}\right)\right]}^{-1}F\left(x^{\left(k\right)}\right),\hfill \\ z^{\left(k\right)}\hfill & =& x^{\left(k\right)}-\left(\frac{1}{2}G\left(x^{\left(k\right)}\right)\right)H\left(G\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 \\ x^{(k+1)}\hfill & =& z^{\left(k\right)}-\left(I+\frac{3}{2}G\left(x^{\left(k\right)}\right)\right){\left[F^{\prime }\left(x^{\left(k\right)}\right)\right]}^{-1}F\left(z^{\left(k\right)}\right),k=0,1,\dots ,\hfill \end{array}\right.$$

(3)

where $G\left(x^{\left(k\right)}\right)=I-{\left[F^{\prime }\left(x^{\left(k\right)}\right)\right]}^{-1}{F}^{\prime }\left(y^{\left(k\right)}\right)$ and $H\left(G\left(x^{\left(k\right)}\right)\right)=I-\frac{1}{4}G\left(x^{\left(k\right)}\right)+\frac{11}{8}{\left(G\left(x^{\left(k\right)}\right)\right)}^2$.

Behl et al. in [29], using the indeterminate parameter procedure, designed a family of iterative sixth-order methods for solving systems of nonlinear equations. One of its members, denoted by PM1, has the following iterative expression:

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

(4)

Finally, Yaseen and Zafar presented in [30] a Jarratt-type scheme of three steps for solving nonlinear systems, denoted by $FS6$, with sixth-order convergence and iterative expression:

$$\left\{\begin{array}{ccc}{y}^{\left(k\right)}\hfill & =& x^{\left(k\right)}-\frac{2}{3}{\left[F^{\prime }\left(x^{\left(k\right)}\right)\right]}^{-1}F\left(x^{\left(k\right)}\right),\hfill \\ z^{\left(k\right)}\hfill & =& x^{\left(k\right)}-\left(\frac{5}{8}{\left[U_k\right]}^{-1}+\frac{3}{8}{U}_k\right){\left[F^{\prime }\left(y^{\left(k\right)}\right)\right]}^{-1}F\left(x^{\left(k\right)}\right),\hfill \\ x^{(k+1)}\hfill & =& z^{\left(k\right)}-\left(\frac{-13}{2}I+\frac{9}{2}{\left[V_k\right]}^{-1}+3V_k\right){\left[F^{\prime }\left(x^{\left(k\right)}\right)\right]}^{-1}F\left(z^{\left(k\right)}\right),k=0,1,\dots ,\hfill \end{array}\right.$$

(5)

where $U_k={\left[F^{\prime }\left(y^{\left(k\right)}\right)\right]}^{-1}{F}^{\prime }\left(x^{\left(k\right)}\right)$ and $V_k={\left[F^{\prime }\left(x^{\left(k\right)}\right)\right]}^{-1}{F}^{\prime }\left(y^{\left(k\right)}\right)$.

The rest of the paper is organized as follows. In Section 2, we describe the nonlinear system obtained for predicting the shear strength of a reinforced concrete beam. The efficient method for estimating its solution is presented in Section 3, as well as its convergence order. Section 4 is devoted to the efficiency analysis. The numerical performance of our proposed methods are studied on academical problems and on the nonlinear shear model described in Section 2. Finally, some conclusions are exposed.

2. Problem Statement

In [3], the authors proposed this stress–strain relationship for concrete cracked in tension:

$${\sigma }_1=\left\{\begin{array}{cc}{E}_c{\epsilon }_1,\hfill & \hfill {\epsilon }_1\le {\epsilon }_{ct},\\ \frac{\alpha f_{ct}}{1+\sqrt{500{\epsilon }_1}},\hfill & \hfill {\epsilon }_1>{\epsilon }_{ct},\end{array}\right.$$

(6)

where ${\sigma }_1$ represents the contribution of tensile stresses in the concrete between the cracks or tension stiffening effect, ${\epsilon }_1$ is the principal tensile strain, $E_c$ being the modulus of elasticity of the concrete, ${\epsilon }_{ct}$ the strain related to the strength of the tensile, $f_{ct}$. Coefficient $\alpha$ is equal to 1.0 in cases of fast and non-cyclic loads and for deformed bars.

Regarding the concrete behaviour in compression, Vecchio and Collins formulated in [9], using the Modified Compression Field Theory (MCFT), the following relationship between diagonal compressive strain, ${\epsilon }_2$, and the diagonal (or principal) compressive stress, ${\sigma }_2$:

$$\begin{array}{c}{\sigma }_2=f_{2max}\left[2\left(\frac{{\epsilon }_2}{{\epsilon }_c}\right)-{\left(\frac{{\epsilon }_2}{{\epsilon }_c}\right)}^2\right],\hfill \\ \mathrm{with}{f}_{2max}=\frac{f_c}{0.8+170{\epsilon }_1}\le f_c,\hfill \end{array}$$

(7)

where ${\epsilon }_c$ is the compressive stress related to the compressive strength of concrete in a cylindrical test $f_c$, $f_{2max}$ is the maximum compressive stress in a diagonally cracked web and ${\epsilon }_1$ is the coexisting principal tensile stress.

In CFT procedures, a perfect bond between concrete and steel is assumed; in consequence, any deformation developed by the reinforcement is identical to the one experienced by the surrounding concrete in the same direction; thus, a single average strain tensor of the composite material is adopted. The following relationship is considered regarding the compatibility of the strains in the reinforcement and the diagonally stressed concrete:

$${tan}^2\theta =\frac{{\epsilon }_x-{\epsilon }_2}{{\epsilon }_t-{\epsilon }_2}=\frac{{\epsilon }_1-{\epsilon }_t}{{\epsilon }_1-{\epsilon }_x},$$

(8)

where ${\epsilon }_x$ is the mean longitudinal strain and ${\epsilon }_t$ is the mean transversal strain on the web of a beam oriented according to the orthogonal $x-t$ direction (see Figure 1). The strain ${\epsilon }_2$ is aligned in the direction of the compressive struts, at angle $\theta$ to the longitudinal axis $\left(x\right)$ of the beam. Moreover, due to strain tensor, the main tensile strain is

$${\epsilon }_1={\epsilon }_x+{\epsilon }_t+{\epsilon }_2.$$

(9)

On the other hand, in CFT models, the equilibrium between the external loads and the internal forces is governed by the following equations:

$${\sigma }_2=\frac{\nu }{b_{w}z}(tan\theta +cot\theta )-{\sigma }_1,$$

(10)

$$2A_{st}{\sigma }_{st}=({\sigma }_2{sin}^2\theta -{\sigma }_1{cos}^2\theta )b_{w}s,$$

(11)

$$4A_{sx}{\sigma }_{sx}+A_p{\sigma }_p=({\sigma }_2{cos}^2\theta -{\sigma }_1{sin}^2\theta )b_{w}z=\frac{\nu }{tan\theta }-{\sigma }_1{b}_{w}z,$$

(12)

where $\theta$ is the angle of the main tensile stress, z is the flexural lever arm, s is the stirrup spacing, $\nu$ is the internal shear force, and $b_w$ is the web width; $A_{sx}$, $A_{st}$ and $A_p$ are the cross-section surfaces for the longitudinal bars, the stirrup legs, and the prestressed reinforcement, respectively, and ${\sigma }_{sx}$, ${\sigma }_{st}$ and ${\sigma }_p$ are the related mean tensile stresses. The angles of the inclination of the principal strains coincide with the angles of the inclination of the principal stresses; this is known as EPA assumption or as Wagner’s hypothesis [31].

Regarding the stress–strain relationship of the steel reinforcement, beyond the type of steel to consider (such as, for example, mild steel or stainless steel), CFT methods mainly differ in terms of the treatment of the steel behavior [9,11,13]. In this work, one of the most recent approaches to steel behaviour is adopted: the RCFT, previously introduced in Section 1, which is based on the concept of an embedded bar model that takes into account the concrete tension stiffening effect between cracks. The latter theory allows us to apply, in the most general case, the following mean stress–strain model for each type of steel reinforcement of the beam (i.e., longitudinal reinforcement and transverse stirrups):

$$\begin{array}{c}{\sigma }_{s,i}=\left\{\begin{array}{ccc}{f}_{y,i}-\frac{{\lambda }_i{A}_{c,i}}{A_{s,i}}\frac{f_{ct}}{1+\sqrt{3.6M_i{\epsilon }_{s,i}}}\hfill & \mathrm{if}& {\epsilon }_{s,i}\ge {\epsilon }_{max,i}\\ E_s{\epsilon }_{s,i}\hfill & \mathrm{if}& {\epsilon }_{s,i}<{\epsilon }_{max,i}\end{array}\right.\hfill \\ \hfill \\ i=\{x,t\}\hfill \\ \mathrm{in}\mathrm{which}\hfill \\ {\epsilon }_{max,i}=\frac{f_{y,i}}{E_s}-\frac{\frac{{\lambda }_i{A}_{c,i}{f}_{ct}}{1+\sqrt{3.6M_i{\epsilon }_{max,i}}}}{E_s{A}_{s,i}}\hfill \\ M_i=\frac{{\lambda }_i{A}_{c,i}}{\sum \pi {\varphi }_i},\hfill \end{array}$$

(13)

where the subscripts x and t refer to the longitudinal and the transverse reinforcement, respectively (then, (13) actually involves two equations); $f_y$ is the steel yield stress, $E_s$ is the elastic modulus of the steel, ${\sigma }_{s,av}$ is the average tensile stress in the steel, ${\epsilon }_{s,av}$ is the average strain in the reinforcing bar, ${\epsilon }_{max,i}$ is the apparent yield strain (cf. [13]), M is the joint parameter, $A_s$ is the cross-section of the steel bars (longitudinal or transverse), $A_c$ is the area of concrete attached to the bar that participates in the tensile stiffening effect; this is usually considered equal to the rectangular area surrounding the bar of diameter $\varphi$ and over a distance no greater than $7.5\varphi$ from the center of the bar, and finally, ${\lambda }_i$ is the coefficient for fixing the numerical solvability of the steel constitutive model.

In the case of prestressed concrete members, the following two additional equations are required:

$${\epsilon }_p={\epsilon }_x+\Delta {\epsilon }_p,$$

(14)

$${\sigma }_p=\left\{\begin{array}{ccc}{E}_p{\epsilon }_p\hfill & ,& {\epsilon }_p\le \frac{f_{py}}{E_p},\\ f_{py}\hfill & ,& {\epsilon }_p>\frac{f_{py}}{E_p},\end{array}\right.$$

(15)

where (14) represents the strain compatibility, $\Delta {\epsilon }_p$ and ${\epsilon }_p$ being the strain imposed by the prestressing system and the strain of the prestressing strand, respectively, and Equation (15) represents the stress–strain relationship for the prestressing steel, $f_{py}$ and $E_p$ being its yield stress and elastic modulus, respectively.

Equation (13) is based on the concept of force equilibrium between a general section (or non-cracked section, where both the steel and the surrounding concrete contribute) and a cracked section (where only the reinforcement resists the internal forces; please see Figure 2). The greatest value of the area $A_c$ in order to preserve the solvability of the embedded steel constitutive model proposed by the RCFT (i.e., in order to preserve the internal equilibrium of forces, in such a way that as the concrete participation increases, the steel stress diminishes) is obtained by the application of the following coefficient [32]:

$${\lambda }_{max,i}=\frac{A_s·f_y}{A_c·f_{ct}}·\left(\frac{2}{3}+\frac{\sqrt{{\left(1+10.8·M·{ϵ}_y\right)}^3}}{48.6·M·{ϵ}_y}\right),$$

(16)

where the coefficient ${\lambda }_{max,i}$ represents the boundary of the solvability region for the embedded steel constitutive model in the i-direction (i.e., the maximum value of the coefficient $\lambda$ in the i-direction in order to preserve the solvability), and ${ϵ}_y$ is the strain corresponding to the steel yield stress (i.e., ${ϵ}_y=f_y/E_s$). For certain design cases, the previous boundary may lay within the design range prescribed by technical codes for the tension stiffening area, $A_c$.

In summary, for a given value of tensile principal strain in concrete, ${\epsilon }_1$, where such strain works as an input parameter, the shear model to predict the load–deformation response of a prestressed concrete beam is derived from the nonlinear system defined by (7)–(15), containing up to 10 equations (notice that (13) is actually two equations in turn) in the 10 unknowns ($\theta$, ${\epsilon }_x$, ${\epsilon }_t$, $\nu$, ${\epsilon }_2$, ${\sigma }_2$, ${\sigma }_{s,x}$, ${\sigma }_{s,t}$, ${\epsilon }_p$, and ${\sigma }_p$).

Two thousand solutions obtained from solving the nonlinear system of Equations (7)–(15) has been obtained from a set of input vectors uniformly generated. The range of the input parameters considered to this aim are presented in

Table 1. Ranges for input parameters of the nonlinear system of Equations (7)–(15), with $E_p$ = 190,000 MPa and $f_{py}=1674$ MPa.

Input	Range
$E_s\left(\mathrm{MPa}\right)$	[195,000, 205,000]
$f_y\left(\mathrm{MPa}\right)$	$\left[350,500\right]$
${\varphi }_x\left(\mathrm{mm}\right)$	$\left[6,40\right]$
${\varphi }_t\left(\mathrm{mm}\right)$	$\left[6,40\right]$
$\Delta {\epsilon }_p\left(-\right)$	$\left[0.10f_{py}/E_p,0.90f_{py}/E_p\right]$
$f_c\left({\mathrm{mm}}^2\right)$	$\left[25,50\right]$
$b_w\left(\mathrm{mm}\right)$	$\left[100,1000\right]$
$s\left(\mathrm{mm}\right)$	$\left[15{\varphi }_t,600\right]$
$A_p\left({\mathrm{mm}}^2\right)$	$\left[300,1200\right]$
${\lambda }_x\left(\mathrm{mm}\right)$	$\left[0.1{\lambda }_{max,x},0.9{\lambda }_{max,x}\right]$
${\lambda }_t\left(\mathrm{mm}\right)$	$\left[0.1{\lambda }_{max,t},0.9{\lambda }_{max,t}\right]$
${\epsilon }_1\left(-\right)$	$\left[0.0001,0.01\right]$

. These solutions were obtained using Newton’s method and considering the same initial approximation of all the cases.

3. Development and Convergence of the Method

By using the weight matrix function procedure, we present a class of three-step iterative methods with the following iterative expression:

$$\begin{array}{ccc}\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 z^{\left(k\right)}& =& y^{\left(k\right)}-G\left({\mu }^{\left(k\right)}\right)b{\left[F^{\prime }\left(x^{\left(k\right)}\right)\right]}^{-1}F\left(y^{\left(k\right)}\right)),\hfill \\ \hfill x^{(k+1)}& =& z^{\left(k\right)}-G\left({\mu }^{\left(k\right)}\right){\left[F^{\prime }\left(x^{\left(k\right)}\right)\right]}^{-1}\left(iF\left(z^{\left(k\right)}\right)+hF\left(y^{\left(k\right)}\right)\right),k=0,1,2,\dots ,\hfill \end{array}$$

(17)

where $\mu ={\left[F^{\prime }\left(x\right)\right]}^{-1}F\left(y\right)$ is the variable of the weight function G, and b, i, and h are free parameters.

On the other hand, with F being a sufficiently differentiable Fréchet function, we can regard $\xi +m\in {\mathbb{R}}^n$ as being in the neighbourhood of the zero of F, $\xi$. Using Taylor developments and $F^{\prime }\left(\xi \right)$ being nonsingular,

$$F(\xi +m)=F^{\prime }\left(\xi \right)\left[h+\sum _{q=2}^{p-1}{C}_q{m}^q\right]+O\left(m^p\right),$$

(18)

where $C_q={\displaystyle \frac{1}{q!}}{\left[F^{\prime }\left(\xi \right)\right]}^{-1}{F}^{\left(q\right)}\left(\xi \right)$ for $q\ge 2$. Also, $C_q{h}^q\in {\mathbb{R}}^n$, as $F^{\left(q\right)}\left(\xi \right)\in \mathcal{L}({\mathbb{R}}^n\times \cdots \times {\mathbb{R}}^n,{\mathbb{R}}^n)$ and ${\left[F^{\prime }\left(\xi \right)\right]}^{-1}\in \mathcal{L}\left({\mathbb{R}}^n\right)$. Therefore,

$$F(\xi +m)=F^{\prime }\left(\xi \right)\left[I+\sum _{q=2}^{p-1}q{C}_q{m}^{q-1}\right]+O\left(m^{p-1}\right),$$

(19)

being $q{C}_q{m}^{q-1}\in \mathcal{L}\left({\mathbb{R}}^n\right)$. For more details of this notation, see [33].

Indeed, following the notation introduced by Artidiello et al. in [34], the matrix function $G:X\to X$ can be defined in such a way that its Fréchet derivatives holds

(a)

$G^{\prime }\left(u\right)\left(v\right)=G_{1}uv,$ being $G^{\prime }:X\to \mathcal{L}\left(X\right)$, $G_1\in \mathbb{R}$

(b)

$G^{″}(u,v)\left(w\right)=G_{2}uvw,$ being $G_2:X\times X\to \mathcal{L}\left(X\right)$, $G_2\in \mathbb{R}$

• when $X={\mathbb{R}}^{n\times n}$ is the Banach space of real $n\times n$ matrices, and $\mathcal{L}\left(X\right)$ is the set of linear operators defined in X.

In the next result, we present the convergence of the family (17).

Theorem 1.

Let $F:D\subseteq {\mathbb{R}}^n\to {\mathbb{R}}^n$ be a sufficiently differential Fréchet function defined on a convex neighborhood D of $\xi \in {\mathbb{R}}^n$, and a zero of F. Also, let $G:{\mathbb{R}}^{n\times n}\to {\mathbb{R}}^{n\times n}$ be a sufficiently differentiable matrix function. Suppose that $F^{\prime }\left(\xi \right)$ is nonsingular, and that $x^{\left(0\right)}$ is a seed sufficiently close to ξ. Therefore, the sequence ${\left\{x^{\left(k\right)}\right\}}_{k\ge 0}$ from (17) converges to ξ with an order of convergence six if $b={\displaystyle \frac{1}{G_0}}$, $h=0$, $i={\displaystyle \frac{1}{G_0}}$, $G_1=-G_0$, and $|G_2|<\infty$, $G_0=G\left(I\right)$ and I being the identity matrix of size $n\times n$. In this case, the error equation is

$$\begin{array}{ccc}\hfill e^{(k+1)}& =& \left(24C_2^5-4G_2(37C_2^5-6C_3-6C_2{C}_3{C}_2^2{C}_3{C}_2+3C_3{C}_2{C}_3){\displaystyle \frac{1}{G_0}}\right.\hfill \\ & & +4G_2(32C_2^5{G}_0+C_2^5{G}_2-6G_0{C}_3{C}_2^3-6G_0{C}_2{C}_3{C}_2^2-6G_0{C}_2^2{C}_3{C}_2+3G_0{C}_3{C}_2{C}_3){\left({\displaystyle \frac{1}{G_0}}\right)}^2\hfill \\ & & \left.+2C_2^5{G}_0^2(3G_0-G_2){\left({\displaystyle \frac{1}{G_0}}\right)}^3\right){e^{\left(k\right)}}^6+O\left({e^{\left(k\right)}}^7\right),\hfill \end{array}$$

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

Proof. 

By means of the Taylor expansion of $F\left(x^{\left(k\right)}\right)$ and $F^{\prime }\left(x^{\left(k\right)}\right)$ about $\xi$, we obtain

$$F\left(x^{\left(k\right)}\right)=F^{\prime }\left(\xi \right)\left[e^{\left(k\right)}+\sum _{i=2}^6{C}_i{e^{\left(k\right)}}^i\right]+\mathcal{O}\left({e^{\left(k\right)}}^7\right),$$

and

$$F^{\prime }\left(x^{\left(k\right)}\right)=F^{\prime }\left(\xi \right)\left[I++\sum _{i=2}^{5}i{C}_i{e^{\left(k\right)}}^{i-1}\right]+\mathcal{O}\left({e^{\left(k\right)}}^6\right).$$

We can deduce that

$${\left[F^{\prime }\left(x^{\left(k\right)}\right)\right]}^{-1}=\left[I++\sum _{i=2}^5{X}_i{e^{\left(k\right)}}^{i-1}\right]{\left[F^{\prime }\left(\xi \right)\right]}^{-1}+\mathcal{O}\left({e^{\left(k\right)}}^6\right),$$

where $X_2=-2C_2$, $X_3=-3C_3+4C_2^2$, $X_4=-4C_4+6C_2{C}_3+6C_3{C}_2-8C_2^3$ and

$$\begin{array}{cc}\hfill X_5& =-5C_5+8C_2{C}_4-12C_2^2{C}_3+9C_3^2+8C_4{C}_2-12C_2{C}_3{C}_2+16C_2^4-12C_3{C}_2^2,\hfill \\ \hfill X_6& =-6C_6+10C_2{C}_5+12C_4{C}_3-18C_2{C}_3^2-18C_3{C}_2{C}_3+24C_2^3{C}_3+12C_3{C}_4\hfill \\ & \phantom{=i}-16C_2^2{C}_4+10C_5{C}_2-16C_2{C}_4{C}_2-18C_3^2{C}_2+24C_2^2{C}_3{C}_2-16C_4{C}_2^2\hfill \\ & \phantom{=i}+24C_2{C}_3{C}_2^2+24C_3{C}_2^3-32C_2^5.\hfill \end{array}$$

Then,

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

and

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

Moreover,

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

So,

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

and the expansion of the variable ${\mu }^{\left(k\right)}$ is

$$\begin{array}{cc}\hfill {\mu }^{\left(k\right)}& =-2C_2{e}^{\left(k\right)}+(6C_2^2-3C_3){e^{\left(k\right)}}^2+(-16C_2^3-4C_4+10C_2{C}_3+6C_3{C}_2){e^{\left(k\right)}}^3\hfill \\ & +(40C_2^4+9C_3^2-5C_5+14C_2{C}_4-28C_2^2{C}_3-15C_3{C}_2^2+8C_4{C}_2-18C_2{C}_3{C}_2){e^{\left(k\right)}}^4\hfill \\ & +(-96C_2^5-6C_6-30C_2{C}_3^2+18C_2{C}_5-40C_2^2{C}_4+72C_2^3{C}_3+36C_3{C}_2^3+12C_3{C}_4\hfill \\ & -12C_3^2{C}_2-24C_4{C}_2^2+12C_4{C}_3+10C_5{C}_2+42C_2{C}_3{C}_2^2\hfill \\ & -24C_2{C}_4{C}_2+48C_2^2{C}_3{C}_2-24C_3{C}_2{C}_3){e^{\left(k\right)}}^5+\mathcal{O}\left({e^{\left(k\right)}}^6\right).\hfill \end{array}$$

For the weight function G, we set

$$\begin{array}{cc}\hfill G\left({\mu }^{\left(k\right)}\right)& =G\left(I\right)+G_1({\mu }^{\left(k\right)}-I)+{\displaystyle \frac{1}{2}}{G}_2{({\mu }^{\left(k\right)}-I)}^2+\mathcal{O}{({\mu }^{\left(k\right)}-I)}^3,\hfill \end{array}$$

that is,

$$\begin{array}{cc}\hfill G\left({\mu }^{\left(k\right)}\right)& =G\left(I\right)-2C_2{G}_1{e}^{\left(k\right)}+(6C_2^2{G}_1-3C_3{G}_1+2C_2^2{G}_2){e^{\left(k\right)}}^2\hfill \\ & +(-16C_2^3{G}_1-4C_4{G}_1-12C_2^3{G}_2+10G_1{C}_2{C}_3+3G_2{C}_2{C}_3+6G_1{C}_3{C}_2\hfill \\ & +3G_2{C}_3{C}_2){e^{\left(k\right)}}^3+\mathcal{O}\left({e^{\left(k\right)}}^4\right).\hfill \end{array}$$

We denote that $S={\left[F^{\prime }\left(x^{\left(k\right)}\right)\right]}^{-1}F\left(y^{\left(k\right)}\right)$. So, its Taylor development can be expressed as

$$\begin{array}{cc}\hfill S& =C_2{e^{\left(k\right)}}^2+(-4C_2^2+2C_3){e^{\left(k\right)}}^3+(13C_2^3+3C_4-8C_2{C}_3-6C_3{C}_2){e^{\left(k\right)}}^4\hfill \\ & +(-38C_2^4-12C_3^2+4C_5-12C_2{C}_4+26C_2^2{C}_3+18C_3{C}_2^2-8C_4{C}_2+20C_2{C}_3{C}_2){e^{\left(k\right)}}^5\hfill \\ & +\mathcal{O}\left({e^{\left(k\right)}}^6\right).\hfill \end{array}$$

So,

$$\begin{array}{cc}\hfill z^{\left(k\right)}& =(C_2-b{C}_2{G}_0{e^{\left(k\right)}}^2+(-2C_2^2+2C_3+4b{C}_3{G}_0){e^{\left(k\right)}}^3\hfill \\ & +(4C_2^3+3G_0-13b{C}_2^3{G}_0-3b{C}_4{G}_0-14b{C}_2^3{G}_2-2b{C}_2^3{G}_2-4C_2{C}_3\hfill \\ & +8bG{C}_2{C}_3+4b{G}_1{C}_2{C}_3-3C_3{C}_2+6b{G}_0{C}_2+3b{G}_1{C}_3{C}_2){e^{\left(k\right)}}^4\hfill \\ & +(-8C_2^4-6C_3^2+38b{C}_2^4{G}_0+12b{C}_3^2{G}_0+66b{C}_2^4{G}_1+6b{C}_3^2{G}_1+20b{C}_2^4{G}_2\hfill \\ & +4C_5-4b{G}_0{C}_5-6C_2{C}_4+12b{G}_0{C}_2{C}_4+6b{G}_1{C}_2{C}_4+8C_2^2{C}_3-26b{G}_0{C}_2^2{C}_3\hfill \\ & -28b{G}_1{C}_2^2{C}_3-4b{G}_2{C}_2^2{C}_3+6C_3{C}_2^2-18b{G}_0{C}_3{C}_2^2-18b{G}_1{C}_3{C}_2^2-3b{G}_2{C}_3{C}_2^2\hfill \\ & -4C_4{C}_2+8b{G}_0{C}_4{C}_2+4b{G}_1{C}_4{C}_2+6C_2{C}_3{C}_2-20b{G}_0{C}_2{C}_3{C}_2-22b{G}_1{C}_2{C}_3{C}_2\hfill \\ & -3b{G}_2{C}_2{C}_3{C}_2){e^{\left(k\right)}}^5+\mathcal{O}\left({e^{\left(k\right)}}^6\right),\hfill \end{array}$$

and

$$\begin{array}{cc}\hfill F\left(z^{\left(k\right)}\right)& =(C_2-b{C}_2{G}_0){e^{\left(k\right)}}^2+(-2C_2^2+2C_3+4b{C}_2{G}_0-2b{C}_3{G}_0-2b{C}_2^2{G}_1){e^{\left(k\right)}}^3\hfill \\ & +(5C_2^3+3C_4-15b{C}_2^3{G}_0-3b{C}_4{G}_0+b^2{c}_2^3{G}_0^2-14b{C}_2^3{G}_1-2b{C}_2^3{G}_2-4C_2{C}_3\hfill \\ & +8b{G}_0{C}_2{C}_3+4b{G}_1{C}_2{C}_3-3C_3{C}_2+6b{G}_0{C}_3{C}_2+3b{G}_1{C}_3{C}_2){e^{\left(k\right)}}^4\hfill \\ & +(-8C_2^4-6C_3^2+38b{C}_2^4{G}_0+12b{C}_3^2{G}_0+66b{C}_2^4{G}_1+6b{C}_3^2{G}_1+20b{C}_2^4{G}_2+4C_5\hfill \\ & -4b{G}_0{C}_5-6C_2{C}_4+12b{G}_0{C}_2{C}_4+6b{G}_1{C}_2{C}_4+8C_2^2{C}_3-26b{G}_0{C}_2^2{C}_3\hfill \\ & -28b{G}_1{C}_2^2{C}_3-4b{G}_2{C}_2^2{C}_3+6C_3{C}_2^2-18b{G}_0{C}_3{C}_2^2-18b{G}_1{C}_3{C}_2^2-3b{G}_2{C}_3{C}_2^2\hfill \\ & -4C_4{C}_2+8b{G}_0{C}_4{C}_2+4b{G}_1{C}_4{C}_2+6C_2{C}_3{C}_2-20b{G}_0{C}_3{C}_2-22b{G}_1{C}_2{C}_3{C}_2\hfill \\ & -3b{G}_2{C}_2{C}_3{C}_2){e^{\left(k\right)}}^5+\mathcal{O}\left({e^{\left(k\right)}}^6\right).\hfill \end{array}$$

Now, we denote $Sc={\left[F^{\prime }\left(x^{\left(k\right)}\right)\right]}^{-1}F\left(z^{\left(k\right)}\right)$. Therefore, its Taylor expansion is

$$\begin{array}{cc}\hfill Sc& =(C_2-b{C}_2{G}_0){e^{\left(k\right)}}^2+(-4C_2^2+2C_3+6b{C}_2^2{G}_0-2b{C}_3{G}_0+2b{C}_2^2{G}_1){e^{\left(k\right)}}^3\hfill \\ & +(13C_2^3+3C_4-27b{C}_2^3{G}_0-3b{C}_4{G}_0+b^2{C}_2^3{G}_0^2-18b{C}_2^3{G}_1-2b{C}_2^3{G}_2-8C_2{C}_3\hfill \\ & +12b{G}_0{C}_2{C}_3+4b{G}_1{C}_2{C}_3-6C_3{C}_2+9bG{C}_3{C}_2+9b{G}_0{C}_3{C}_2+3b{G}_1{C}_3{C}_2){e^{\left(k\right)}}^4\hfill \\ & +(-34C_2^4-12C_3^2+92b{C}_2^{4}G+18b{C}_3^2{G}_0-2b^2{C}_2^4{G}_0^2+102b{C}_2^4{G}_1+6b{C}_3^2{G}_1\hfill \\ & +24b{C}_2^4{G}_2+4C_5-4b{G}_0{C}_5-12C_2{C}_4+18b{G}_0{C}_2{C}_4+6b{G}_1{C}_2{C}_4+24C_2^2{C}_3\hfill \\ & -50b{G}_0{C}_2^2{C}_3-36b{G}_1{C}_2^2{C}_3-4b{G}_2{C}_2^2{C}_3+18C_3{C}_2^2-36b{G}_0{C}_3{C}_2^2-24b{G}_1{C}_3{C}_2^2\hfill \\ & -3b{G}_2{C}_3{C}_2^2-8C_4{C}_2+12b{G}_0{C}_4{C}_2+4b{G}_1{C}_4{C}_2+18C_2{C}_3{C}_2-38b{G}_0{C}_2{C}_3\hfill \\ & -28b{G}_1{C}_2{C}_3{C}_2-3b{G}_2{C}_2{C}_3{C}_2){e^{\left(k\right)}}^5+\mathcal{O}\left({e^{\left(k\right)}}^6\right).\hfill \end{array}$$

If $Ss=i·Sc+h·S$, then it is expanded as

$$\begin{array}{cc}\hfill Ss& =(C_{2}h+C_{2}i-b{C}_{2}Gi){e^{\left(k\right)}}^2\hfill \\ & +(-4C_2^{2}h+2C_{3}h-4C_2^{2}i+2C_{3}i+6b{C}_2^2{G}_{0}i-2b{C}_3{G}_{0}i+2b{C}_2^2{G}_{1}i){e^{\left(k\right)}}^3\hfill \\ & +(13C_2^{3}h+3C_{4}h+13C_2^{3}i+3C_{4}i-27b{C}_2^3{G}_{0}i-3b{C}_4{G}_{0}i+b^2{C}_2^3{G}_0^{2}i-18b{C}_2^3{G}_{0}i\hfill \\ & -2b{C}_2^3{G}_{2}i-8h{C}_2{C}_3-8i{C}_2{C}_3+12b{G}_{0}i{C}_2{C}_3+4b{G}_{1}i{C}_2{C}_3-6h{C}_3{C}_2-6i{C}_3{C}_2\hfill \\ & +9b{G}_{0}i{C}_3+3b{G}_{1}i{C}_3{C}_2){e^{\left(k\right)}}^4\hfill \\ & +(-38C_2^{4}h-12C_{3}h-34C_2^{4}i-12C_3^{2}i+92b{C}_2^4{G}_{0}i+18b{C}_3^2{G}_{0}i-2b{C}_2^4{G}_0^{2}i+102b{C}_2^4{G}_{0}i\hfill \\ & +6b{C}_3^2{G}_{0}i+24b{C}_2^4{G}_{2}i+4h{C}_5+4i{C}_5-4bGii{C}_5-12h{C}_2{C}_4-12i{C}_2{C}_4+18bGi{C}_2{C}_4\hfill \\ & +6b{G}_{1}i{C}_2{C}_4+26h{C}_2{C}_3+24i{C}_2^2{C}_3-50b{G}_{0}i{C}_2^2{C}_3-36b{G}_{1}i{C}_2^2{C}_3-4b{G}_{2}i{C}_2^2{C}_3\hfill \\ & +18G_0{C}_3{C}_2^2+18i{C}_3{C}_2^2-36bGi{C}_3{C}_2^2-24b{G}_{1}i{C}_3{C}_2^2-3b{G}_{2}i{C}_3{C}_2^2-8h{C}_4{C}_2-8i{C}_4{C}_2\hfill \\ & +12b{G}_{0}i{C}_4{C}_2+4b{G}_{1}i{C}_4{C}_2+20h{C}_2{C}_3{C}_2+18i{C}_2{C}_3{C}_2-38bGi{C}_2{C}_3{C}_2-28b{G}_{0}i{C}_2{C}_3{C}_2\hfill \\ & -3b{G}_{2}i{C}_2{C}_3{C}_2){e^{\left(k\right)}}^5+\mathcal{O}\left({e^{\left(k\right)}}^6\right).\hfill \end{array}$$

Then, the error equation is

$$\begin{array}{cc}\hfill e^{(k+1)}& =(C_2-b{C}_2{G}_0-C_2{G}_{0}h-C_2{G}_{0}i+b{C}_2{G}_0^{2}i){e^{\left(k\right)}}^2\hfill \\ & +(-2C_2^2+2d+4b{C}_2^2{G}_0-2bd{G}_0+2b{C}_2^2{G}_0+4C_2^2{G}_{0}h-2d{G}_{0}h\hfill \\ & +2C_2^2{G}_{0}h+4C_2^2{G}_{0}i-2d{G}_{0}i-6b{C}_2^2{G}_0^{2}i+2bd{G}_0^{2}i+2C_2^2{G}_{1}i-4b{C}_2^2{G}_0{G}_{1}i){e^{\left(k\right)}}^3\hfill \\ & +(4C_2^3+3C_4-13b{C}_2^3{G}_0-3b{C}_4{G}_0-14b{C}_2^3{G}_1-2b{C}_2^3{G}_2-13C_2^3{G}_{0}h-3C_4{G}_{0}h\hfill \\ & -14C_2^3{G}_{1}h-2C_2^3{G}_{2}h-13C_2^{3}Gi-3C_4{G}_{0}i+27b{C}_2^3{G}_0^{2}i+3b{C}_4{G}_0^{2}i-b^2{C}_2^3{G}_0^{3}i\hfill \\ & -14C_2^3{G}_{1}i+36b{C}_2^3{G}_0{G}_{1}i+4b{C}_2^3{G}_1^{2}i-2C_2^3{G}_{2}ii+4b{C}_2^3{G}_0{G}_{2}ii-4C_2{C}_3+8b{G}_0{C}_2{C}_3\hfill \\ & +4b{G}_1{C}_2{C}_3+8G_{0}h{C}_2{C}_3+4G_{1}h{C}_2{C}_3+8G_{0}i{C}_2{C}_3-12b{G}_0^{2}i{C}_2{C}_3+4G_{1}i{C}_2{C}_3\hfill \\ & -8b{G}_0{G}_{1}i{C}_2{C}_3-3C_3{C}_2+6b{G}_0{C}_3{C}_2+3b{G}_1{C}_3+6G_{0}h{C}_3{C}_2+3G_{1}h{C}_3{C}_2\hfill \\ & +6G_{0}i{C}_3{C}_2-9b{G}_0^{2}i{C}_3{C}_2+G_{1}i{C}_3{C}_2-6b{G}_0{G}_{1}i{C}_3{C}_2){e^{\left(k\right)}}^4+M5{e^{\left(k\right)}}^5\hfill \\ & +M6{e^{\left(k\right)}}^6+\mathcal{O}\left({e^{\left(k\right)}}^7\right).\hfill \end{array}$$

By fixing $b={\displaystyle \frac{1}{G_0}}$, $h=0$, $i={\displaystyle \frac{1}{G_0}}$, and $G_1=-G_0$, the error equation becomes

$$\begin{array}{cc}\hfill e^{(k+1)}& =\left(24C_2^5-4G_2(37C_2^5-6C_3-6C_2{C}_3{C}_2^2{C}_3{C}_2+3C_3{C}_2{C}_3){\displaystyle \frac{1}{G_0}}\right.\hfill \\ & +4G_2(32C_2^5{G}_0+C_2^5{G}_2-6G_0{C}_3{C}_2^3-6G_0{C}_2{C}_3{C}_2^2-6G_0{C}_2^2{C}_3{C}_2\hfill \\ & +3G_0{C}_3{C}_2{C}_3){\left({\displaystyle \frac{1}{G_0}}\right)}^2\left.+2C_2^5{G}_0^2(3G_0-G_2){\left({\displaystyle \frac{1}{G_0}}\right)}^3\right){e^{\left(k\right)}}^6+\mathcal{O}\left({e^{\left(k\right)}}^7\right).\hfill \end{array}$$

With this, the proof is finished. □

Let us notice that the order of convergence of this class of iterative methods can be increased up to 7 for specific values of $G_2$, depending on $G_0$ also being free. However, in order to reduce the computational cost, we set $G_2=0$, $G_0=I$, and therefore the matrix weight function to be used in the iterative expression is

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

Therefore, the family is reduced to an iterative method of only order 6, denoted by O6, the iterative expression of which is

$$\begin{array}{ccc}\hfill y^{\left(k\right)}& =& x^{\left(k\right)}-{\left[F^{\prime }\left(x^{\left(k\right)}\right)\right]}^{-1}F\left(x^{\left(k\right)}\right),k=0,1,2,\dots ,\hfill \\ \hfill z^{\left(k\right)}& =& y^{\left(k\right)}-\left[2I-{\left[F^{\prime }\left(x^{\left(k\right)}\right)\right]}^{-1}{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 \\ \hfill x^{(k+1)}& =& z^{\left(k\right)}-\left[2I-{\left[F^{\prime }\left(x^{\left(k\right)}\right)\right]}^{-1}{F}^{\prime }\left(y^{\left(k\right)}\right)\right]{\left[F^{\prime }\left(x^{\left(k\right)}\right)\right]}^{-1}F\left(z^{\left(k\right)}\right).\hfill \end{array}$$

Let us remark that this scheme has especially good properties, due to the existence of only one inverse operator. This yields that all the linear systems to be solved per iteration have the same coefficient matrix and therefore, the computational cost can be reduced by its LU factorization and the solution of several triangular linear systems. This is discussed in depth in the next section, in comparison with the introduced known procedures.

4. Efficiency Indices

To compare the iterative methods used, we use the computational efficiency index, $CI$, defined as [35]

$$\begin{array}{c}\hfill CI={\rho }^{{\displaystyle \frac{1}{d+op}}},\end{array}$$

where d is the number of functional evaluations and $op$ is the number of products/quotients per iteration.

In each iteration, five linear systems are solved with the same coefficient matrix, there are two matrix–vector products and, with respect to functional evaluations, we have two evaluations of Jacobian matrices and three of functions. The computational cost of method O6 is

$${\displaystyle \frac{1}{3}}{n}^3+9n^2+{\displaystyle \frac{8}{9}}n,$$

(20)

In

Table 2. Comparisons of CI.

Method	CI
Newton	$2^{\frac{1}{\frac{1}{3}{n}^3+2n^2+\frac{2}{3}n}}$
O6	$6^{\frac{1}{\frac{1}{3}{n}^3+9n^2+\frac{8}{9}n}}$
PM1	$6^{\frac{1}{\frac{2}{3}{n}^3+11n^2+\frac{4}{3}n}}$
$M_{2,6}$	$6^{\frac{1}{\frac{2}{3}{n}^3+10n^2+\frac{4}{3}n}}$
$M_{6,2}(1/2,0)$	$6^{\frac{1}{\frac{1}{3}{n}^3+12n^2+\frac{5}{3}n}}$
FS6	$6^{\frac{1}{\frac{2}{3}{n}^3+10n^2+\frac{4}{3}n}}$

, the rest of the CI corresponding to the comparison methods are exposed. The way in which they have been calculated is similar to that of the O6 method.

The results are represented in the semi-logarithmic scale; see Figure 3 for a better visualization of the differences between the indices (CI) for the methods used and several sizes $\left(n\right)$ of the systems.

In Figure 3a, we can observe that, for $2\le n\le 7$, the best CI index corresponds to the Newton method, O6 being the best for $n\ge 8$. In Figure 3b, we can check that for bigger systems, $n\ge 10$; the best CI remains as O6.

5. Numerical Performance

We analyze the performance of the methods described above to check their efficiency and compare it with other known methods. The results from Table 3, Table 4, Table 5 and Table 6 correspond to the calculations made with Matlab R2022b, by using variable precision arithmetics with 1200 digits of mantissa, on a PC equipped with an $Intel$ Core™i5-5200U CPU 2.20GHz. In all the tables, we show the residual errors $\parallel x^{(k+1)}-x^{\left(k\right)}\parallel$ and $\parallel F\left(x^{(k+1)}\right)\parallel$ of the last iteration satisfying the stopping criterium $\parallel x^{(k+1)}-x^{\left(k\right)}\parallel <{10}^{-300}$ or $\parallel F\left(x^{(k+1)}\right)\parallel <{10}^{-300}$, and the CPU time obtained as the mean of 20 executions (e-time). Moreover, a computational estimation of the order of convergence is obtained by the means of ACOC, introduced as

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

(21)

5.1. Example

We consider the nonlinear system, $F_1\left(x\right)={(f_1\left(x\right),f_2\left(x\right),\dots ,f_n\left(x\right))}^T=0$, such that

$$f_i\left(x\right)=x_i-cos\left(2x_i-\sum _{j=1}^4{x}_j-x_i\right),i=1,2,3,4\dots ,20,$$

(22)

with seed $x^{\left(0\right)}={(0.75,0.75,\dots ,0.75)}^T$, and in this case, $\alpha \approx {(0.519,0.519,\dots ,0.519)}^T$.

Table 3. Numerical results for Example 5.1.

Method	Iteration	$\parallel {\mathit{x}}^{(\mathit{k}+1)}-{\mathit{x}}^{\left(\mathit{k}\right)}\parallel$	$\parallel \mathit{F}\left({\mathit{x}}^{(\mathit{k}+1)}\right)\parallel$	$\mathit{\rho }$	e-Time
Newton	8	$3.1586\times {10}^{-160}$	$2.297\times {10}^{-320}$	2.0	0.98
O6	4	$3.4133\times {10}^{-217}$	0.0	6.0	0.99
FS6	4	$1.682\times {10}^{-201}$	$1.614\times {10}^{-1207}$	6.0	1.02
PM1	4	$6.0584\times {10}^{-186}$	$4.036\times {10}^{-1115}$	6.0	0.99
$M{2}_6$	4	$4.7636\times {10}^{-127}$	$1.891\times {10}^{-635}$	5.0	1.00
$M_{6,2}(1/2,0)$	4	$1.4631\times {10}^{-189}$	$4.603\times {10}^{-1137}$	6.0	1.02

Table 3. Numerical results for Example 5.1.

Method	Iteration	$\parallel {\mathit{x}}^{(\mathit{k}+1)}-{\mathit{x}}^{\left(\mathit{k}\right)}\parallel$	$\parallel \mathit{F}\left({\mathit{x}}^{(\mathit{k}+1)}\right)\parallel$	$\mathit{\rho }$	e-Time
Newton	8	$3.1586\times {10}^{-160}$	$2.297\times {10}^{-320}$	2.0	0.98
O6	4	$3.4133\times {10}^{-217}$	0.0	6.0	0.99
FS6	4	$1.682\times {10}^{-201}$	$1.614\times {10}^{-1207}$	6.0	1.02
PM1	4	$6.0584\times {10}^{-186}$	$4.036\times {10}^{-1115}$	6.0	0.99
$M{2}_6$	4	$4.7636\times {10}^{-127}$	$1.891\times {10}^{-635}$	5.0	1.00
$M_{6,2}(1/2,0)$	4	$1.4631\times {10}^{-189}$	$4.603\times {10}^{-1137}$	6.0	1.02

In Table 3, it can be observed that the number of iterations of all the sixth-order schemes are equal and the time is very similar in all of the methods; however, the best residual is obtained by the proposed scheme, O6. The ACOC estimates the theoretical order of the convergence accurately in all the cases.

5.2. Example

The second example is given by $F_2\left(x\right)={(g_1\left(x\right),g_2\left(x\right),\dots ,g_n\left(x\right))}^T=0$, such that

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

(23)

with seed $x^{\left(0\right)}={(1,1,\dots ,1)}^T$ and $\alpha \approx {(9.376,9.376,\dots ,9.376)}^T$.

Table 4. Numerical results for Example 5.2.

Method	Iterations	$\parallel {\mathit{x}}^{(\mathit{k}+1)}-{\mathit{x}}^{\left(\mathit{k}\right)}\parallel$	$\parallel \mathit{F}\left({\mathit{x}}^{(\mathit{k}+1)}\right)\parallel$	$\mathit{\rho }$	e-Time (Sec)
Newton	11	$1.1642\times {10}^{-199}$	$3.409\times {10}^{-401}$	2.0	9.99
O6	5	$3.3111\times {10}^{-100}$	$1.032\times {10}^{-608}$	6.0	10.24
FS6	5	$5.1171\times {10}^{-73}$	$1.310\times {10}^{-445}$	6.0	11.00
PM1	6	$2.169\times {10}^{-291}$	$3.749\times {10}^{-1755}$	6.0	10.58
$M{2}_6$	6	$5.3026\times {10}^{-198}$	$4.327\times {10}^{-996}$	6.0	11.07
$M_{6,2}(1/2,0)$	6	$6.2633\times {10}^{-289}$	$1.325\times {10}^{-1206}$	6.0	10.79

Table 4. Numerical results for Example 5.2.

Method	Iterations	$\parallel {\mathit{x}}^{(\mathit{k}+1)}-{\mathit{x}}^{\left(\mathit{k}\right)}\parallel$	$\parallel \mathit{F}\left({\mathit{x}}^{(\mathit{k}+1)}\right)\parallel$	$\mathit{\rho }$	e-Time (Sec)
Newton	11	$1.1642\times {10}^{-199}$	$3.409\times {10}^{-401}$	2.0	9.99
O6	5	$3.3111\times {10}^{-100}$	$1.032\times {10}^{-608}$	6.0	10.24
FS6	5	$5.1171\times {10}^{-73}$	$1.310\times {10}^{-445}$	6.0	11.00
PM1	6	$2.169\times {10}^{-291}$	$3.749\times {10}^{-1755}$	6.0	10.58
$M{2}_6$	6	$5.3026\times {10}^{-198}$	$4.327\times {10}^{-996}$	6.0	11.07
$M_{6,2}(1/2,0)$	6	$6.2633\times {10}^{-289}$	$1.325\times {10}^{-1206}$	6.0	10.79

5.3. Example

Let us define now the nonlinear system $F_3\left(x\right)={(h_1\left(x\right),h_2\left(x\right),\dots ,h_n\left(x\right))}^T=0$, such that

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

(24)

with seed $x^{\left(0\right)}={(0.5,0.5,\dots ,0.5)}^T$, $n=20$, and $\alpha \approx {(0.1758,0.1758,\dots ,0.1758)}^T$.

Of note, in Table 5, O6 and FS6 provide a solution satisfying the stopping criterium in a lower or equal number of iterations than the rest of the schemes. Indeed, the value of the residual errors in O6 and FS6 highly improve that of Newton’s. This is the reason why their residuals are not as close to zero as those of the other schemes.

Table 5. Numerical results for Example 5.3.

Method	Iterations	$\parallel {\mathit{x}}^{(\mathit{k}+1)}-{\mathit{x}}^{\left(\mathit{k}\right)}\parallel$	$\parallel \mathit{F}\left({\mathit{x}}^{(\mathit{k}+1)}\right)\parallel$	$\mathit{\rho }$	e-Time (Sec)
Newton	10	$1.2449\times {10}^{-154}$	$1.322\times {10}^{-307}$	2.0	1.22
O6	5	$1.3563\times {10}^{-218}$	$2.414\times {10}^{-1207}$	6.0	1.24
FS6	4	$4.4455\times {10}^{-58}$	$4.283\times {10}^{-344}$	6.0	1.20
PM1	5	$1.2252\times {10}^{-218}$	$8.687\times {10}^{-1208}$	6.0	1.27
$M{2}_6$	5	$6.2256\times {10}^{-173}$	$5.016\times {10}^{-861}$	5.0	1.39
$M_{6,2}(1/2,0)$	5	$2.5983\times {10}^{-252}$	$3.704\times {10}^{-1208}$	6.0	1.39

Table 5. Numerical results for Example 5.3.

Method	Iterations	$\parallel {\mathit{x}}^{(\mathit{k}+1)}-{\mathit{x}}^{\left(\mathit{k}\right)}\parallel$	$\parallel \mathit{F}\left({\mathit{x}}^{(\mathit{k}+1)}\right)\parallel$	$\mathit{\rho }$	e-Time (Sec)
Newton	10	$1.2449\times {10}^{-154}$	$1.322\times {10}^{-307}$	2.0	1.22
O6	5	$1.3563\times {10}^{-218}$	$2.414\times {10}^{-1207}$	6.0	1.24
FS6	4	$4.4455\times {10}^{-58}$	$4.283\times {10}^{-344}$	6.0	1.20
PM1	5	$1.2252\times {10}^{-218}$	$8.687\times {10}^{-1208}$	6.0	1.27
$M{2}_6$	5	$6.2256\times {10}^{-173}$	$5.016\times {10}^{-861}$	5.0	1.39
$M_{6,2}(1/2,0)$	5	$2.5983\times {10}^{-252}$	$3.704\times {10}^{-1208}$	6.0	1.39

Regarding the applied problem described in Section 2, the underlying data of the nonlinear shear model of a reinforced concrete beam are provided by random values with few digits inside the prescribed ranges by technical standards for structural concrete; moreover, some of these parameters were fixed taking into consideration the solvability region of the adopted steel constitutive model. The stopping criterium is $\parallel x^{(k+1)}-x^{\left(k\right)}\parallel <{10}^{-6}$ or $\parallel F\left(x^{(k+1)}\right)\parallel <{10}^{-6}$. The initial estimation used is $\theta =34$, ${\epsilon }_x=0.0001$, ${\epsilon }_t=200,000$, $\nu =0.0001$, ${\epsilon }_2=0.0001$, ${\sigma }_2=0.0001$, ${\sigma }_{s,x}=200$, ${\sigma }_{s,t}=7$, ${\epsilon }_p=200$ and ${\sigma }_p=100$. The results provided by the new and existing schemes appear in Table 6. The ACOC does not appear in this table, as it yields to unstable data in all cases.

Table 6. Problem statement Section 2.

Method	Iterations	$\parallel {\mathit{x}}^{(\mathit{k}+1)}-{\mathit{x}}^{\left(\mathit{k}\right)}\parallel$	$\parallel \mathit{F}\left({\mathit{x}}^{(\mathit{k}+1)}\right)\parallel$	e-Time (Sec)
Newton	5	0.0342	$2.448\times {10}^{-10}$	18.3242
O6	3	9.0862	$2.253\times {10}^{-29}$	21.6703
FS6	3	218.18	$4.071\times {10}^{-16}$	22.1258
PM1	4	3.8554	$3.138\times {10}^{-28}$	21.9594
$M{2}_6$	4	7.3926	$9.621\times {10}^{-21}$	28.0727
$M_{6,2}(1/2,0)$	3	0.0423	$3.710\times {10}^{-31}$	22.4797

Table 6. Problem statement Section 2.

Method	Iterations	$\parallel {\mathit{x}}^{(\mathit{k}+1)}-{\mathit{x}}^{\left(\mathit{k}\right)}\parallel$	$\parallel \mathit{F}\left({\mathit{x}}^{(\mathit{k}+1)}\right)\parallel$	e-Time (Sec)
Newton	5	0.0342	$2.448\times {10}^{-10}$	18.3242
O6	3	9.0862	$2.253\times {10}^{-29}$	21.6703
FS6	3	218.18	$4.071\times {10}^{-16}$	22.1258
PM1	4	3.8554	$3.138\times {10}^{-28}$	21.9594
$M{2}_6$	4	7.3926	$9.621\times {10}^{-21}$	28.0727
$M_{6,2}(1/2,0)$	3	0.0423	$3.710\times {10}^{-31}$	22.4797

However, the best methods in terms of the number of iteration are O6, FS6, and $M_{6,2}(1/2,0)$, all with three iterations. Among the sixth-order methods, the lowest e-time corresponds to our proposed scheme, O6. Although with this initial estimation, the e-time of Newton’s method is the best, small changes in some of the coordinates of the seed yields to better results of O6 than Newton’s scheme. This good performance allows us to assure the reliance and robustness of our proposed procedure.

6. Conclusions

In this article, we have developed a vectorial parametric family of numerical methods of the sixth order to solve nonlinear systems. In particular, it is applied on a constitutive equation of reinforced concrete (6). The order of the convergence of the new class (O6) is proven, and a particular member of the family is selected with better computational properties, as only one inverse operator is needed. Its efficiency is compared to other existing methods with the same order of convergence, and also with Newton’s scheme, in terms of the computational efficiency index. For the size of the system $n\ge 8$, the proposed method, O6, gives the best results. In the numerical tests, all the comparison procedures need the same or more iterations and achieve lower precision results in the same or shorter execution time to achieve the required tolerance. This confirms the accuracy, robustness, and applicability of the proposed scheme.