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

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 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 preexisting 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. Shear design procedures for reinforced concrete that determine the inclination of such cracks by considering the strains in the diagonally stressed concrete, as well as in the longitudinal and transverse reinforcement, are known as Compression Field Theories (CFT) [1,2].

As justified at Section 2, CFT mechanical models involve several types of nonlinearities, among other reasons, due to the constitutive relationships of the reinforced concrete. The implementation of these models usually requires the application of numerical methods for solving the corresponding nonlinear equations, such as, for example, Newton-type methods [3]. In fact, the correct solver for solving a nonlinear problem is often a choice between computational cost and accuracy [4,5,6,7]. 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 at Section 2.

Solving systems of nonlinear equations is an important problem in science and engineering. The objective is to find the roots of the nonlinear system $F\left(x\right)=0$, being F 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 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 a Jacobian matrix of F at k-th iteration.

Several Newton type procedures, by using different techniques, have been published in the last years. Their main aim is accelerating the convergence or increasing their efficiency. In the last section 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, convergence order 6. We will compare these methods, from the point of view of results, convergence order and computational efficiency, with the methods proposed in this paper that also have order 6 and use $F^{\prime }\left(x\right)$ in their expressions.

In [8], by using the weight function procedure, the authors designed a Jarratt-type method for solving nonlinear systems, denoted by $M{2}_6$, whose iterative expression 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 needs 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 [9] from a Chebyshev-Halley-type family, constructed a class of iterative schemes of 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 [10], using the indeterminate parameter procedure, designed a family of iterative sixth-order methods for solving systems of nonlinear equations, one of whose members, denoted by PM1, has the 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 [11] 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.

2. Problem Statement

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

