Convergence and dynamical study of a new sixth order convergence iterative scheme for solving nonlinear systems

1.   Introduction

From the early ages, real-world and engineering problems have required solving nonlinear equations or systems of nonlinear equations. These solutions usually can only be estimated, as their exact values cannot be calculated. This is the area where iterative methods show their power: to obtain approximations of the solutions of these kinds of problems in the most efficient way.

Therefore, our main goal is to solve the equation $f(x) = 0$, where $f : D\subset \mathbb{R} \rightarrow \mathbb{R}$ is a scalar function, or vectorial $F(x) = 0$, where $F : D\subset \mathbb{R}^n \rightarrow \mathbb{R}^n$, $n > 1$, defined on a convex set $D$. From a seed $x^{(0)}\in \mathbb{R}^n$ ($x_0\in \mathbb{R}$ if $n = 1$), the most employed iterative scheme for estimating the zeroes of these kinds of functions is Newton's scheme, with iterative expression

where $F'(x^{(k)})$ is the Jacobian matrix of $F$ evaluated at $x^{(k)}$.

Some classical results about the extension of scalar one-point iterative methods to systems can be found in texts from Traub ^[29], Ortega and Rheimboldt ^[22] or Ostrowski ^[23]. Also, a good and more recent review can be found in ^[15]. Multipoint methods arrived with the need of improving the efficiency of the schemes without using second-order derivatives. Gerlach in ^[18] and Weerakoon and Fernando in ^[31] used quadrature formulas for this aim, and this task was generalized by Frontini and Sormani in ^[17], Cordero and Torregrosa in ^[12,13], Ozban in ^[24] and Wang et al. in ^[30].

Other techniques were developed to design vectorial methods, such as the pseudocomposition ^[14], still directly related with quadrature formulas. However, Artidiello et al., in ^[2], developed the matrix weight function technique to get higher order schemes. Other authors have also used these and other procedures to construct multidimensional methods, as in ^{[3,4,20,25,26,27,32,34]}, but the amount of this kind of scheme in the literature is still low, compared with their scalar partners. Moreover, it is known that the higher the order of convergence is, the smaller the set of converging initial estimations.

To assure the convergence of Newton's procedure, a seed near to the solutions is necessary, but in general, the modeling of most technical and scientific problems does not hold this condition. To enlarge the set of converging initial estimations, a modification of classical Newton's scheme was designed for the scalar case $f(x) = 0$ in ^[16,33] as

where $\gamma_k$ is defined as a sequence of real numbers depending on previous and actual iterations.

Let us remark that, for fixed values of $\gamma_k$ along the iterative process, the order of convergence is only linear. Nevertheless, if it is calculated at each iteration from the available information, it holds the quadratic convergence, and the sets of converging initial estimations are wider (see ^[5]).

One of the possibilities for $\gamma_k$ is due to Kalitkin and Ermakov (see ^[16]),

where $\|\cdot\|$ denotes any norm. In ^[5], Budzko et al. designed a two-step scheme, and also its multidimensional extension was proposed. In it, the predictor step is a damped Newton type ($\gamma_k = \alpha$) and the corrector is Kalitkin-Ermakov type. In the scalar case,

For systems, coefficient ${\gamma_k}_{eq}$ becomes

with $I$ the identity matrix and $\alpha, b, c$ free parameters.

In this paper, a novel class $(PM_{KE})$ for solving nonlinear systems is proposed and studied. This class has frozen Jacobian and sixth-order convergence. With the aim of widening the regions of starting points that assure the convergence, this family has been designed with the classical Newton's procedure at the first step. Moreover, the second step has a damping factor of Kalitkin-Ermakov type, and the third step is based on a weight function procedure. For this class, the damped coefficient in the second step for systems has the form $\Gamma^{(k)}_{sys}(1, 1, 1) = \Gamma^{(k)}$.

As we construct a whole class of iterative procedures, depending on the weight function to be selected and also on any parameter, we apply techniques from discrete dynamical systems in order to select those members with better stability properties. We understand these improvements in terms of the wideness of the sets of initial guesses able to converge to the roots of nonlinear problems. In ^[11], the authors developed a procedure to arrange this kind of analysis by means of multidimensional real dynamics. From that moment on, different authors have used this technique (or complex dynamics if they return to the scalar case) in order to study the qualitative behavior of the iterative schemes designed. See, for example, ^{[1,6,7,8,9,10,19]}, among others. These techniques are useful in solving nonlinear partial differential equations, as well as stochastic differential equations. In the numerical section, proposed methods are applied to an equation of heat conduction for a non-homogeneous medium.

The information of this paper is organized as follows: In Section 2, we introduce, after some basic definitions, the new three step iterative scheme, and its convergence study is presented in Section 3. Next, the dynamical analysis of the proposed method is realized in Section 4. Then, some numerical experiments are shown, and some conclusions are stated, in Sections 5 and 6, respectively.

2.   Notation and concepts

For the sake of completeness, we first introduce some main definitions about iterative processes in numerical analysis. After that, we define some concepts and properties of discrete dynamical systems, which are necessary to follow the results obtained in Section 4.

