A new family of fourth-order Ostrowski-type iterative methods for solving nonlinear systems

1.   Introduction

Solving nonlinear systems of equations is a classic problem in engineering and experimental science. Many practical problems can be transformed into integral equations, partial derivative equations, matrix equations, and other methods for nonlinear systems. The iterative method is the most common method for solving nonlinear systems.

In order to solve the problem of finding a solution ${{x}^{*}}\in {{\mathbb{R}}^{n}}$ to a nonlinear system $F\left(x \right) = 0$, where $F:E\subseteq {{\mathbb{R}}^{n}}\to {{\mathbb{R}}^{n}}$ and $F$ is a differentiable function in an open convex set $E$ of a sufficient Fréchet, an approximation to the exact solution ${{x}^{*}}$ can be obtained through the iterative method. The iterative method can be described as

Starting from the initial estimate ${{x}^{\left(0 \right)}}\in {{\mathbb{R}}^{n}}$, as an approximation to the exact solution ${{x}^{*}}$, an iteration sequence $\left\{ {{x}^{\left(j \right)}} \right\}$ of approximations to ${{x}^{*}}$ is obtained by the iterative method (1.1). In the last decades, the most commonly used form of iterative is Newton's method ^[1] with second-order of convergence. In order to improve its convergence order, King ^[2] introduced a fourth-order convergence iterative method for solving nonlinear systems. Starting from a fourth-order iterative family for solving nonlinear systems, Cordero et al. ^[3] proposed a new method by replacing the Jacobi matrix with the division difference. There are many higher-order iterative methods that have been proposed; see ^{[4,5,6,7,8,9,10]}.

Nowadays, an effective way to analyze stability of iterative methods is to use real dynamics tools. Cordero et al. ^[11] studied the real dynamical analysis of fourth-order iterative methods on quadratic polynomials. Bakhtiari et al. ^[12] used a parameter for acceleration to obtain a new method with memory and analyzed its dynamic behavior. Amat et al. ^[13] performed a kinetic study of a classical third-order Newton-type iterative method when it is applied to polynomials of the second and third degree.

In this paper, a family of one-parameter Ostrowski-type iterative methods are proposed for solving a nonlinear system. In order to select the method with better stability, real dynamical tools are used to analyze the stability. Theoretical proof of the convergence order of the new iterative family is given in Section 2. In Section 3, the singular fixed points and critical points of the new method are calculated and its dynamical planes are plotted. We use real dynamics to analyze stability of the new iterative method, comparing the attractive basins of different parameter values so as to determine which parameter values of the iterative method is most stable. In Section 4, some nonlinear problems are subjected to numerical experiments. Numerical results are obtained to confirm the theoretical results. A short conclusion is derived in Section 5.

2.   A new family of iterative methods

Now, focus on the Ostrowski's fourth-order method ^[14], which can be written as

where $j = 0, 1, 2, \cdots.$ Now, we use Ostrowski's method to construct a new family of iterative methods. The iteration format is as follows:

where $j = 0, 1, 2, \cdots$ and $\lambda$ is any real number. When $\lambda = 0$, the iterative method (2.2) is the classical Ostrowski's method. The convergence order of family (2.2) is proved below.

Theorem 1. Let $F:E\subseteq {{\mathbb{R}}^{n}}\to {{\mathbb{R}}^{n}}$ be a sufficiently differentiable function in the open convex set $E$, and use ${{x}^{*}}\in E$ to represent the solution of $F\left(x \right) = 0$ such that ${F}'$ is continuous and non-singular in ${{x}^{*}}$. If the initial estimate ${{x}^{\left(0 \right)}}$ of ${{x}^{*}}$ is close enough, then for any value of $\lambda \in \mathbb{R}$, a member of family (2.2) converges to ${{x}^{*}}$ with fourth-order, and its error equation is as follows:

where ${{e}^{\left(j \right)}} = {{x}^{\left(j \right)}}-{{x}^{*}}$ is the error in $j$-th iteration and ${{C}_{i}} = \frac{1}{i!}{{\left[{F}'\left({{x}^{*}} \right) \right]}^{-1}}{{F}^{\left(i \right)}}\left({{x}^{*}} \right), i\ge 2.$

Proof. The Taylor's expansion of $F\left({{x}^{\left(j \right)}} \right)$ around ${{x}^{*}}$ can be written by

and its derivative is expressed as

From the equality ${{\left[{F}'\left({{x}^{\left(j \right)}} \right) \right]}^{-1}}{F}'\left({{x}^{\left(j \right)}} \right) = I$, the inverse of ${F}'\left({{x}^{\left(j \right)}} \right)$ can be calculated:

where,

${{H}_{1}} = -2{{C}_{2}}, \\ {{H}_{2}} = 4C_{2}^{2}-3{{C}_{3}}, \\ {{H}_{3}} = -8C_{2}^{3}+6{{C}_{2}}{{C}_{3}}+6{{C}_{3}}{{C}_{2}}-4{{C}_{4}}, \\ {{H}_{4}} = 16C_{2}^{4}+9C_{3}^{2}+8{{C}_{2}}{{C}_{4}}+8{{C}_{4}}{{C}_{2}}-12C_{2}^{2}{{C}_{3}}-12{{C}_{2}}{{C}_{3}}{{C}_{2}}-12{{C}_{3}}C_{2}^{2}-5{{C}_{5}}.$

Using (2.4) and (2.6), we can get

Similar to (2.5), we get

From (2.4) and (2.8), we can obtain

where,

${{P}_{1}} = {{C}_{2}}, \\ {{P}_{2}} = C_{2}^{2}+{{C}_{3}}, \\ {{P}_{3}} = -2C_{2}^{3}+3{{C}_{2}}{{C}_{3}}+{{C}_{4}}, \\ {{P}_{4}} = 4C_{2}^{4}-8C_{2}^{2}{{C}_{3}}+2C_{3}^{2}+4{{C}_{2}}{{C}_{4}}+{{C}_{5}}.$

From $\left[{{x}^{(j)}}, {{y}^{(j)}}; F \right]^{-1}\left[{{x}^{(j)}}, {{y}^{(j)}}; F \right] = I$, the inverse of $\left[{{x}^{(j)}}, {{y}^{(j)}}; F \right]$ can be calculated

where,

${{P}_{1}} = -{{C}_{2}}, \\ {{P}_{2}} = -{{C}_{3}}, \\ {{P}_{3}} = 3C_{2}^{3}-{{C}_{2}}{{C}_{3}}-{{C}_{4}}, \\ {{P}_{4}} = -9C_{2}^{4}+13C_{2}^{2}{{C}_{3}}-C_{3}^{2}-2{{C}_{2}}{{C}_{4}}-{{C}_{5}}.$

Finally, substituting (2.8), (2.9), and (2.10) into the family (2.2), the error equation can be obtained as

So, family (2.2) is fourth-order convergence for any $\lambda \in \mathbb{R}.$ □

3.   Stability analysis

This section analyzes the dynamical behavior of family (2.2) to investigate whether the stability of iterative families is related to the choice of parameter values $\lambda$. To this end, the stability of this family is analyzed by real multidimensional dynamical tools. Dynamics in the multidimensional case is a relatively new field of research ^[15,16,17]. More similar studies and information can be found in a number of studies; see ^[18,19,20]. In order to analyze the stability of the family (2.2), we use a real multidimensional dynamics analysis method, where the convergence is judged from initial estimates. Some basics are presented first to help understand the study.

3.1.   Fundamentals of multidimensional real dynamics

It is usually a rational function $X:{{\mathbb{R}}^{n}}\to {{\mathbb{R}}^{n}}$, obtained by applying the iterative method to the $n$-variable polynomial $q\left(x \right)$, $q:{{\mathbb{R}}^{n}}\to {{\mathbb{R}}^{n}}$, such that it represents the rational vector function related to the iterative method of the polynomial $q\left(x \right)$. The orbit of a point ${{x}^{\left(0 \right)}}\in {{\mathbb{R}}^{n}}$ is defined to be the set of consecutive applications of $\mathbb{R}:$

Obviously, there are different progressive behaviors when performing a dynamical analysis of initial points in ${{\mathbb{R}}^{n}}$, and they can be classified as such. When $X\left(\widetilde{x} \right) = \widetilde{x}$, $\widetilde{x}\in {{\mathbb{R}}^{n}}$ is the fixed point of $X$. A fixed point that is different from the root of the polynomial $q\left(x \right)$ is called a strange fixed point. Indeed, if a point $\widetilde{x}\in {{\mathbb{R}}^{n}}$ satisfies ${{X}^{k}}\left(\widetilde{x} \right) = \widetilde{x}$, then ${{X}^{p}}\left(\widetilde{x} \right)\ne \widetilde{x}$ ($p = 1, 2, \cdots, k-1$) is called a $k$-periodic point. It can be noted that if $k = 1$, then the $k$-periodic point $\widetilde{x}$ is a fixed point.

Theorem 2. Let $X:{{\mathbb{R}}^{n}}\to {{\mathbb{R}}^{n}}$ be ${{C}^{2}}$. Let us also assume that $\widetilde{x}$ is a $k$-periodic point, and let ${{\lambda }_{1}}, {{\lambda }_{2}}, \ldots, {{\lambda }_{m}}$ be the eigenvalues of ${X}'\left(\widetilde{x} \right)$, then

(1) $\widetilde{x}$ is attracting if $\left| {{\lambda }_{j}} \right| < 1$, for all $j = 1, 2, \cdots, m.$

(2) If $\exists {{j}_{0}}\in \left\{ 1, 2, \cdots, m \right\}$ such that $\left| {{\lambda }_{{{j}_{0}}}} \right| > 1$, then $\widetilde{x}$ is unstable (repelling or saddle).

(3) $\widetilde{x}$ is repelling if $\left| {{\lambda }_{j}} \right| > 1$, for all $j = 1, 2, \cdots, m.$

Furthermore, fixed points are defined as hyperbolic points when all the eigenvalues ${{\lambda }_{j}}$ satisfy $\left| {{\lambda }_{j}} \right|\ne 1$. The fixed points are called saddle points. In particular, this hyperbolic point is defined as a saddle point if there are eigenvalues ${{\lambda }_{j}}$ leading to $\left| {{\lambda }_{j}} \right| > 1$ and eigenvalues ${{\lambda }_{i}}$ leading to $\left| {{\lambda }_{i}} \right| < 1$. This fixed point is the superattractive point, when all eigenvalues are equal to zero.

Let the basin that attracts the fixed point ${{x}^{*}}$ be the set of points whose orbits converge to ${{x}^{*}}$, expressed by $A\left({{x}^{*}} \right)$,

Finally, if the point $x\in {{\mathbb{R}}^{n}}$ of $X$ satisfies $\frac{\partial {{r}_{i}}\left(x \right)}{\partial {{x}_{j}}} = 0$, then $x$ is a critical point, for all $i, j = 1, 2, \cdots, m$, being ${{r}_{i}}\left(x \right)$ the coordinate functions of the vectorial rational operator $X$. The superattracting fixed point will also be the critical point. A critical point that does not match the root of $q\left(x \right)$ is called a free critical point.

3.2.   Dynamical analysis of new family

Now, the actual dynamical analysis of the family (2.2) is applied to the polynomial system. It will be studied in two dimensions in order to graphically visualize the basin of attraction on the actual plane. In addition, the dynamical properties of the two-dimensional polynomial iterative method can be expanded to similar $n$-dimensional polynomials. The systems considered in the stability analysis are the following polynomials:

Applying the family (2.2) to a square polynomial $q\left(x \right)$ yields vector rational operators with the following expression:

where

We can note that for some specific values of the parameter $\lambda$, this rational function simplifies highly. For example, for $\lambda = 0$, the resulting component is:

A summary of the fixed point stability study of $X\left(x, \lambda \right)$ appears in Theorem 3 below.

Theorem 3. The theoretical function $X\left(x, \lambda \right)$ related to the iterative method family (3) has ${{2}^{n}}$ superattractive fixed points whose components are the roots of $q\left(x \right)$. It has also a number of real fixed points distinct from the roots, and the components of the roots are combined from the roots of the polynomials $l\left(t \right) = -\lambda +\left(2+4\lambda \right){{t}^{2}}+\left(18-6\lambda \right){{t}^{4}}+\left(38+4\lambda \right){{t}^{6}}+\left(6-\lambda \right){{t}^{8}}$, depending on $\lambda$, and the roots of $q\left(x \right)$:

(1) If $\lambda < -\frac{1}{2}$, there is no real root here.

(2) If $\lambda \geqslant\frac{1}{2}$, there exist two real roots of $l\left(t \right)$, denoted by ${{l}_{i}}\left(a \right)$, $i = 1, 2.$

Proof. The fixed point of $X\left(x, \lambda \right)$ can be obtained by solving for $X\left(x, \lambda \right) = x$, that is, solving

That is, when $t \ne 0,$ the component ${{x}_{j}} = \pm 1$ of the fixed point is also one root of the polynomial $l\left(t \right)$.

The roots of $l\left(t \right)$ are denoted as ${{l}_{i}}\left(\lambda \right)$, $i = 1, 2, \cdots, 8.$ Up to six roots of $l\left(t \right)$ are real when the parameter $\lambda$ takes different values. The stability of the fixed point $X\left(x, \lambda \right)$ is provided by the integral value of the relevant Jacobi matrix eigenvalues at the fixed point. Due to the polynomial system nature, the values of the eigenvalues coincide to the coordinate functions of rational operators:

being ${{\nu }_{j}}\left(\lambda \right) = -1+{{l}_{j}}{{\left(\lambda \right)}^{2}}$, ${{\alpha }_{1}}\left(\lambda \right) = 1+4{{l}_{j}}{{\left(\lambda \right)}^{2}}+3{{l}_{j}}{{\left(\lambda \right)}^{4}}$, ${{\beta }_{1}}\left(\lambda \right) = {{l}_{j}}\left(\lambda \right)+3{{l}_{j}}{{\left(\lambda \right)}^{3}}$ and ${{\beta }_{2}}\left(\lambda \right) = 3+25{{l}_{j}}{{\left(\lambda \right)}^{2}}+33{{l}_{j}}{{\left(\lambda \right)}^{4}}+3{{l}_{j}}{{\left(\lambda \right)}^{6}}$. Based on the consideration of the absolute values of the eigenvalues at different immovable points in the real interval, it can be shown that immovable points with component $\pm 1$ are superattractive. □

The dynamical analysis in this paper is to select the optimal value of parameter $\lambda$ from the stability aspect, and the corresponding iterative method has good stability. Now, the Jacobi matrix $X'\left(x, \lambda \right)$ of the rational function, and its critical points are analyzed. A critical point is known to be a solution of $det\left(X'\left(x, \lambda \right)\right) = 0$. A critical point is defined as a free critical point when the critical point is not a solution of $q\left(x \right) = 0$.

Theorem 4. The free critical points of operator $X\left(x, \lambda \right)$, denoted by $\left(c{{r}_{1}}, c{{r}_{2}}, \ldots, c{{r}_{n}} \right), n\le 12$, are the components that make the Jacobian matrix ($j = 1, 2, \cdots, n$), in which all components different from the roots of $q\left(x \right)$ are zero, that is:

(1) If $-\infty < \lambda < 0$, there are no free critical points.

(2) If $0 \le \lambda < 6$, then $c{{r}_{1}}\left(\lambda \right)$ and $c{{r}_{2}}\left(\lambda \right)$ are the different components of the free critical points, being ${z}^{*}$ is the positive root of polynomial.

(3) If $6\le\lambda$, there are no free critical points.

Proof. Different from zero components of the Jacobian matrix $X'\left(x, \lambda \right)$ are:

Therefore, it is obvious that the roots of $2{{d}^{2}}+12{{d}^{4}}+18{{d}^{6}}-\lambda \left(3+25{{d}^{2}}+33{{d}^{4}}+3{{d}^{6}} \right)$, when they are real, are the components of the critical points. Solving for $k\left(t \right)$ by making $t = {{d}^{2}}$ defines the free critical point. □

To visualize the attractive basin of the family (2.2), the parameter planes are plotted with different parameter $\lambda$ values. In Figure 1, the set of initial iterations in family (2.2) is represented by the ${{X}_{1}}$-axis and ${{X}_{2}}$-axis. The initial estimate corresponds to each point $\left({{x}_{1}}, {{x}_{2}} \right)$ on the plane. In the dynamical plane, it can be observed that the set of normal mappings consists of all points whose orbits converge to an attractive fixed point. Theorems 3 and 4 state that the singular fixed point asymptotic behavior and the amount of free critical points depend on the parameter $\lambda$. For the purpose of comparing the stability of different family (2.2), different parameter values are chosen. Figure 1 shows the dynamical planes of a set of initial estimates $\left({{x}_{1}}, {{x}_{2}} \right)\in \left[-1.5, 1.5 \right]\times \left[-1.5, 1.5 \right].$ To realize the dynamical planes, we define convergence as when the orbit of $\left\| \left({{x}_{1}}, {{x}_{2}} \right)-\left(\pm1, \pm1 \right) \right\| < {{10}^{-3}}$ or a point reaches 50 iterations without converging to the roots of the polynomial, and denote this by a white star. When a point's orbit converges to different roots of $q\left(x \right)$, the points in the plane are represented in blue, green, orange, and red, and otherwise in black.

First, the basin of attraction displayed in the dynamical planes are shown in Figure 1. In this case, Figure 1(f) represents the dynamical plane associated with $\lambda = 0$. Moreover, Figure 1(a) and 1(e) represent the dynamical plane associated with $\lambda = -1$ and $\lambda = 10$, when the black region of non-convergence is very large. Figure 1(b) and 1(d) represent the dynamical planes associated with $\lambda = \frac{1}{2}$ and $\lambda = 8$, respectively, where there are fewer black regions that do not converge. Figure 1(c) represents the dynamical plane associated with $\lambda = \frac{5}{7}$, where there are no non-convergence black regions.

As can be seen from Figure 1, considering the dynamical plane related to the parameters, when the parameter is $\lambda < 0$ and $\lambda > 6$, the non-convergence black region is large and the stability is poor. In the case of parameter $0\le\lambda \le 6$, the non-convergence black region is small and the stability is better; where the parameter is $\lambda = \frac{5}{7}$, the stability is the best.

4.   Numerical experiments

In this section, a number of numerical problems are used to demonstrate the computational efficiency and convergence behavior of the method proposed in this paper. In order to test its performance, several of the iterative family (2.2) are chosen to solve the nonlinear problems. It can be noticed that the family (2.2) for $\lambda = 0$ is Ostrowski's method. In Tables 1–3, iter represents the number of iterations, the errors values are represented by $\left\| {{x}^{\left(k+1 \right)}}-{{x}^{\left(k \right)}} \right\|$, and $\left\| F\left({{x}^{\left(k+1 \right)}} \right) \right\|$ represents the errors of function values at the last step. The approximate computational order of convergence (ACOC) presented in ^[21] is used to numerically estimate the orders of the iterative method, defining it as:

All numerical calculations use a stopping criterion of $\left\| {{x}^{\left(k+1 \right)}}-{{x}^{\left(k \right)}} \right\| < {{10}^{-100}}$ or $\left\| F\left({{x}^{\left(k+1 \right)}} \right) \right\| < {{10}^{-100}}$ to minimize rounding errors. Three nonlinear problems are selected here for numerical comparison.

Example 4.1. Hammerstein equations ^[22] are a class of nonlinear equations that are very important in practical applications and are formulated as follows:

where $x\in \mathbb{C}\left[0, 1 \right], r, s\in \left[0, 1 \right],$ with the kernel $N$ as

Equation (4.2) is transformed into a nonlinear system by a discretization process. The Gauss-Legende quadratic is used to approximate the integrals appearing in Eq (4.2),

where ${{s}_{j}}$ and ${{w}_{j}}$ are the nodes and weights of the Gauss-Legan polynomial, respectively. ${{x}_{i}}$ ($i = 1, \ldots, 7$) is an approximate representation of $x\left({{s}_{j}} \right)$. We then use a system of nonlinear equations to estimate (4.2)

where

So, the system can be written as

In the above equation, $F$ is the nonlinear operator, which is in the Banach space ${{\mathbb{R}}^{L}}$, while ${F}'$ is its Fréchet derivative in $L = \left({{\mathbb{R}}^{L}}, {{\mathbb{R}}^{L}} \right)$. By method (2.2), the solution of the nonlinear system will be found. By taking ${{x}^{\left(0 \right)}} = {{\left(1.8, 1.8, \cdots, 1.8 \right)}^{T}}$, Table 1 can be obtained.

Example 4.2. Boundary-value problem:

We discretize the first and second derivatives in this problem by using the difference method

and

Split the interval $\left[0, 1 \right]$ into $n$ smaller intervals with endpoints $0 = {{x}_{0}} < {{x}_{1}} < \ldots < {{x}_{n-1}} < {{x}_{n}} = 1$. The splitting is regular, that is, $\vartriangle {{x}_{t}} = {1}/{n}\;$ for all $t$. The following nonlinear system can be found:

The solution ${{\left(0.01073, 0.02147, 0.03220, \ldots, 0.9928 \right)}^{T}}$ is built from the initial value ${{x}^{\left(0 \right)}} = {{\left(0.5, 0.5, \ldots, 0.5 \right)}^{T}}$ for $n = 101$. Table 2 displays the numerical results.

Example 4.3. Nonlinear system:

The solution $\varsigma = {{\left\{ 1, 1, \ldots, 1 \right\}}^{t}}$ is sought by choosing the initial guess value ${{x}^{\left(0 \right)}} = {{\left\{ 5, 5, \cdots, 5 \right\}}^{t}}$. Parameter $n = 99$ is chosen. The numerical results are shown in Table 3 below.

5.   Conclusions

In this paper, a new iterative family (2.2) is proposed. When $\lambda = 0$, the new iterative method is the Ostrowski's method (2.1). It can be calculated theoretically that every member of this family is fourth-order convergence. Since different parameter values lead to different situations, their stability is analyzed using real multidimensional dynamical tools. The number of strange fixed points and critical points for different parameters $\lambda$ are analyzed and the dynamical planes are plotted. This is used to select the parameter values with better performance. In numerical experiments, some specific members are chosen to be applied to the nonlinear problems. It can be concluded from the dynamical plane and numerical experiments that the stability is better when the parameter is $\lambda = \frac{5}{7}$. In this paper, the stability of the new iterative method is only analyzed by real dynamical tools. We will analyze stability with complex dynamics tools in the future.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

This research was supported by the National Natural Science Foundation of China (No. 61976027), the Natural Science Foundation of Liaoning Province (Nos. 2022-MS-371, 2019-ZD-0502), Educational Commission Foundation of Liaoning Province of China (Nos. LJKMZ20221492, LJKMZ20221498), and the Key Project of Bohai University (No. 0522xn078).

Conflict of interest

The authors declare no conflict of interest.