$${\sigma }_1=\left\{\begin{array}{ccc}\hfill & E_c{\epsilon }_1,\hfill & \hfill {\epsilon }_1\le {\epsilon }_{ct},\\ \hfill & \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, being $E_c$ 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 case of fast and non-cyclic loads and for deformed bars.

Regarding to the concrete behaviour in compression, Vecchio and Collins formulated in [12], inside 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 strain related to the compressive strength of concrete in a cylinder test $f_c$, ${\epsilon }_1$ is the coexisting principal tensile strain, and $f_{2max}$ is the maximum compressive stress in a diagonally cracked web.

In CFT procedures, 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 by 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$ directions (see Figure 1). The strain ${\epsilon }_2$ is aligned in the direction of the compressive struts, at angle $\theta$ to the longitudinal axis (x) 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 inclination of principal strains coincide with the angles of inclination of principal stresses; this is known as EPA assumption or as Wagner’s hypothesis [2,13].

The main difference among CFT methods lies in the treatment of the steel behavior [12,14,15]. In this work one of the most recent approaches to the steel behaviour is adopted, the so-called Refined Compression Field Theory (RCFT) [15,16], which is based on the concept of an embedded bar model that takes into account the concrete tension stiffening effect between cracks. This last theory allows us to implement, in the most general case, the following average stress-strain model for each kind of steel reinforcement in the beam (that is, the transverse stirrups and the longitudinal reinforcement):

$$\begin{array}{c}{\sigma }_{s,i}=\left\{\begin{array}{ccc}{f}_{y,i}-\frac{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\in \{x,t\},\mathrm{and\; in\; which}\hfill \\ \hfill \\ {\epsilon }_{max,i}=\frac{f_{y,i}}{E_s}-\frac{\frac{A_{c,i}{f}_{ct}}{1+\sqrt{3.6M_i{\epsilon }_{max,i}}}}{E_s{A}_{s,i}},\hfill \\ M_i=\frac{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, expression (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. [15]), $A_s$ is the cross-section of the longitudinal or transverse steel bars, M is the bond parameter and $A_c$ is the area of concrete bonded to the bar playing in the tension stiffening effect (that usually is considered equal to the rectangular surface around the bar and on a distance not exceeding 7.5$\varphi$ from its center, where $\varphi$ is the diameter of the bar).

In case of prestressed concrete elements, the following two more equations are needed:

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

(14)

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

(15)

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

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 for the prediction of the load-deformation behavior of a prestressed concrete beam is based on the nonlinear system of Equations (3-11), with to 10 equations (notice that Equation (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$).

Equation (13) is based on the concept of forces balance between a general section (or non-cracked section, where both steel and the surrounding concrete contribute) and a cracked section (where only the reinforcement resist the internal forces; 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., to preserve the internal equilibrium of forces, so that as the concrete participation increases, the steel tension decreases) is obtained by applying the following coefficient [17]:

$$\lambda =\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$ represents the boundary of the solvability region for the embedded steel constitutive model 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$.

3. Development and Convergence of the Method

By using the weight matrix function procedure we present a class a 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 free parameters.

On the other hand, since F is a sufficiently differentiable Fréchet function, we can regard $\xi +m\in {\mathbb{R}}^n$ in a neighbourhood of $\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 we can see [18].

Indeed, following the terms introduced by Artidiello et al. in [19], a matrix function $H: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^{\prime \prime }(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.

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

Theorem 1.

Let us $F:D\subseteq {\mathbb{R}}^n\to {\mathbb{R}}^n$ be a sufficiently Fréchet differentiable function defined on an open neighborhood D of $\xi \mathbb{R}$, a zero of F. Let also be $G:{\mathbb{R}}^{n\times n}$ 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 ${x^{\left(k\right)}}_{k\ge 0}$ from (17) converges to ξ with 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$, being $G_0=G\left(I\right)$ and I 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 get

$$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-16C_2^2{C}_4+10C_5{C}_2\hfill \\ & \phantom{=i}-16C_2{C}_4{C}_2-18C_3^2{C}_2+24C_2^2{C}_3{C}_2-16C_4{C}_2^2+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

$$\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,

$$\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+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 $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+4C_5-4b{G}_0{C}_5\hfill \\ & -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-28b{G}_1{C}_2^2{C}_3-4b{G}_2{C}_2^2{C}_3\hfill \\ & +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-4C_4{C}_2+8b{G}_0{C}_4{C}_2\hfill \\ & +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-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-28b{G}_1{C}_2^2{C}_3\hfill \\ & -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-4C_4{C}_2+8b{G}_0{C}_4{C}_2\hfill \\ & +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-3b{G}_2{C}_2{C}_3{C}_2){e^{\left(k\right)}}^5\hfill \\ & +\mathcal{O}\left({e^{\left(k\right)}}^6\right).\hfill \end{array}$$

Now, we denoted $Sc={\left[F^{\prime }\left(x^{\left(k\right)}\right)\right]}^{-1}F\left(z^{\left(k\right)}\right)$. So 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+12b{G}_0{C}_2{C}_3\hfill \\ & +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+24b{C}_2^4{G}_2+4C_5\hfill \\ & -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-50b{G}_0{C}_2^2{C}_3-36b{G}_1{C}_2^2{C}_3\hfill \\ & -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-3b{G}_2{C}_3{C}_2^2-8C_4{C}_2+12b{G}_0{C}_4{C}_2\hfill \\ & +4b{G}_1{C}_4{C}_2+18C_2{C}_3{C}_2-38b{G}_0{C}_2{C}_3-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=iSc+hS$, then

$$\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}$$

Therefore, 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+6G_{0}i{C}_3{C}_2\hfill \\ & -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+M6{e^{\left(k\right)}}^6+\mathcal{O}\left({e^{\left(k\right)}}^7\right)\hfill \end{array}$$

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+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+\mathcal{O}\left({e^{\left(k\right)}}^7\right).\hfill \end{array}$$

With this the proof is finished. □

4. Efficiency Indices

To compare the iterative methods used, we will use the efficiency index $EI={\rho }^{1/d}$, introduced by Ostrowski [20], being $\rho$ the convergence order and d is the amount of functional evaluations, per iteration. Let us also remark that it is necessary to evaluate n scalar functions for each F and $n^2$ for each $F^{\prime }$.

Another index for comparing different iterative schemes was introduced in [21]

In the Figure 3 we compare the computational efficiency index, $CI$, of several methods with Newton’s method.

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

where $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)

Table 1. Comparisons of CE.

Method	CE
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}}$

Table 1. Comparisons of CE.

Method	CE
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 results are represented in 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 3 we can observe that, for $2\le n\le 7$, the best CI index corresponds to Newton method, being O6 the best for $n>\ge 8$. In Figure 3b we can check that for bigger systems, $n\ge 10$, the best CI remains to be O6, our proposed scheme.

5. Numerical Performance

We analyse the performance of the methods described above, to check their efficiency and compare it with other known methods The results from Table 2 to Table 5 correspond to the calculations made with Matlab $R2022b$, by using variable precision arithmetics with 1200 digits of mantissa, on a PC equipped with a $Intel$ Core™ i5-5200U CPU 2.20GHz. In all the tables we show the residuals $\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$ or $\parallel F\left(x^{(k+1)}\right)\parallel$ being lower than ${10}^{-300}$. Moreover, a computational estimation of the order of convergence is obtained by means of ACOC introduced in [22] 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 [23]:

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

(22)

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

In Table 2, 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 convergence accurately, in all cases.

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

5.2. Example

The second example is given by [24]

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

(23)

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

Table 3. 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 3. 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 [24],

$$F_3\left(x_i\right)=arctan\left(x_i\right)+1-2\left(\sum _n^{j=1}{x}_j^2-x_i^2\right),$$

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

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

Let us notice in Table 4 that O6 and FS6 provide a solution satisfying the stopping criterium in less iterations as the rest of schemes. This is the reason why their residuals are not as close to zero as those of other schemes. The e-time of O6 is the second best-one.

Regarding the applied problem described in Section 2, the nonlinear shear model of a reinforced concrete beam, as its underlying data are provided by random values with few digits inside the prescribed ranges by technical standards for structural concrete; moreover, some of this 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$ or $\parallel F\left(x^{(k+1)}\right)\parallel$ being lower than ${10}^{-6}$. The results provided by the new and existing schemes appear in Table 5. The ACOC does not appear in this table, as it yields to unstable data in all cases.

Table 5. 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 5. 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 number of iteration are O6, FS6 and $M_{6,2}(1/2,0)$, all with 3 iterations. Among them, the lowest e-time corresponds to our proposed scheme O6. 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 order six, to solve a constitutive equation of reinforced concrete (6). The new class (O6) is compared to other existing methods with the same order of convergence and also with Newton’s scheme. All the comparison procedures need the same or more, iterations and achieve lower precision results in the same or lower execution time, to achieve the required tolerance. This proves the accuracy, robustness and applicability of the proposed scheme.