Let us consider a sequence $\{x^{(k)}\}_{k\geq0}$ in ${\mathbb{R}}^{n}$ converging to the root of $F(x) = 0$. Then, convergence has order $p\geq1$, if there exist $N > 0$ ($0 < N < 1$ with $p = 1$) and $k_0$ satisfying

or

Assuming three consecutive iterates $x^{(k-1)}, \; x^{(k)}, \; x^{(k+1)}$ close to $\xi$, $\rho$ was introduced in ^[12] as the Approximated Computational Order of Convergence (ACOC),

On the other hand, as $F$ is a sufficiently Fréchet differentiable function, we can consider $\xi + h \in {\mathbb{R}}^n$ in a neighborhood of $\xi$. By using Taylor expansions and $F'(\xi)$ being non-singular,

In this expression, $C_q = \frac{1}{q!}[F'(\xi)]^{-1} F^{(q)}(\xi)$, $q\geq2$. Let us remark that $C_qh^q\in {\mathbb{R}}^n$, since $[F'(\xi)]^{-1}\in \mathcal{L}({\mathbb{R}}^n)$, and $F^{(q)}(\xi) \in \mathcal{L}({\mathbb{R}}^n\times\cdots\times {\mathbb{R}}^n, \, {\mathbb{R}}^n)$. Then, we can express $F'$ as

where $qC_qh^{q-1} \in \mathcal{L}({\mathbb{R}}^n)$.

Moreover, following the notation defined by Artidello et al. in ^[2], if $X = \mathbb{R}^{n \times n}$ is the Banach space of real $n \times n$ matrices, a matrix function $H:X \rightarrow X$ can be defined such that its Fréchet derivatives satisfy

(a) $H'(u)(v) = H_1uv, \; \mathrm{being}\; H':X\rightarrow \mathcal{L}(X), \ H_1 \in \mathbb{R}$,

(b) $H''(u, v)(v) = H_2uvw, \; \mathrm{being}\; H'':X \times X\rightarrow \mathcal{L}(X)\ H_2 \in \mathbb{R}$.

Let $G$ be a vectorial rational operator related to an iterative scheme applied on a polynomial function $p : D \subseteq \mathbb{R}^n \rightarrow \mathbb{R}^n$, and the orbit of a point $x\in \mathbb{R}^n$ is defined as the successive images $G(x^{0})$ in the form $\{x, G(x), \dots, G^m(x)\}$.

If $x^*$ is a such point that $G(x^*) = x^*$, then it is denominated a fixed point of $G$. This point is attracting if all the eigenvalues $\lambda_j$ of $G'(x^*)$ satisfy the condition $|\lambda_j| < 1$, it is unstable if one eigenvalue satisfies the condition $|\lambda_{j0}| > 1$, and it is repelling if all $\lambda_j$ satisfy $|\lambda_j| > 1$. On the other hand, the basin of attraction $\mathcal{B}(x^{*})$ of a certain $x^*$ is defined as the set of pre-images of any order such that

If $x_c$ is such a point that $G(x_c) = 0$, then it is called a critical point. Those critical points different from the roots of $F(x) = 0$ are called free critical points. In each basin of attraction of a periodic point, there exists one critical point, so the study of those points and their orbits is crucial.

All these notations are widely used in the demonstration of the convergence and dynamical study in the following sections.

3.   Convergence analysis of the proposed class $PM_{KE}$

Let us consider the three step iterative method $PM_{KE}$:

