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

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

In [1], the authors proposed this stress-strain relationship for concrete cracked in tension: 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$: 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: 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

On the other hand, in CFT models the equilibrium between the external loads and the internal forces is governed by the following equations: 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): 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: 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]: 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$.

By using the weight matrix function procedure we present a class a three-step iterative methods with the following iterative expression: 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, 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, 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

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