A single parameter fourth-order Jarrat-type iterative method for solving nonlinear systems

1.   Introduction

In the field of engineering and scientific computing, solving systems has become a core task. Iterative algorithms have become the preferred method for solving such problems due to their efficiency and flexibility in estimating solutions to nonlinear systems. The iterative method solves the nonlinear problem by gradually approaching the estimated value of the iterative adjustment solution until it converges to the exact solution or meets the predetermined error tolerance.

The Jarrat-type iterative method has been optimized and improved based on the Newton iteration method, and is an effective method for solving nonlinear systems ^[1]. The classic Jarrat-type iterative method takes the form of (1.1):

where $\Gamma$ is a function defined on $R^{n}$, $\nu \in R^{n}$ and $S = \left [6\Gamma '(\hbar^{(k)})-2\Gamma '(\nu^{(k)}) \right] ^{-1}\left [3\Gamma '(\hbar^{(k)})+\Gamma' (\nu^{(k)}) \right]$, comparing the convergence order of Newton's iteration method and the classical Jarrat-type iteration method, it is evident that the Jarrat-type iteration method has better convergence.

In order to further improve the convergence speed and order of the Jarrat-type iteration method, and expand their scope of use, Cordero et al. ^[2], denoted by S1, proposed a two-step fourth-order Jarrat-type iterative method:

Francisco et al. ^[3] combined the Newton method with the Traub method and proposed two new methods for solving nonlinear systems, which were named S2 and S3, respectively:

and

where $\sigma _{k} = \left [\Gamma'\left (\nu^{(k)} \right) \right]^{-1} \left (\Gamma'\left (\nu^{(k)} \right)-\left [\hbar ^{\left (k \right) }, \nu ^{\left (k \right) }; F \right] \right).$ Inspired by the above methods, this article proposes a fourth-order iterative algorithm based on weight functions. Compared to traditional Jarrat-type iterative methods, this new algorithm exhibits superior convergence performance in solving nonlinear system problems.

This article proposes a single-parameter Jarrat-type iterative method for solving nonlinear systems. In Section 2, the convergence of iterative methods has been rigorously proven within the specific mathematical framework of Banach spaces. Based on the proven convergence conclusions of the iterative method, the stability of the iterative method was analyzed in depth in Section 3 of this article. At the same time, the number of fixed points and critical points was accurately calculated through corresponding mathematical operations. In Section 4, the dynamic plane under different parameters is drawn. The stability of the iterative method under different parameters is analyzed by using the real dynamics method. Determine the parameter values with better stability by comparing different attractive basins. In Section 5, the efficiency of the iterative method designed in this paper is compared with two iterative methods of the same order. In Section 6, four numerical experiments were conducted on the iterative method, and the obtained numerical results further proved that the parameter values of the iterative method have better stability under $0\le\beta < \frac{423}{8}$. Finally, in Section 7 of the article, a summary was made.

2.   Convergence analysis of Jarrat-type iterative method

According to ^[4], a new method for solving nonlinear systems has been proposed:

where $\lambda = \left [\Gamma'\left ({{\hbar}^{(k)}} \right) \right] ^{-1}\Gamma'\left({{\nu}^{(k)}} \right).$ And

where $k = 0, 1, 2, \cdots$, and $G(\lambda)$ is a weight function. In particular, when parameter $\beta = \frac{9}{8}$, the iterative method will be as shown in (1.2). When parameter $\beta = 0$, the iterative method is in the form of S4:

and iterative method (2.1) has a fourth order of convergence, which is proved by the following theorem.

Theorem 2.1. Let $\Gamma:I\subseteq {{\mathbb{R}}^{n}}\to {{\mathbb{R}}^{n}}$ be a sufficiently differentiable function in the open convex set $I$. Let us consider that ${{\nu}^{*}}\in I$ is a root of $\Gamma\left(\nu \right) = 0$ and the initial estimate ${{\nu}^{\left(0 \right)}}$ is sufficiently close to ${{\nu}^{*}}$, then the family (2.1) is fourth-order convergent, and the error equation is as follows:

where ${{C}_{w}} = \frac{1}{w!}{{[\Gamma'({{\nu}^{*}})]}^{-1}}{{\Gamma}^{(w)}}({{\nu}^{*}})$, $w\ge 2$, and ${{\kappa }^{\left(k \right)}} = {{\nu}^{\left(k \right)}}-{{\nu}^{*}}$ is the error in $k$-th iteration.

Proof. The Taylor expansion of $\Gamma\left({{\nu}^{\left(k \right)}} \right)$ and $\Gamma'\left({{\nu}^{\left(k \right)}} \right)$ at point ${{\nu}^{*}}$ is as follows:

and its derivative is expressed as

From the equality ${{\left[{\Gamma}'\left({{\nu}^{\left(k \right)}} \right) \right]}^{-1}}{\Gamma}'\left({{\nu}^{\left(k \right)}} \right) = I$, ${{\left[{\Gamma}'\left({{\nu}^{\left(k \right)}} \right) \right]}^{-1}}$ can be calculated:

where

Using (2.4) and (2.6), we obtain

Thus,

The first derivative of $\Gamma\left({{\hbar}^{\left(k \right)}} \right)$ can be found as

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

where

According to (2.8) and (2.9), we obtain

Thus,

where $e\left(\beta\right) = 2\left(\left(-\frac{10}{27}\beta-\frac{1}{6} \right)C_{2}^{3}-\frac{3}{4}{{C}_{2}}{{C}_{3}}-\frac{3}{12}{{C}_{4}} \right).$ Substituting the above equations into the family (2.1), the error equation can be obtained:

3.   Stability analysis

In this section of the research, the focus is on exploring the iterative method family (2.1), with the core focus on whether its stability is related to the value of parameter $\beta$. To deeply analyze this key issue, a targeted analysis method was specifically selected, which is to systematically and meticulously analyze the stability of high-order iterative methods in real dynamic planes using dimensional quadratic polynomials. It is worth emphasizing that when conducting stability analysis for the iterative method family (2.1), we strictly follow the rigorous multidimensional dynamic analysis method, which is a scientific and standardized path. Under the analytical framework constructed by this method, the determination of convergence is carried out by referring to the preliminary estimation process detailed in academic works ^[5,6,7], in order to ensure the logical rigor and reliability of the entire analysis process, and provide solid theoretical support for stability-related research in the iterative method family (2.1).

3.1.   The basic definition of real dynamics analysis

Applying the n-dimensional rational vector polynomial $p\left (\nu\right)$, $p:\mathbb{R}^{n}\to \mathbb{R}^{n}$, and the iterative method (2.1), the rational function $E$ can be obtained: $\mathbb{R}^{n}\to \mathbb{R}^{n}$, a orbit point is defined as ${{\nu}^{\left(0 \right)}}\in {{\mathbb{R}}^{n}}$, defining a set of $\mathbb{R}:$$\left\{ {{\nu}^{\left(0 \right)}}, E\left({{\nu}^{\left(0 \right)}} \right), {{E}^{2}}\left({{\nu}^{\left(0 \right)}} \right), \cdots, {{E}^{n}}\left({{\nu}^{\left(0 \right)}} \right) \right\}$.

It is evident that distinct progressive dynamics are exhibited when conducting a dynamical analysis of initial points within the n-dimensional Banach space. These behaviors can be systematically categorized as follows: A point $\bar{\nu}$ is a fixed point of $E$ if $E\left({\bar{\nu}} \right) = \bar{\nu}$. Particularly, if $\bar{\nu}$ is the solution of the equation $E\left({\bar{\nu}} \right)\ne \bar{\nu}$, $\bar{\nu}$ is a strange fixed point. Furthermore, a point $\bar{\nu}\in {{R}^{n}}$ satisfies ${{E}^{k}}\left({\bar{\nu}} \right) = \bar{\nu}$, then ${{E}^{i}}\left({\bar{\nu}} \right)\ne \bar{\nu}$ ($i = 1, 2, \cdots, k-1$) is called an $i$-periodic point; see ^[8,9,10].

Theorem 3.1. Consider the function $E:{{\mathbb{R}}^{n}}\to {{\mathbb{R}}^{n}}$. Let us presuppose that $\bar{\nu}$ is an $i$-periodic point that considers the function $E:{{\mathbb{R}}^{n}}\to {{\mathbb{R}}^{n}}$. Denote the eigenvalues of the Jacobian matrix $E'\left({\bar{\nu}} \right)$ by ${{\xi }_{1}}, {{\xi }_{2}}, \ldots{{\xi }_{n}}$, then:

(1) $\bar{\nu}$ is superattracting if $\left| {{\xi }_{i}} \right| = 0$, for all $i = 1, 2, \cdots, n.$

(2) $\bar{\nu}$ is an attracting if $\left| {{\xi }_{i}} \right| < 1$, for all $i = 1, 2, \cdots, n.$

(3) $\bar{\nu}$ is neutral if $\left| {{\xi }_{i}} \right| = 1$, for all $i = 1, 2, \cdots, n.$

(4) $\bar{\nu}$ is repelling if $\left| {{\xi }_{i}} \right| > 1$, for all $i = 1, 2, \cdots, n.$

(5) $\bar{\nu}$ is a hyperbolic point or saddle point if $\left| {{\xi }_{i}} \right|\ne 1$, for all $i = 1, 2, \cdots, n.$

(6) $\bar{\nu}$ is a superattractive point if $\left| {{\xi }_{i}} \right| = 0$, for all $i = 1, 2, \cdots, n.$

Allow the basin that draws the fixed point $\tilde{\nu}$ to be the collection of points whose track tends towards $\tilde{\nu}$, denoted by $A\left({\tilde{\nu}} \right)$,

Finally, critical points satisfy $\frac{\partial {{r}_{i}}\left(\nu \right)}{\partial {{\nu}_{i}}} = 0$, $\nu\in {{\mathbb{R}}^{n}}$ of $V$, for all $i = 1, 2, \cdots, n$. According to (6) in Theorem 3.1: It can be concluded that the superattractive fixed points also satisfy the definition of critical points.

3.2.   Stability analysis of fixed points

This part mainly relates the family of iterative methods (2.1) to polynomial systems. From the perspective of two-dimensional space, the graphical visualization of the basins of attraction on the real dynamic plane is analyzed. And generalize the dynamical properties of two-dimensional polynomials to $n$-dimensional polynomials. The system under scrutiny within the range of stability analysis is delineated by the ensuing polynomial expression:

By employing the family (2.1) to a quadratic polynomial $p\left(\nu \right)$, one derives vector rational operators with the ensuing formulation:

where

The stability of the fixed point $X\left(\nu, \beta \right)$ will be summarized in Theorem 3.2 next.

Theorem 3.2. There will be a combination of $V\left(\nu, \beta \right)$ and the representation method (2.1), with ${{2}^{n} }$ super-attractive fixed points, whose component is the root of $p\left(\nu \right)$. Let the function $V\left(\nu, \beta\right) = \nu$, and the real roots obtained are called fixed points. At this point, the fixed points are also the roots of $p\left(\nu \right)$ and the polynomials $l\left(t \right) = 8\beta {{\left(-1+{{t}^{2}} \right)}^{3}}-9{{t}^{2}}\left(5+20{{t}^{2}}+47{{t}^{4}} \right)$ that depend on $\beta$:

(1) If $\beta < 0$, there are two real roots, ${{l}_{1}}\left(\beta \right)$ and ${{l}_{2}}\left(\beta \right)$.

(2) If $0\le \beta < \frac{423}{8}$, there is no real root here.

(3) If $\beta \ge \frac{423}{8}$, there are two real roots, ${{l}_{3}}\left(\beta \right)$ and ${{l}_{4}}\left(\beta \right)$.

Proof. The solution of the rational function $V\left(\nu, \beta \right) = \nu$ is the value of a fixed point:

It is actually solving for $\left(-1+{{\nu}_{j}}^{2} \right)\left(8\beta {{\left(-1+{{\nu}_{j}}^{2} \right)}^{3}}-9\nu_{j}^{2}\left(5+20\nu_{j}^{2}+47{{\nu}_{j}}^{4} \right) \right)$. Apparently, ${{\nu}_{j}} = \pm 1$ is a fixed-point polynomial of ${{p}_{j}}\left(\nu \right) = \nu_{j}^{2}-1$. Meanwhile, when $t\ne 0$, there are six roots on the $l\left(t \right)$ component, but only ${{l}_{1}}\left(\beta \right)$ and ${{l}_{2}}\left(\beta \right)$ are real roots.

Since the highest degree term of polynomial $l\left(t \right)$ is sixth, let ${{l}_{1}}\left(\beta \right), {{l}_{2}}\left(\beta \right), \cdots, {{l}_{6}}\left(\beta \right)$ be the six real roots of $l\left(t \right)$ as the root of polynomial $l\left(t \right) = 0$, where ${{l}_{i}}\left(\beta \right)$, $i = 1, 2, \cdots, 6$ as the six real roots of $l\left(t \right)$ is related to parameter $\beta$. The equilibrium of the fixed point \(V\left(\nu, \beta\right) \) is ascertained by the sum of the integral eigenvalues of the corresponding Jacobian matrix at the fixed point. Owing to the intrinsic properties of the polynomial system, the eigenvalues align with the coordinate functions of the rational operators. The equilibrium state of a fixed point $V\left(\nu, \beta \right)$ is determined by the sum of the integral eigenvalues of the corresponding Jacobian matrix at that fixed point. Given the inherent properties of polynomial systems, eigenvalues are obtained when the coordinate functions of rational operators exhibit relative homogeneity.

where ${{u}_{j}}\left(\beta \right) = -1+{{l}_{j}}{{\left(\beta \right)}^{2}}$, ${{\alpha }_{j}}\left(\lambda \right) = 5+20{{\left({{l}_{j}}\left(\beta \right) \right)}^{2}}+47{{\left({{l}_{j}}\left(\beta \right) \right)}^{4}}$ and ${{h}_{j}}\left(\beta \right) = 1+2{{l}_{j}}{{\left(\beta \right)}^{2}}$. In real space, consider the absolute values of fixed point eigenvalues. Therefore, fixed points with $\pm 1$ components exhibit superattractive fixed points. □

3.3.   Stability analysis of critical point

In this section, if we want to achieve optimal stability in the iterative method of dynamic analysis, we must determine the most favorable value for the parameter $\beta$. Therefore, the first-order derivative of the rational function $V\left(\nu, \beta \right)$ can be obtained to analyze the sum critical point of the Jacobian matrix $V'\left(\nu, \beta \right)$. The free critical point is the root that satisfies both $V'\left(\nu, \beta \right) = 0$ and $p\left(\nu\right)\ne 0$ simultaneously.

Theorem 3.3. The free critical point $\left(c{{r}_{1}}, c{{r}_{2}}, c{{r}_{3}} \right)$ constitutes a solution to the equation $V'\left(\nu, \beta \right) = 0$, but it does not satisfy the polynomial equation $p\left(\nu\right) = 0$, that is:

(1) If $-\infty < \beta < 0$, the expressions $cr_1 = -\sqrt{r^{\ast}}$ and $cr_2 = \sqrt{r^{\ast}}$ delineate the distinct components of the free critical points. Here, $r^{\ast}$ signifies the positive root of the polynomial

(2) If $\beta = 0$, there is a free critical point with $cr_{3} = 0$.

(3) If $\beta > 0$, there dose not exist free critical points.

Proof. If $V'\left(\nu, \beta \right) = 0$, two components can be multiplied, and the free critical point is the component other than $p\left(\nu\right)$, as follows:

Therefore, it is obvious that the free critical point is the root of $24\beta +\left (45+152\beta \right) \sigma +\left (306+16\beta \right)\sigma ^{2}$, when the values of the roots are real numbers, they are the components of the critical point. Let $t = {{\sigma}^{2}}$, and the free critical point can be defined by solving $k\left(t \right)$. □

4.   Stable and unstable dynamic plane under parameter $\beta$

In order to conduct in-depth visualization analysis of the dynamic behavior of iterative algorithms, this study adopted the method of drawing the attractive basin map of the method family under different parameter values $\beta$. By careful observation on the dynamic plane, we can evaluate the stability characteristics of specific members in the iterative algorithm family. This method not only reveals the sensitivity of the algorithm to parameter changes, but also provides an intuitive basis for identifying parameter configurations with high stability. In Figure 1, we show the initial iteration set of family (2.1), which is defined by the $X_1$-axis and the $X_2$-axis in a two-dimensional coordinate system. On this coordinate plane, each specific point $(x_1, x_2)$ represents a guess value during the initial iteration process. These points constitute the geometric representation of the initial iteration set, and their distribution on the plane provides us with an intuitive understanding of the starting point of the iteration process. Each point corresponds to a potential starting position in the iterative algorithm, from which the iterative process will evolve according to the defined mapping rules. By systematically analyzing these initial guess points, we can further explore the dynamic behavior and convergence characteristics of each member in the family (2.1). Within the framework of dynamic plane analysis, it is worth noting that the set of standard mappings encompasses all points whose trajectories tend towards attractive fixed points. Specifically, the dynamic evolution of these points will ultimately stabilize at a fixed point in the system, which is a key indicator for system stability analysis. Determine the existence of singular fixed points and the correlation between the free critical point and the value of parameter $\beta$ through Theorems 3.2 and 3.3. To assess the stability across various members of the family (2.1), distinct parameter values are selected. We use graphical visualization tools to construct a $400\times 400$ grid to display the dynamic plane, where each point on the grid is an initial estimated set $\left({{x}_{1}}, {{x}_{2}} \right)\in \left[-2, 2 \right]\times \left[-2, 2 \right].$

In order to achieve intuitive visualization in the dynamic plane, we define specific convergence criteria: When the orbit of a point satisfies the condition $\left\| \left({{x}_{1}}, {{x}_{2}} \right)-\left(\pm1, \pm1 \right) \right\| < {{10}^{-3}}$, or when the point has not converged to the root of the polynomial $p\left(\nu \right)$ after 50 iterations, we mark the point with a white asterisk. This standard ensures that we can identify points that fail to achieve the expected convergence behavior. Furthermore, when the trajectory of a point converges to different roots of a polynomial $p\left(\nu \right)$, these points will be represented in blue, green, orange, and red on the plane to distinguish their corresponding roots. Points that fail to converge to any specific root are represented in black, which helps us to conduct detailed analysis and interpretation of dynamic behavior, as described in references ^{[11,12,13,14]}. Through this color coding, we can clearly identify points with different convergence properties, thereby gaining a deeper understanding of the structure and behavior of dynamic systems.

Next, the attractive basin in the dynamic plane is shown in Figure 1. Figures 1(a) and 1(b) respectively represent the dynamic planes of parameters $\beta = -5$ and $\beta = -1$, indicating that the convergence domains in these two figures are irregular. The convergence effect is not as good as shown in Figures 1(c) and 1(d). Figure 1(c) corresponds to the dynamic plane with parameter $\beta = 1$, and Figure 1(d) corresponds to the dynamic plane with parameter $\beta = 3$, both of which do not have black convergence regions. The last two figures, Figures 1(e) and 1(f), correspond to the dynamic planes of parameter $\beta = 55$ and parameter $\beta = 60$, respectively. There are obvious black areas of non-convergence in these two planes, especially in the $\left({{x}_{1}}, {{x}_{2}} \right)\in \left[-0.5, 0.5 \right]\times \left[-0.5, 0.5 \right]$ area in Figure 1(f). Therefore, according to Figure 1, when the parameters $\beta < 0$ and $\beta\ge \frac{423}{8}$, there are other points of attraction and non-convergence regions, showing poor stability. Finally, when the parameter $0\le\beta < \frac{423}{8}$ is within the range, there is no non-convergence region and no other attraction point in the dynamic plane, and the iterative method has the best stability.

5.   Computational efficiency index

According to the calculation cost and function evaluation times, in order to select a more effective method, we compare the efficiency of several fourth-order iterative methods proposed in this paper. We can use the efficiency index $I = R^{\frac{1}{k} }$ defined by Cordero ^[15], where $p$ is the order of the method and $k$ is the number of new function evaluations required for each iteration of the method. The computational cost of solving linear equations depends on its size. Since the proposed methods can be used to solve large-scale equations, we compare the performance of these methods with the computational efficiency index introduced in ^[16], that is, $CT = R^{\frac{1}{k+op} }$, where OP is the number of operations per iteration (product sum quotient) and the term $CL$ denotes the computational cost. Next, in Table 1, the calculation cost of each operation of n-dimensional linear equations is explained.

Figure 2 shows the efficiency index of iterative methods S1–S4 in Table 1 when n takes different values. It is obvious from Figure 2 that the efficiency of iterative methods S1 and S4 is significantly higher than that of other iterative methods of the same order. Although the efficiency of these four iterative methods decreases with the increase of the n value, the efficiency of iterative methods S1 and S4 is still higher than that of the other two iterative methods.

6.   Numerical experiments

In this section of the study, we adopted the iterative method (2.1) to solve three nonlinear systems. Through this method, we conducted a detailed analysis of the computational accuracy and convergence of the iterative process under different parameter values $\beta$. This not only verifies the optimal range of parameters for the stability of the iterative method in the fourth part of the paper, but also better evaluates the performance of the iterative method (2.1). The variables $iter$, $\left \| \nu^{\left (k+1 \right) } -\nu^{\left (k \right) } \right \|$, $time$, and $\left \| \Gamma \left (\nu^{\left (k+1 \right) } \right) \right \|$ in Tables 1–3 represent the number of iterations, the error value, computing time, and the residual error related to the final function value, respectively. Using the approximate computational order of convergence (ACOC) described in ^[17], the convergence rate of the iterative method is numerically determined to be

Every numerical calculation follows the termination criteria $\left \| \nu^{\left (k+1 \right) } -\nu^{\left (k \right) } \right \|$ or $\left \| \Gamma \left (\nu^{\left (k+1 \right) } \right) \right \| < 10^{-100}$.

Example 6.1. Hammerstein-type integral equation and demonstrate the application of theoretical results. The Hammerstein-type nonlinear integral equation is a specific version of the Urysohn-type integral equation. The following is the form of the Hammerstein equation:

where $\nu\in\mathbb{C}\left(0, 1\right), \iota, \omega\in\left[0, 1\right],$ with the kernel $\Theta$ as

By using discretization techniques, Eq (5.1) is transformed into a nonlinear equation form for easy solution. Subsequently, this study used the Gauss-Legendre numerical integration method to approximate the integration in Eq (5.1).

In this study, variables $\omega_{i}$ and $\Psi_{i}$ were designated as nodes and weights of Gaussian Legendre polynomials, respectively. By constructing a nonlinear system, we represent the approximate value of $\nu_{i}$ and estimate Eq (5.2) accordingly. Specifically, for $i = 1, 2, \ldots, 7$, we use $\nu(\omega_{i})$ to approximate $\nu_{i}$, where $\nu_{i}$ is approximated by

where

Rewriting the system in one method is

According to (2.1), we will use $\Phi'$ to solve the corresponding nonlinear system. This process involves calculating the derivatives of nonlinear operators and applying them to the solution process of nonlinear systems. By taking $\nu_{0} = (1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5)^{T}$, we obtain numerical results in Table 2.

Example 6.2. Boundary-value problem:

This study uses the finite difference method to numerically discretize the first and second derivatives in the discussed problem. Through this method, we convert continuous derivative operations into discrete algebraic expressions, making numerical calculations and solutions easier.

and

In order to split the interval $\left[0, 1 \right]$ into $n$ smaller intervals, we first define the endpoints as $0 = {{\vartheta}_{0}} < {{\vartheta}_{1}} < \ldots < {{\vartheta}_{n-1}} < {{\vartheta}_{n}} = 1$. This splitting is regular, meaning that the length of each subinterval is equal, which is $\bigtriangleup_{{\vartheta}_{k}} = \frac{1}{n}$, for all $t$. On the basis of this rule splitting, we can discover the following nonlinear systems:

In practical applications, we may need to choose appropriate numerical methods and computational strategies based on the nature and requirements of specific problems to ensure that the solving process of nonlinear systems is both efficient and accurate. The solution ${{\left(0.01073, 0.02147, 0.03220, \ldots, 0.9928 \right)}^{T}}$ is built from the initial value ${{ \vartheta}^{\left(0 \right)}} = {{\left(1.5, 1.5, \ldots, 1.5 \right)}^{T}}$ for $n = 90$. Table 3 displays the numerical results.

Example 6.3. Let us examine the ensuing system of nonlinear equations (excerpted from Grau-Sánchez), which is articulated by the function:

we choose $n = 3$ and use $\nu^{\left (0 \right) }$ to represent the initial guess. The solution to this problem is $\omega ^{\ast } \approx \left (0.5671433, 0.5671433, 0.5671433 \right) ^{T}$, where we select the initial estimates $\left (0.5, 0.5, 0.5 \right) ^{T}$. The numerical performance of the iterative method under different parameter values is shown in Table 4.

Example 6.4. Nonlinear PDE problem ^[18]:

and the initial condition is

In numerical analysis, this problem is discretized into a nonlinear system via the finite difference. The partial differential equation (PDE) problem under consideration is a heat conduction problem. Specifically, the spatial domain $\left [0, 1 \right]$ and the temporal domain $\left [0, T \right]$ are partitioned into N uniform subintervals, yielding step sizes of $\left(\varepsilon = \frac{1}{N} \right)$ in the spatial direction and $\left(\theta = \frac{T}{N} \right)$ in the temporal direction. Let $\eta = \eta \left (x, \varsigma \right)$ denote the exact solution of the problem, where $\eta \left (x, \varsigma \right)$ is a function of the spatial variable $\left(x \right)$ and the temporal variable $\left(\varsigma \right)$. Utilizing the finite difference method, we obtain approximations $\eta_{i, j} \approx \eta\left (x_{i}, \varsigma_{i} \right)$, $\eta_{x} \left (x, \varsigma \right)\approx \frac{\eta \left (x+\varepsilon, \varsigma \right)-\eta \left (x-\varepsilon, \varsigma \right) }{2\varepsilon}$, $\eta_{\varsigma} \left (x, \varsigma \right)\approx \frac{\eta \left (x, \varsigma \right)+\eta \left (x, \varsigma-\theta \right) }{\theta }$, and $\eta_{xx} \left (x, \varsigma \right)\approx \frac{\eta \left (x+\varepsilon, \varsigma \right)-2\eta \left (x, \varsigma \right)+\eta \left (x-\varepsilon, \varsigma \right) }{\varepsilon^{2} }$ at discrete points within the domain. This process results in the formation of a nonlinear system, which is presented as follows:

for $i = 1, 2, \ldots, N-1$ and $j = 1, 2, \ldots, N$, a series of nonlinear systems of size $\left(N-1\right)$$\times$$\left(N-1\right)$ are generated for each fixed $j$. By varying the values of $N$ and $T$, the problem is solved using distinct numerical methods. The results are presented in Table 5, where $iter$ denotes the number of iterations required for convergence, and $Time$ represents the computational CPU time in seconds. The exact solution and the approximate solution of the problem are depicted in Figures 3 and 4, respectively, for the case where $T = 0.01$ and $N = 10$. Additionally, Figure 5 illustrates the absolute error for the problem.

According to the observations in Tables 1–5, we can conclude that the convergence order of the iterative method (2.1) can reach or approach the fourth-order convergence under different parameter values $\beta$. Especially in the range of $0 < \beta < \frac{423}{8}$, the iterative method shows high accuracy. In order to further confirm the optimal performance parameters of the iterative method (2.1), we compared the calculation time of this method and the other two fourth-order iterative methods. The experimental results show that when $\beta = \frac{9}{8}$ is the iterative method S1, the calculation time of the iterative method is the shortest, the error is the smallest, and the calculation efficiency is the best.

7.   Conclusions

In this paper, a fourth-order Jarrat-type iterative method for solving nonlinear systems is proposed. The algorithm optimizes the iterative formula (2.1) by introducing a parametric weight function, ensuring that the iterative process can achieve fourth-order convergence regardless of the parameter value. In the stability analysis part, we calculate the fixed point and critical point, and use the critical point to analyze the dynamic plane. In numerical experiments, we solve Hammerstein-type integral equations, boundary value problems, and partial differential equations. The experimental results show that when the parameters meet $0\le\beta < \frac{423}{8}$, the iterative method shows excellent stability. In addition, through the comparative experiment with two kinds of fourth-order iterative methods, we found that the iterative method showed obvious advantages in calculation time, cost, and efficiency when the parameter $\beta = \frac{9}{8}$.

Author contributions

Jia Yu: Writing–original draft, Writing–review and editing; Xiaofeng Wang: Writing-review and editing. All authors have read and agreed to the published version of the manuscript.

Use of Generative-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 Educational Commission Foundation of Liaoning Province of China (No. LJ212410167008) and the Key Project of Bohai University (No. 0522xn078).

Conflict of interest

The authors declare no conflicts of interest.