where $\Gamma^{(k)} = \left[I-\left(I-[F'(x^{(k)})]^{-1}[y^{(k)}, x^{(k)}; F] \right)^2\right]^{-1}$ and $\tau^{(k)} = I+ [F'(x^{(k)})]^{-1}[z^{(k)}, y^{(k)}; F] \:\Gamma^{(k)}$. In the following result, we state the conditions to be satisfied to get fifth- and sixth-order convergence of class (3.1).

Theorem 3.1. Let us consider a sufficiently Fréchet differentiable function $F : D \subseteq \mathbb{R}^n \rightarrow \mathbb{R}^n$ defined in an open neighborhood $D$ of its zero $\xi \in \mathbb{R}$, and $H:\mathbb{R}^{n \times n} \rightarrow \mathbb{R}^{n \times n}$ is a sufficiently differentiable matrix function. Let us also assume that $F'(\xi)$ is non-singular, and $x^{(0)}$ is a seed near enough to $\xi$. Then, the sequence $\{x^{(k)} \}_{k \geq 0}$ obtained from (3.1) converges to $\xi$ with order 5 if $H_0 = (2+H_2)I, H_1 = -(1+2 H_2)$, and $|H_2| < \infty$, where $I$ is the identity matrix, and $H_0 = H(I)$. In this case, the error equation is

Moreover, if $H_2 = 1$, the order of convergence is six, and the error equation is

where $C_q = \frac{1}{q!} [F'(\xi)]^{-1} F^{(q)}(\xi)$, $q = 2, 3, ...$

Proof. By means of Taylor series of $F(x)$ around $\xi$, we get the expansion of $F(x^{(k)})$, $F'(x^{(k)})$, $F''(x^{(k)})$ and $F'''(x^{(k)})$ around $\xi$. Indeed,

where $X_2 = -2C_2$, $X_3 = -3C_3+4C^2_2$, $X_4 = -4C_4+6C_2C_3+6C_3C_2-8C^3_2$, $X_5 = -5C_5+8C_2C_4-12C^2_2C_3+9C^2_3+8C_4C_2-12C_2C_3C_2+16C^4_2-12C_3C^2_2$, and

These values of $X_i$, for $i = 2, 3, \ldots, 6$, have been obtained by assuming that $[F'(x^{(k)})]^{-1}[F'(x^{(k)})] = I$ is satisfied and solving the subsequent linear system. Then,

and

Therefore,

Now, we can express the last product of the second step as

Using the Hermite-Genocchi formula (see ^[22]) for the operator of first divided difference $[x+h, x; F] = \int_{0}^1 F'(x+t h) dt = F'(x)+\frac{1}{2}F''(x)h+\frac{1}{6}F'''(x)h^2+O(h^3)$, $\forall (x, h)\in {\mathbb{R}}^n \times {\mathbb{R}}^n$, we can obtain the operator $[y, x; F]$ related to step two of the iterative scheme (3.1):

and

The Taylor expansion around the solution of the term $B = I+\left(I-[F'(x^{(k)})]^{-1}[x^{(k)}, y^{(k)}; F]\: \right)^2$ in the second step of (3.1) can be expressed as

Let us now find the inverse of $B$ that can be defined as $B^{-1} = \Gamma^{(k)} = I+K_1e^{(k)}+K_2{e^{(k)}}^{2}+K_3{e^{(k)}}^{3}+K_4{e^{(k)}}^{4}+O({e^{(k)}}^{5})$. Fixing the condition $B{B}^{-1} = {B}^{-1}B = I$, we can find the coefficients:

Therefore,

and with the help of expressions (3.3), (3.5) and (3.8), we obtain the error in the second step:

To get the error equation, we obtain first the first order divided difference $[z^{(k)}, y^{(k)}; F]$, knowing that

and

Using (3.10) and the Hermite-Genocchi formula, we get

By using (3.2), (3.8) and (3.12), the variable $\tau^{(k)}$ can be expanded as

Then, employing (3.9), (3.11)–(3.13) and the Taylor expansion of the weight function $H(\tau^{(k)}) = H_0+H_1(\tau^{(k)}-I)+H_2 {(\tau^{(k)}-I)}^2+O({(\tau^{(k)}-I)}^3)$, the error equation takes the form

Fixing $H_0 = (2+H_2)I$ and $H_1 = -(1+2H_2)$, the error equation becomes

Finally, making $H_2 = 1$, we obtain

and the proof is complete.

3.1.   Specific iterative subfamilies of $PM_{KE}$

Once the conditions that the weight function must meet have been obtained, there are several ways to define the matrix weight function $H$, and each of them generates a different iterative family.

Family 1. The weight function defined by the polynomial

where $\beta\in\mathcal{R}$ satisfies the hypothesis of Theorem 3.1 and $\tau^{(k)} = I+ [F'(x^{(k)})]^{-1}[z^{(k)}, y^{(k)}; F]\left[I+\left(I-[F'(x^{(k)})]^{-1}[y^{(k)}, \ x^{(k)}; F]\right)^2\right]^{-1}$, generates a new parametric class of fifth order for $\beta\neq1$ and a simple scheme of sixth order for $\beta = 1$,

This family is denoted by ${PM(\beta)_{KE}}_p$ depending on the $\beta$ parameter.

Family 2. The weight function defined by

where $\beta\in\mathcal{R}$ also satisfies the hypothesis of Theorem 3.1, generates a new class of fifth and sixth order with $\beta\neq1$ and $\beta = 1$, respectively.

with $\tau^{(k)} = I+ [F'(x^{(k)})]^{-1}[z^{(k)}, y^{(k)}; F]\left[I+\left(I-[F'(x^{(k)})]^{-1}[y^{(k)}, \ x^{(k)}; F]\right)^2\right]^{-1}$. In the following we will use the family 1 $({PM(\beta)_{KE}}_p)$ for all the numerical studies. The algorithm of this class is presented.

Algorithm:

● Step 1: Knowning the nonlinear system to be solved, defined by function $F$, and its related Jacobian matrix $F'$, an initial estimation $x^{(0)}$, a small constant $tol$ as the accuracy, let $k = 0$, and start the iteration process.

● Step 2: Apply iterative expression 3.17 to get the next iterate, $x^{(1)}$.

● Step 3: If $\|F(x^{(1)})\| < tol$, we stop and output $x^{(1)}$ as the solution. Otherwise, set $k = k+1$, go to Step 2 and continue the iteration process.

4.   Dynamical analysis of the family ${PM(\beta)_{KE}}_p$

In this section, we study the performance of the vectorial rational function obtained by applying iterative class ${PM(\beta)_{KE}}_p$ on $p (x) = \{x ^ 2_1-1, \, x ^ 2_2-1 \}$ for $x = (x_1, x_2)\in\mathbb{R}^2$. Because of the dependence of this family on the beta parameter, it is important to study the existence of fixed points different from the roots of the system $p(x) = 0$, mainly those with attractive character and other attractive structures that can be interesting from the dynamical point of view. Through this analysis, those members whose only behavior is the convergence to the roots will be derived.

So, by applying ${PM(\beta)_{KE}}_p$ on $p(x)$, we get its rational vectorial operator $G(x, \beta) = \{g_1(x, \beta), g_2(x, \beta)\}$, whose $j$-th coordinate is

Taking into account the last result, it is possible to formulate the following theorem about the stability of fixed points of $G(x, \beta)$.

Theorem 4.1. The duples $(1, 1), \; (1, -1), \; (-1, 1), \; (-1, -1)$ are superattractive fixed points of the rational function $G(x, beta)$ which are zeroes of $p(x)$. In addition, denoting by $\mathcal{B}$ the set of strange fixed points of $G(x, \beta)$, it is defined by $(l_i, l_j)$, $(\pm 1, l_j)$ and $(l_i, \pm 1)$ for $i, j\leq26$, whose elements different from $\pm1$ are real zeroes of

Hence, the amount of fixed points in $\mathcal{B}$ depends on $\beta$:

i) If $\beta\in (-\infty, -61.956\, ]$, $\mathcal{B}$ is composed by thirty-six repulsive and twenty-four saddle strange fixed points.

ii) If $\beta \in (-61.956, -5.4247)$, $\mathcal{B}$ has the same composition as the previous item.

iii) If $\beta \in (-5.4247, -4.8788)$, $\mathcal{B}$ is composed by sixteen repulsive and sixteen saddle strange fixed points.

iv) If $\beta \in [-4.8788, -2)$, the composition of $\mathcal{B}$ is the same as the previous item.

v) For $\beta \in [-2, -0.9328)$, there exist sixty strange fixed points, whose characters depend on $\beta$. Because of their stability, two different situations can be meet:

a) If $\beta \in [-2, -0.9578)$, $\mathcal{B}$ has the same composition as the first item.

b) When $\beta \in (-0.9578, -0.9329)$, $\mathcal{B}$ has twelve attractor, sixteen repulsive and thirty-two saddle strange fixed points.

vii) If $\beta \in (-0.9329, +\infty)$, $\mathcal{B}$ has four repulsive and eight saddle strange fixed points.

Proof. $G(x, \beta)$ is symmetric with respect to its component functions. Considering that the fixed points are solutions of the equation $g_j(x, \beta) = x_j$, $j = 1, 2$, we get

where

Considering (4.2), it is observed that $x_j = \pm1$ are solutions of this equation, and therefore, $(1, 1)$, $(1, -1)$, $(-1, 1)$, $(-1, -1)$ are fixed points of the rational operator $G(x, \beta)$ and at the same time roots of the polynomial $p(x)$. In order to study the stability of the rational operator, the diagonal matrix form of the associated Jacobian $G'(x, beta)$ is computed:

where

and

It is straightforward that $J_j(\pm1, \beta) = 0$, $j = 1, 2$, for any $\beta$, and therefore the zeroes of $p(x)$ are superattracting fixed points, given by the fact that the eigenvalues of $G'((\pm 1, \pm 1), \beta)$ are null.

Fixed points can be calculated through $l(t, \beta)$ by means of $s = t^2$. A $13$-th degree polynomial is then obtained,

The operator $G(x, beta)$ has real fixed points that are symbolized as $L_i$, with $L_i$ any real and positive root of $L(s)$. Therefore, the amount of components of $\mathcal{B}$ is related to the number of real and positive roots of $L(s)$, as well as their combinations with $\pm1$. Then,

$\rm i)$ Three real and positive roots $L_{1}$, $L_{2}$, $L_{3}$ are found if $\beta\in (-\infty, -61.956\, ]$, and the roots $\{+\sqrt{L_{1}}, +\sqrt{L_{2}}, +\sqrt{L_{3}}, -\sqrt{L_{1}}, -\sqrt{L_{2}}, -\sqrt{L_{3}}\}$ are denoted by $l_i$ for $i\in\{1, 2, 3, 4, 5, 6\}$. The set of all strange fixed points is obtained combining in pairs $l_i$, $i = 1, 2, \ldots, 6$ with themselves and with $1$ and $-1$. In the first case, we have that all combined pairs $(l_i, l_k)$ are 36, in the second the kind of $(l_i, \pm 1), (\pm 1, l_i)$ pair are $12+12 = 24$. The information about stability of these fixed points of $G(x, \beta)$ can be inferred from the analysis of the absolute value of the eigenvalues of matrix $G'((l_i, l_k), (\beta))$, $i, k\in\{1, 2, 3, 4, 5, 6\}$; these functions of $\beta$ are called stability functions of the respective fixed points. Due to the nature of the polynomial system, the eigenvalues $\lambda((l_{i}, l_k), (\beta)) = J_1((l_i, l_k), (\beta)) = J_2((l_i, l_k), (\beta))$ for $i, k\in\{1, 2, 3, 4, 5, 6\}$; however, if any of the components of the fixed point is $\pm1$, the corresponding eigenvalue is always null.

Figure 1 shows the values of $|\lambda|$ associated to the Jacobian matrix and evaluated at $l_i$ for $i\in\{1, 2, 3, 4, 5, 6\}$. For all of them the corresponding absolute values of $\lambda$ are greater than one, so the behavior of the strange fixed points $(l_i, l_k)$ is repulsive for the current interval. Indeed, fixed points $(l_i, \pm 1), (\pm 1, l_i)$, for $i\in\{1, 2, 3, 4, 5, 6\}$, are saddle points, as one of the eigenvalues is zero, and the another one is greater than one.

$\rm ii)$ The amount of fixed points for $\beta \in (-61.956, -5.4247)$ is the same as in the previous case, due to the existence of three real and positive roots. The information about qualitative behavior of fixed points can be extracted from the stability functions represented in Figure 2. For all of them, the corresponding absolute values of $\lambda$ are greater than one, and therefore the pair $(l_i, \pm 1), (\pm 1, l_i)$ are saddle points for $i\in\{1, 2, 3, 4, 5, 6\}$.

$\rm iii)$ If $\beta \in (-5.4247, -4.8788)$ the polynomial $L(s)$ has two real positive roots $L_1$, $L_2$, and the elements $\{+\sqrt{L_{1}}, +\sqrt{L_{2}}, -\sqrt{L_{1}}, -\sqrt{L_{2}}\}$ are denoted by $l_i$ for $i\in\{1, 2, 3, 4\}$. The analysis of stability functions associated to the fixed points in Figure 3, shows that the pairs $(l_i, l_k)$ for $i, k\in\{1, 2, 3, 4\}$ are repulsive, because their corresponding $|\lambda|$ are greater than one, so $(l_i, \pm 1)$ and $(\pm 1, l_i)$ are saddle fixed points.

$\rm iv)$ For $\beta \in [-4.8788, -2)$, the same structure and the same reasoning as in the previous case take place; see Figure 4.

$\rm v)$ As the first item, in this case $L(s)$ has three real positive roots, and the elements $l_i$ are denoted in the same way for $i\in\{1, 2, 3, 4, 5, 6\}$. For some values of $\beta$, there are attractive fixed points, and such values can be calculated by solving the equation $|lambda_{j}((l_{i}, l_k), (\beta))| = 1$, for $i, k\in\{1, 2, 3, 4, 5, 6\}$ and $j = 1, 2$. As a result we have found that the values of $\beta_1\approx-0.9578$ and $\beta_2\approx-0.9328$ satisfy the above equation, related to the strange fixed points $(l_{i}, l_k)$ with $i, k$, Figure 5c, d.

a) For $\beta \in [-2, \beta_1)$ all thirty-six pairs $(l_i, l_j)$ are repulsive fixed points, as absolute values of $|\lambda_j((l_{i}, l_k), (\beta))| > 1$ for $i, k\in\{1, 2, 3, 4, 5, 6\}$ and $j = 1, 2$. On the other hand, pairs like $(l_i, \pm 1)$, $(\pm 1, l_i)$ are saddle fixed points, and their quantity is twenty-four for $i\in\{1, 2, 3, 4, 5, 6\}$. The stability of any of these fixed points can be inferred from Figure 5.

b) There are not hyperbolic points for the values $\beta = \beta_1, \beta_2$. When $\beta\in (\beta_1, \beta_2)$ the stability of fixed points $(l_{i}, \pm1)$, $(\pm1, l_k)$ and $(l_i, l_k)$ for $i, k\in\{2, 5\}$ correspond to a Jacobian matrix whose eigenvalues take values $|\lambda_{j}((l_{i}, \pm 1), (\beta))| < 1$, $|\lambda_{j}((l_{i}, \pm 1), (\beta))| = 0$ and $|\lambda_{j}((l_{i}, l_k), (\beta))| < 1$, $|\lambda_{j}((l_{i}, l_k), (\beta))| = 0$ for $j = \{1, 2\}$, so they are attracting, and their quantity is twelve. On the other hand, all forty-eight pairs $(l_2, l_i)$, $(l_i, l_2)$, $(l_5, l_i)$, $(l_i, l_5)$, $(\pm1, l_i)$ and $(l_i, \pm1)$ are saddle fixed points due to the fact that $|\lambda_{j}((l_{i}, l_k), (\beta))| > 1$, for $i\in\{1, 3, 4, 6\}$, Figure 5.

vi) Finally, when $\beta \in (-0.9328, +\infty)$ the polynomial $L(s)$ has one real positive root $L_1$, and the elements $\{+\sqrt{L_{1}}, -\sqrt{L_{1}}\}$ are denoted by $l_i$ for $i\in\{1, 4\}$. The four fixed points $(l_i, l_k)$ for $i, k\in\{1, 4\}$ and $j = 1, 2$ are repulsive, as absolute values of $|\lambda_j((l_{i}, l_k), (\beta))| > 1$. On the other hand, pairs such as $(l_1, \pm1)$, $(\pm1, l_1)$, $(l_4, \pm1)$ (eight in total) are saddle fixed points due to the fact that $\lambda_j((\pm1, \pm1), \beta) = 0$ and $|\lambda_j((l_{i}, l_k), \beta)| > 1$ (see Figure 6).

The stability of strange fixed points is already known, and it has been shown that attractive fixed points only appear for $\beta \ at (-0.9578, -0.9329)$. However, in addition to the latter, there are other attractive structures whose values of the $\beta$ parameter should also be avoided to achieve stable performance of the iterative methods. These structures must be sought in basins of attraction in which the free critical points, those different from the roots of $p(x) = 0$, are also found. Related to the latter, we can make the following statement about the location of the critical points.

Theorem 4.2. Let $\mathcal{K}$ denote the set of free real critical points of $G(x, \beta)$ related to iterative family 2 and let the values $\{\alpha_1 = -0.46806, \alpha_2 = 6.28814, \alpha_3 = 66.1728\}$ for the real roots collection of a 23th-degree polynomial.

Then, $\mathcal{K}$ is the collection of the duples $(q_i, q_j)$, $(q_i, \pm 1)$ and $(\pm 1, q_j)$ for $i, j\leq10$ whose entries different from $\pm1$ are the real zeroes of polynomial

Then, $\mathcal{K}$ is made up as follows:

i) If $\beta \in (-\infty, -2]$, $\mathcal{K}$ contains 12 free critical points.

ii) If $\beta \in (-2, -0.5]$, $\mathcal{K}$ contains 32 free critical points.

iii) If $\beta \in (-0.5, \alpha_1]$, $\mathcal{K}$ contains 100 free critical points.

iv) If $\beta \in (\alpha_1, 1.40576)$, $\mathcal{K}$ contains 32 free critical points.

v) If $\beta \in (1.40576, \alpha_2)$, $\beta \in [\alpha_2, \alpha_3)$ or $\beta \in (\alpha_3, +\infty)$, $\mathcal{K}$ contains 12 free critical points.

Proof. Is known from Theorem 4.1 that the eigenvalues of Jacobian matrix $G'(x, \beta)$ are $\lambda_j(x, \beta) = J_j(x, \beta)$, for $j\in\{1, 2\}$, that is,

where

On the other hand, the critical points are the solutions of equation $\lambda_j(x, \alpha) = 0$ for $j\in\{1, 2\}$ that are found by using $t^2 = s$ on $r(t)$. Its roots (real and positive, symbolized as $Q_j$) will be the entries of the free critical points, as $q_k = \pm \sqrt{Q_j}$. The quantity of positive roots depends on the value of parameter $\beta$, forcing $Q_j$ to be real:

(i) When $\beta \in (-\infty, -2]$, $r(s)$ has one real positive root, denoted by $Q_1$. Then, $\mathcal{K}$ is composed by $(q_1, q_1), \; (q_1, q_5), \; (q_5, q_1)$ and $(q_5, q_5)$, where $q_{1, 5} = \pm \sqrt{Q_1}$, and also by eight critical points of the kinds $(q_j, \pm1)$ and $(\pm1, q_j)$, $j\in\{1, 5\}$.

(ii) If $\beta \in (-2, -0.5]$, $r(s)$ has only two positive real roots $Q_1$ and $Q_2$. Therefore, $\mathcal{K} = \{(q_{i}, q_{j}), (q_i, \pm 1), (\pm 1, q_{j}) \ : \ i, j\in\{ 1, 2, 5, 6 \} \}$, holding 32 different free critical points.

(iii) For $\beta \in (-0.5, \alpha_1]$, there exist four real positive roots of polynomial $r(s)$: $Q_1$, $Q_2$, $Q_3$ and $Q_4$. So, $\mathcal{K}$ has the elements of kind $(q_{i}, q_{j})$, with $i, j\in \{1, 2, \ldots, 8\}$, where $q_{1, 5} = \pm \sqrt{Q_1}$, $q_{2, 6} = \pm \sqrt{Q_2}$, $q_{3, 7} = \pm \sqrt{Q_3}$ and $q_{4, 8} = \pm \sqrt{Q_4}$. In addition, mixed points $(\pm 1, q_{j})$ and $(q_i, \pm 1)$ where $i, j\in \{1, 2, \ldots, 8\}$ also belong to $\mathcal{K}$. The total quantity of critical points is 100.

iv) When $\beta \in (\alpha_1, 1.40576)$, $r(s)$ has two real positive roots, and the $\mathcal{K}$ composition is the same as in (ii).

v) If $\beta \in (1.40576, \alpha_2)$, $\beta \in [\alpha_2, \alpha_3)$ or $\beta \in (\alpha_3, +\infty)$, $r(s)$ has one real positive root (different in each case), and the $\mathcal{K}$ composition is the same as in (i).

4.1.   Study of pathological behavior

Theorems 4.1 and 4.2 give us information about stability of the iterative class acting on a certain test polynomial. In general, elements of the strange fixed and critical point collections should be avoided as initial estimates.

However, from the dynamic point of view, the analysis of the long-term performance of free critical points can give us the key about other strange attractive structures, periodic orbits or even chaotic behavior.

To take a first glance to the study of the critical points we are going to use the tool named parameter line, first used in ^[21].

This tool consists of a mesh of $500 \times 500$ points in a specific interval of the $\beta$ parameter. Each parameter-dependent free critical point used as initial estimate is colored red, if after 200 iterations it converges to some root of the polynomial $p(x)$, and black in any other case. Each of these points is thickened by a multiplication by $[0, 1]$. The error tolerance is $10^{-3}$. The red color is more clear when the number of iterations needed to converge is lower.

Figure 7 shows the limit performance of free critical points for $\beta \in (-2, -0.5]$. In Figure 7a, all the free critical points $(q_i, q_j)$ with $i, j\in\{1, 5\}$ converge to any of the roots of $p(x)$; moreover, in Figure 7b, unstable behavior is discovered for critical points $(q_i, q_j)$ with $i, j\in \{2, 6\}$, $(q_1, q_m)$ for $m\in \{2, 6\}$ and $(q_2, q_n)$ for $n\in \{1, 5\}$ around the values $\beta\approx -1.75$ and $\beta\approx-0.95$.

We continue the study of the unstable performance in the two dark zones using bifurcation diagrams. They are calculated using each free critical point of the rational operator as a seed, for $\beta\in (-2, 0.5\, ]$, in a mesh of $3000$ subintervals. In this interval they are real, and we observe their performance in the last $100$ of $1000$ iterations. In Figure 8, detailed bifurcation diagrams of critical point $(q_2, q_6)$ are shown. They are related with the dark zones of no convergence in the parameter lines.

Because of polynomial $p(x)$ has separate variables, the coordinate functions of the rational operator are symmetrical. Taking it into account, in the bifurcation diagrams the $x_1$ and $x_2$ dimensions are plotted in blue and red color, respectively. For the details, only positive coordinates were chosen.

Observing the details of bifurcation diagrams given by Figure 8b, e, we can find patterns of bifurcations, such as period doubling cascade, periodic orbits and chaotic behaviors for $\beta \approx-1.76, -1.762$ and $-1.765$ in the first case and for $\beta \approx-0.96, -0.97$ and $-0.986$ in the second case, respectively. In the following we will find and show some of these pathological attractive structures previously observed in the bifurcations diagrams.

To take a full view of dynamic behavior of the rational operator $G(x, \beta)$ we plotted a dynamical plane, and these graphics are constructed by iterating each point of a mesh with a 0.01 step and painting according to the root to which they converge, with an error tolerance lower than $10^{-3}$ or with a maximum of 50 iterations as a stopping criterion. The brighter the color is the smaller the amount of iteration needed to converge. In case of non-convergence after the number of iterations is completed, the point is colored black. The roots of the polynomial are represented with circles.

Related with the first dark zone on parameter line (see Figure 7b), in Figure 9a, we can see that the amount of connected components in the basins of attraction of the roots for $\beta = -1.7665$ is infinite. Green and white diamonds represent the repulsive and saddle strange fixed points, respectively. All of them lie in the Julia set. On the other hand, free critical points are represented with yellow squares. Some of them lie in the basin of attraction of the roots of $p(x)$, and the rest are located in the black zone. Observing the bifurcation diagram in Figure 8b and knowing the position of critical point ($q_2, q_6$), we can explore the regions around points $(\mp 0.48, \pm 0.48)$ to find attractive elements. In Figure 9a, we can observe an attractive structure (see small white region around the points $(\mp 0.48, \pm 0.48)$). Any initial estimation in those zones tends to densely fill these same symmetrical regions. The details of these small regions can be observed in Figure 9b.

On the other hand, for the value of $\beta = -1.7621$ (see the bifurcation diagram, Figure 8b), we find periodic orbits around critical point $(q_2, q_6)$. One of them, with period 8, is shown in Figure 9c, and its elements are

Its details can be seen in Figure 9d.

Finally, by studying the second dark zone of the parameter line in Figure 7b, it is important to note in Figure 8d, e the same color of the bifurcation patterns. This means that the $x_1$ and $x_2$ coordinates maintain the same sign in any attractive structure, that is, each of those attractive elements must be confined in small areas where the coordinates sign of its points do not change. In Figure 10a, for $\beta = -0.986$, this behavior can be seen in a small area around critical point $(q_2, q_6)$ and the coordinates $(\mp 0.48, \mp 0.48)$ where every iteration of initial estimates tend densely to fill this area. Details are shown in Figure 10b. Meanwhile, observing the bifurcation diagram in Figure 8e for $\beta = -0.97$ we find periodic orbits of period four, and one example can be seen in Figure 10c, d whose components are ${(0.4589, -0.4589), (0.4891, -0.4891), (0.4543, -0.4543), (0.5010, -0.5010)}$.

5.   Numerical performance

In what follows, we show the behavior of ${PM(\beta)_{KE}}_p$ (described by (3.15)) on some real-life examples and in comparison with the Newton method. The $\beta$ values used ($\beta = 1, 2$) provides convergence of order 6 and 5, respectively, in accordance with the obtained results of Theorem 3.1. To develop the numerical test, we use MATLAB software with variable precision arithmetic and 2000 digits of mantissa. The experiments have been performed on a computer with Intel(R) Core(TM) i7-1065G7 CPU and 16 GB of RAM. The stopping criterion used is $\|F(x^{(k+1)})\| < 10^{-200}$ or $\|x^{(k+1)}-x^{(k)}\| < 10^{-200}$. The first order divided difference operator used has the elements given by the expression (see ^[22])

being $f_i$, $i = 1, 2, \ldots, n$, the coordinate functions of $F$. For each system, one table of results is displayed, where $x^{(0)}$ is the initial approximation, $k$ is the number of iterations needed to converge ("nc" appears in the table if the method does not converge), the value of residuals $\|F(x^{(k+1)})\|$ or $\|x^{(k+1)}-x^{(k)}\|$ at the last iterate ("-" if there is no convergence), and $ACOC$ ^[12] is

In relation to the latter, it is necessary to determine when the root($\|F(x^{(k+1)})\| < 10^{-200}$ is achieved) is reached or if it is only a very slow convergence with a no significant difference between the two last iterates ($\|x^{(k+1)}-x^{(k)}\| < 10^{-200}$ but $\|F(x^{(k+1)})\| > 10^{-200}$), or both criteria are satisfied. The nonlinear test systems used are the following:

$\bf{Example\, 1.}$ In this example, with the use of ${PM(1)_{KE}}_p$, we will study the equation of heat conduction ^[28] for a non-homogeneous medium :$\frac{\partial}{\partial x}\left(k(u)\frac{\partial u}{\partial x}\right)+ F(x, t) = c\rho\frac{\partial u}{\partial t}$, with $F(x, t)$ the thermal source density and $k$, $c$ and $\rho$ the material quantitative properties of the thermal conductivity coefficient, specific heat capacity and medium density, respectively.

For the sake of simplicity we are taking $c$, $\rho$ as constants and the thermal conductivity coefficient as a linear function of temperature $k(u) = k_0+k_1u$ with $k_0 = 1$. The thermal sources have been distributed equispaced throughout the space-time grid, and their signs have been chosen randomly. After a few transformations and adding some boundary and initial conditions, we state the problem to be studied:

where $f = f(x, t)$ is obtained from F(x, t) (the thermal source density) after the transformations leading to the definition of the resulting boundary problem. By applying an implicit finite differences scheme, we transform the equation of heat conduction in a family of nonlinear systems to be solved. The estimated solution for each instant $t_k$ is obtained from the approximation calculated at $t_{k-1}$. The space-time mesh is constructed taking the spatial step $h = \frac{8}{nx}$ and the timing step $l = \frac{T_{max}}{nt}$, with $nx$, $nt$ and $T_{max}$ the amount of subintervals for variables $x$, $t$ and final instant of the numerical study, respectively.

Making $x_0 = -4$ and $x_f = 4$, we get the mesh domain as $[x_0, x_f] \times [0, T_{max}]$ with $(nx+1)\times(nt+1)$ equispaced nodes $(x_i, t_k)$. Denoting by $u_{i, k}$ the estimation of the solution at the $(x_i, t_k)$ and using the estimations of the derivatives

we construct the nonlinear system with $\alpha = 1$,

The results for different discretization and thermal sources are shown in Table 1. The mean of the number of iterations employed for solving the nonlinear systems for each instant is denominated as "iter" and $\|F(x^{(k+1)})\|$ is the norm of the solution of the last system. As the number of iterations employed for solving each system is very low, the ACOC cannot be calculated.

In Figure 11 are shown the approximated solutions of the head equation for different thermal sources. See the influence of $f(x, t)$ in Figure 11b–d, compared with Figure 11a.

Table 1 and Figure 11 show that the scheme is robust enough to reproduce with effectiveness the behavior of the model for different values of thermal sources density.

$\bf{Example\, 2.}$ This system has variable size and has the solution $\xi\approx(0.5149, ..., 0.5149)^T$. It is defined by

with $n = 5$, and the initial estimation is $x^{(0)} = (0.75, ..., 0.75)^T$. The results appear in Table 2.

In there, it can be observed that the ${PM(\beta)_{KE}}_p$ converges to the solution in fewer iterations than Newton's method, as we might expect. In the case of ${PM(2)_{KE}}_p$ the convergence occurs with null (for the fixed precision of the machine at 2000 digits) or almost null residual $\|F(x^{(k+1)})\|$ for ${PM(2)_{KE}}_p$ and the parameter values involved. In addition, the ACOC is in good accordance with the theoretical one, as we can expect due to Theorem 3.1.

$\bf{Example\, 3.}$ Next system of size $n = 10$ is defined as

whose solution is $\xi\approx(0.1758, ..., 0.1758)^T$, and the seed used is $x^{(0)} = (1.0, ..., 1.0)^T$. Results are given in Table 3. The ${PM(\beta)_{KE}}_p$ family converges with null residuals $\|F(x^{(k+1)})\|$ in fewer iterations than Newton's method, and the ACOC is in total agreement with the theoretical one.

$\bf{Example\, 4.}$ Finally, we test the ${PM(\beta)_{KE}}_p$ family with a nonlinear system of variable size described as

with $n = 20$. The initial estimation $x^{(0)} = (0.25, ..., 0.25)^T$, and the solution is $\xi = (0.0, ..., 0.0)^T$. The results are given in Table 4. The results are similar to those obtained in Example 2, with fewer iterations than Newton's method, null (for the fixed precision of the machine at 2000 digits) or almost null residual $\|F(x^{(k+1)})\|$ and the computational approximation order of convergence in good accordance with the theoretical one.

6.   Conclusions

In this manuscript, a class of iterative schemes with sixth-order convergence has been designed. Its construction has been made using the Ermakov-Kalitkin approach, which improves the wideness of the set of converging initial guesses. The analysis of the stability of the resulting class has proven the inclusion of very stable members, which have been used in the numerical section to solve real-life and academical problems.

Acknowledgments

This research was partially supported by Grant PGC2018-095896-B-C22 funded by MCIN/AEI/31000.13039/ "ERDF A way to making Europe", European Union.

Conflict of interest

There is no conflict of interest regarding publication of this manuscript.