On a Novel Fourth-Order Algorithm for Solving Systems of Nonlinear Equations

Abstract

This paper focuses on solving systems of nonlinear equations numerically. We propose an efficient iterative scheme including two steps and fourth order of convergence. The proposed method does not require the evaluation of second or higher order Frechet derivatives per iteration to proceed and reach fourth order of convergence. Finally, numerical results illustrate the efficiency of the method.

1. Introduction

Many relationships in nature are inherently nonlinear, which according to these effects are not in direct proportion to their cause. In fact, a large number of such real-world applications are reduced to solve a system of nonlinear equations numerically. Solving such systems has become one of the most challenging problems in numerical analysis, (see [1–4]), for example in Chaos (chaos theory studies the behavior of dynamical systems that are highly sensitive to initial conditions), or in climate studies (the weather is famously chaotic, where simple changes in one part of the system produce complex effects throughout).

Some robust and efficient methods for solving nonlinear systems are brought forward (e.g., see [5–7]). The Newton’s iteration for nonlinear equations is a basic method, which converges quadratically under some conditions; this method has also been extended for systems of nonlinear equations. In this method, if the derivative of the multivariate function 𝐹(𝑥) at an iterate point is singular or almost singular, namely |𝐹^′(𝑥)|≈0, the iteration can be aborted due to the overflow or leads to divergence, which restrict its applications to some extent. Recently, there has been some modified Newton’s methods with cubic convergence for systems of nonlinear equations in order to overcome on this flaw; see for example [8]. In fact, there exist fourth-order schemes that avoid the condition of nonsingular Jacobian matrix. For example, Noor et al. in [9] design different algorithms for solving nonlinear systems of equations, obtained by means of variational techniques. One of the iterative schemes proposed in this paper, that we use in the numerical section, has the following expression: 𝑦^(𝑘)=𝑥^(𝑘)−𝐹^′𝑥^(𝑘)−diag𝛽₁𝑓₁,𝛽₂𝑓₂,…,𝛽_𝑛𝑓_𝑛𝑥^(𝑘)−1𝐹𝑥^(𝑘),𝑥^(𝑘+1)=𝑦^(𝑘)−𝐹^′𝑦^(𝑘)−diag𝛽₁𝑓₁,𝛽₂𝑓₂,…,𝛽_𝑛𝑓_𝑛𝑦^(𝑘)−1𝐹𝑦^(𝑘), (1.1) denoted by NR, where 𝑓_𝑖,𝑖=1,…,𝑛 are the coordinate functions of 𝐹.

However, the condition of |𝐹^′(𝑥)|≠0 cannot be always avoided, though as discussed in [10], the generalized inverse could be used instead of 𝐹^′(𝑥)⁻¹, and the convergence still remains quadratically. Moreover, the efficiency index of these solvers is not satisfactory for practical problems, which provide us a way to suggest high order iterative methods with high efficiency index.

In this paper, we suggest a novel two-step iterative scheme for finding the solution of nonlinear systems whose Jacobian matrix is required to be nonsingular. The main merit of our proposed scheme is that per iteration it does not require the evaluation of the second, or higher order Frechet derivatives, while the fourth-order of convergence is attainable using only one function and two first-order Frechet derivatives per computing step.

The paper is organized as follows. In the next section, we give a short review on basic notes in this field of study, while Section 3 proposes the new iterative algorithm to solve systems of nonlinear equations and analyze its local convergence. In Section 4, numerical examples will be further considered to show that the new method performs well. The detailed results for numerical examples are analyzed and discussed in this section. And finally, some conclusions are drawn in Section 5.

2. Preliminary Notes

Let us consider the problem of finding a real zero of a multivariate nonlinear function 𝐹(𝑥)=(𝑓₁(𝑥),𝑓₂(𝑥),…,𝑓_𝑛(𝑥))^𝑇,𝑥=(𝑥₁,𝑥₂,…,𝑥_𝑛)^𝑇,𝑓_𝑖∶ℝ^𝑛→ℝ,forall𝑖=1,2,…,𝑛 and 𝐹∶𝐷⊂ℝ^𝑛→ℝ^𝑛 is a smooth map and 𝐷 is an open and convex set. We assume that 𝛼=(𝛼₁,𝛼₂,…,𝛼_𝑛)^𝑇 is a zero of the system, and 𝑦=(𝑦₁,𝑦₂,…,𝑦_𝑛)^𝑇 is an initial guess sufficiently close to 𝛼.

Consider Taylor expansion of 𝐹(𝑥) about 𝑦 as follows: 𝐹(𝑥)=𝐹𝑦+𝐹^′𝑦𝑥−𝑦+𝐺(𝑥), (2.1) where 𝐺(𝑥)=𝑂((𝑥−𝑦)²). Replace 𝑦 by 𝑥^(𝑘) into (2.1) to get 𝐹(𝑥)≈𝐹𝑥^(𝑘)+𝐹^′𝑥^(𝑘)𝑥−𝑥^(𝑘). (2.2) Solve (2.2) for 𝑥 to get 𝑥≈𝑥^(𝑘)−𝐹^′𝑥^(𝑘)−1𝐹𝑥^(𝑘)+𝐹^′𝑥^(𝑘)−1𝐹(𝑥). (2.3) If 𝑥 is an approximation of 𝛼, we denote it by 𝑥^(𝑘+1) and get 𝑥^(𝑘+1)≈𝑥^(𝑘)−𝐹^′𝑥^(𝑘)−1𝐹𝑥^(𝑘)+𝐹^′𝑥^(𝑘)−1𝐹𝑥^(𝑘+1), (2.4) that is, 𝑥^(𝑘+1)−𝑥^(𝑘)≈𝐹^′𝑥^(𝑘)−1𝐹𝑥^(𝑘+1)−𝐹𝑥^(𝑘). (2.5)

Waziri et al. [11] proposed a technique to approximate the Jacobian inverse [𝐹^′(𝑥^(𝑘))]⁻¹ by a diagonal matrix 𝐷_𝑘=diag(𝑑_𝑖^(𝑘)), which is updated at each iteration, and where the 𝑖th diagonal entry of the matrix 𝐷_𝑘 is defined to be 𝑑_𝑖^(𝑘)=𝑥_𝑖^(𝑘)−𝑥_𝑖^(𝑘−1)𝑓_𝑖𝑥^(𝑘)−𝑓_𝑖𝑥^(𝑘−1), (2.6) where 𝑥_𝑖^(𝑘) is the 𝑖th component of the vector 𝑥^(𝑘), for all 𝑖=1,2,…,𝑛 and 𝑘=0,1,2,…. We here remark that this procedure is similar to the Secant method in one-dimensional case. The approximation is valid only if 𝑓_𝑖(𝑥^(𝑘))−𝑓_𝑖(𝑥^(𝑘−1)) is not equal to zero or |𝑓_𝑖(𝑥^(𝑘))−𝑓_𝑖(𝑥^(𝑘−1))|>𝛿, where 𝛿 is a positive real number for all 𝑖=1,2,…,𝑛. If else, then set 𝑑_𝑖^(𝑘)=𝑑_𝑖^(𝑘−1), for more details, please refer to [11].

The order of convergence and efficiency index of such solvers are indeed not good for practical problems in which iterations with high efficiencies are required. This was just the motivation of producing higher-order iterative nonlinear systems solvers. Third-order methods free from second derivatives were proposed from Quadrature rules for solving systems of nonlinear equations. These methods require one functional evaluation and two evaluations of first derivatives at two different points. An example of such iterations is the Arithmetic Mean Newton method (3rd AM) derived from Trapezoidal rule [12] 𝑥^(𝑘+1)=𝐺_3rdAM𝑥^(𝑘)=𝑥^(𝑘)−𝐴₁𝑥^(𝑘)−1𝐹𝑥^(𝑘), (2.7) where 𝐴₁(𝑥)=𝐹^′(𝑥)+𝐹^′𝑦2, (2.8) wherein 𝑦 is the Newton’s iteration.

The convergence analysis of these methods using point of attraction theory can be found in [13]. This third-order Newton-like method is more efficient than Halley’s method because it does not require the evaluation of a third-order tensor of 𝑛³ values.

In the next section, we present our new algorithm which reaches the highest possible order four by using only three functional evaluations of the multivariate function for solving systems of nonlinear equations.

3. A Novel Computational Iterative Algorithm

Let 𝑝 be the order of a method and 𝑑 be defined as 𝑑=𝑑₀𝑛+^𝑞_𝑗=1𝑑_𝑗𝑛^𝑗+1, (3.1) where 𝑑₀ and 𝑑_𝑗 represent the number of times 𝐹 and 𝐹^(𝑗) are to be evaluated, respectively. The definition of the Logarithms Informational Efficiency or Efficiency Index for nonlinear systems (IE) [1] is given by IE=ln𝑝𝑑₀𝑛+∑^𝑞_𝑗=1𝑑_𝑗𝑛^𝑗+1. (3.2) The efficiency indices of the Newton method (2nd  NM) and the third-order methods free from second derivatives (3rd  AM) are given by IE_2ndNM=ln2𝑛+𝑛²,IE_3rdAM=ln3𝑛+2𝑛², (3.3) respectively. We observe that IE_3rdAM>IE_2ndNM, if 𝑛=1. That is, the third-order methods free from second derivatives are less efficient than Newton’s method for systems of nonlinear equations. Let us also note that NR scheme has the same efficiency index as Newton’s method (IE_4thNR = ln4/(2𝑛+2𝑛²) = IE_2ndNR). Thus, it is important to develop fourth-order methods from these third-order methods to improve the efficiency. Soleymani et al. [5] have recently improved the 3rd  AM method to get a class of fourth-order Jarratt-type methods for solving the nonlinear equation 𝑓(𝑥)=0. A simplified example is 𝑦_𝑘=𝑥_𝑘−23𝑓𝑥_𝑘𝑓^′𝑥_𝑘,𝑥_𝑘+1=𝑥_𝑘−2𝑓𝑥_𝑘𝑓^′𝑥_𝑘+𝑓^′𝑦_𝑘11−4𝑓^′𝑦_𝑘𝑓^′𝑥_𝑘+3−14𝑓^′𝑦_𝑘𝑓^′𝑥_𝑘−1². (3.4)

We extend this method to the multivariate case. The improved fourth-order Arithmetic Mean Newton method (4th  AM) for systems of nonlinear equations can now be suggested as follows: 𝑥^(𝑘+1)=𝐺_4thAM𝑥^(𝑘)=𝑥^(𝑘)−𝐻₁𝑥^(𝑘)𝐴₁𝑥^(𝑘)−1𝐹(𝑥), (3.5) where 𝐻₁(𝑥)=(𝐼−(1/4)(𝜏(𝑥)−𝐼)+(3/4)(𝜏(𝑥)−𝐼)²),𝜏(𝑥)=[𝐹^′(𝑥)]⁻¹𝐹^′(𝑥−(2/3)𝑢(𝑥)), and 𝐼 is the 𝑛×𝑛 identity matrix and 𝑢(𝑥)=[𝐹^′(𝑥)]⁻¹𝐹(𝑥). Thus, by simplifying, we can produce the following fourth-order convergent algorithm.

Proposed Algorithm

Step 1. Solve the linear system 𝐹^′(𝑥^(𝑘))𝑢(𝑥^(𝑘))=𝐹(𝑥^(𝑘)),

Step 2. Calculate(i)𝑦^(𝑘)=𝑥^(𝑘)−(2/3)𝑢(𝑥^(𝑘)), (ii)𝜏(𝑥^(𝑘))=[𝐹^′(𝑥^(𝑘))]⁻¹𝐹^′(𝑦^(𝑘)),

Step 3. Solve the linear system 𝐴₁(𝑥^(𝑘))ℎ(𝑥^(𝑘))=𝐹(𝑥^(𝑘)), where 𝐴₁(𝑥^(𝑘))=(1/2)(𝐹^′(𝑥^(𝑘))+𝐹^′(𝑦^(𝑘))),

Step 4. Calculate 𝑥^(𝑘+1)=𝑥^(𝑘)−(𝐼−(1/4)(𝜏(𝑥^(𝑘))−𝐼)+(3/4)(𝜏(𝑥^(𝑘))−𝐼)²)ℎ(𝑥^(𝑘)).

The main theorem is going to be demonstrated by means of the 𝑛-dimensional Taylor expansion of the functions involved. Let 𝐹∶𝐷⊆ℝ^𝑛→ℝ^𝑛 be sufficiently Frechet differentiable in 𝐷. By using the notation introduced in [14], the 𝑞th derivative of 𝐹 at 𝑢∈ℝ^𝑛,𝑞≥1, is the 𝑞-linear function 𝐹^(𝑞)(𝑢)∶ℝ^𝑛×⋯×ℝ^𝑛→ℝ^𝑛 such that 𝐹^(𝑞)(𝑢)(𝑣₁,…,𝑣_𝑞)∈ℝ^𝑛. It is easy to observe that(1)𝐹^(𝑞)(𝑢)(𝑣₁,…,𝑣_𝑞−1,⋅)∈ℒ(ℝ^𝑛)(2)𝐹^(𝑞)(𝑢)(𝑣_𝜎(1),…,𝑣_𝜎(𝑞))=𝐹^(𝑞)(𝑢)(𝑣₁,…,𝑣_𝑞), for all permutation 𝜎 of {1,2,…,𝑞}.So, in the following, we will denote:(a)𝐹^(𝑞)(𝑢)(𝑣₁,…,𝑣_𝑞)=𝐹^(𝑞)(𝑢)𝑣₁⋯𝑣_𝑞, (b)𝐹^(𝑞)(𝑢)𝑣^𝑞−1𝐹^(𝑝)𝑣^𝑝=𝐹^(𝑞)(𝑢)𝐹^(𝑝)(𝑢)𝑣^{𝑞+𝑝−1}.

It is well known that for 𝛼+ℎ∈ℝ^𝑛 lying in a neighborhood of a solution 𝛼 of the nonlinear system 𝐹(𝑥)=0, Taylor’s expansion can be applied (assuming that the Jacobian matrix 𝐹^′(𝛼) is nonsingular), and 𝐹(𝛼+ℎ)=𝐹^′(𝛼)ℎ+^𝑝−1_𝑞=2𝐶_𝑞ℎ^𝑞+𝑂ℎ^𝑝, (3.6) where 𝐶_𝑞=(1/𝑞!)[𝐹^′(𝛼)]⁻¹𝐹^(𝑞)(𝛼),𝑞≥2. We observe that 𝐶_𝑞ℎ^𝑞∈ℝ^𝑛 since 𝐹^(𝑞)(𝛼)∈ℒ(ℝ^𝑛×⋯×ℝ^𝑛,ℝ^𝑛) and [𝐹^′(𝛼)]⁻¹∈ℒ(ℝ^𝑛).

In addition, we can express 𝐹^′ as 𝐹^′(𝛼+ℎ)=𝐹^′(𝛼)𝐼+^𝑝−1_𝑞=2𝑞𝐶_𝑞ℎ^𝑞−1+𝑂ℎ^𝑝, (3.7) where 𝐼 is the identity matrix. Therefore, 𝑞𝐶_𝑞ℎ^𝑞−1∈ℒ(ℝ^𝑛). From (3.7), we obtain 𝐹^′(𝛼+ℎ)⁻¹=𝐼+𝑋₁ℎ+𝑋₂ℎ²+𝑋₃ℎ³+⋯𝐹^′(𝛼)⁻¹+𝑂ℎ^𝑝, (3.8) where 𝑋₁=−2𝐶₂,𝑋₂=4𝐶²₂−3𝐶₃,𝑋₃=−8𝐶³₂+6𝐶₂𝐶₃+6𝐶₃𝐶₂−4𝐶₄,⋮ (3.9) We denote 𝑒^(𝑘)=𝑥^(𝑘)−𝛼 as the error in the 𝑘th iteration. The equation 𝑒^(𝑘+1)=𝐿𝑒^𝑝(𝑘)+𝑂𝑒^(𝑘)𝑝+1, (3.10) where 𝐿 is a 𝑝-linear function 𝐿∈ℒ(ℝ^𝑛×⋯×ℝ^𝑛,ℝ^𝑛), is called the error equation, and 𝑝 is the order of convergence. Observe that 𝑒^𝑝(𝑘) is (𝑒^(𝑘),𝑒^(𝑘),…,𝑒^(𝑘)).

Theorem 3.1. Let 𝐹∶𝐷⊆ℝ^𝑛→ℝ^𝑛 be sufficiently Frechet differentiable at each point of an open convex neighborhood 𝐷 of 𝛼∈ℝ^𝑛, that is, a solution of the system 𝐹(𝑥)=0. Let us suppose that 𝐹^′(𝑥) is continuous and nonsingular in 𝛼, and 𝑥⁽⁰⁾ close enough to 𝛼. Then the sequence {𝑥^(𝑘)}_𝑘≥0 obtained using the iterative expression (3.5) converges to 𝛼 with order 4, being the error equation 𝑒^(𝑘+1)=19𝐶₄+113𝐶³₂−𝐶₃𝐶₂𝑒^4(𝑘)+𝑂𝑒^5(𝑘). (3.11)

Proof. From (3.6) and (3.7), we obtain 𝐹𝑥^(𝑘)=𝐹^′(𝛼)𝑒^(𝑘)+𝐶₂𝑒^2(𝑘)+𝐶₃𝑒^3(𝑘)+𝐶₄𝑒^4(𝑘)+𝑂𝑒^5(𝑘),𝐹^′𝑥^(𝑘)=𝐹^′(𝛼)𝐼+2𝐶₂𝑒^(𝑘)+3𝐶₃𝑒^2(𝑘)+4𝐶₄𝑒^3(𝑘)+5𝐶₅𝑒^4(𝑘)+𝑂𝑒^5(𝑘), (3.12) where 𝐶_𝑘=(1/𝑘!)[𝐹^′(𝛼)]⁻¹𝐹^(𝑘)(𝛼),𝑘=2,3,…, and 𝑒^(𝑘)=𝑥^(𝑘)−𝛼.
From (3.8), we have 𝐹^′𝑥^(𝑘)−1=𝐼+𝑋₁𝑒^(𝑘)+𝑋₂𝑒^2(𝑘)+𝑋₃𝑒^3(𝑘)𝐹^′(𝛼)⁻¹+𝑂𝑒^4(𝑘), (3.13) where 𝑋₁=−2𝐶₂,𝑋₂=4𝐶²₂−3𝐶₃ and 𝑋₃=−8𝐶³₂+6𝐶₂𝐶₃+6𝐶₃𝐶₂−4𝐶₄.
Then, 𝐹^′𝑥^(𝑘)−1𝐹𝑥^(𝑘)=𝑒^(𝑘)−𝐶₂𝑒^2(𝑘)+2𝐶²₂−𝐶₃𝑒^3(𝑘)+−4𝐶³₂+4𝐶₂𝐶₃+3𝐶₃𝐶₂−3𝐶₄𝑒^4(𝑘)+𝑂𝑒^5(𝑘), (3.14) and the expression for 𝑦^(𝑘) is 𝑦^(𝑘)1=𝛼+3𝑒^(𝑘)+23𝐶₂𝑒^2(𝑘)−43𝐶²₂−𝐶₃𝑒^3(𝑘)+2𝐶₄−83𝐶₂𝐶₃−2𝐶₃𝐶₂+8𝐶³₂𝑒^4(𝑘)+𝑂𝑒^5(𝑘). (3.15)
The Taylor expansion of Jacobian matrix 𝐹^′(𝑦^(𝑘)) is 𝐹^′𝑦^(𝑘)=𝐹^′(𝛼)𝐼+2𝐶₂𝑦^(𝑘)−𝛼+3𝐶₃𝑦^(𝑘)−𝛼²+4𝐶₄𝑦^(𝑘)−𝛼³+5𝐶₅𝑦^(𝑘)−𝛼⁴+𝑂𝑒^5(𝑘)=𝐹^′(𝛼)𝐼+𝑁₁𝑒^(𝑘)+𝑁₂𝑒^2(𝑘)+𝑁₃𝑒^3(𝑘)+𝑂𝑒^4(𝑘), (3.16) where 𝑁₁=23𝐶₂,𝑁₂=43𝐶²₂+13𝐶₃,𝑁₃8=−3𝐶₂𝐶²₂−𝐶₃+43𝐶₃𝐶₂+4𝐶27⁴. (3.17) Therefore, 𝜏𝑥^(𝑘)=𝐹^′𝑥^(𝑘)−1𝐹^′𝑦^(𝑘)=𝐼+𝑁₁+𝑋₁𝑒^(𝑘)+𝑁₂+𝑋₁𝑁₁+𝑋₂𝑒^2(𝑘)+𝑁₃+𝑋₁𝑁₂+𝑋₂𝑁₁+𝑋₃𝑒^3(𝑘)+𝑂𝑒^4(𝑘)4=𝐼−3𝐶₂𝑒^(𝑘)+4𝐶²₂−83𝐶₃𝑒^2(𝑘)+−323𝐶³₂+8𝐶₂𝐶₃+163𝐶₃𝐶₂−104𝐶27⁴𝑒^3(𝑘)+𝑂𝑒^4(𝑘), (3.18) and then, 𝐻₁𝑥^(𝑘)1−𝐼=−4𝜏𝑥^(𝑘)+3−𝐼4𝜏𝑥^(𝑘)−𝐼²=13𝐶₂𝑒^(𝑘)+13𝐶²₂+23𝐶₃𝑒^2(𝑘)+−163𝐶³₂+23𝐶₂𝐶₃+43𝐶₃𝐶₂+26𝐶27⁴𝑒^3(𝑘)+𝑂𝑒^4(𝑘). (3.19)
On the other hand, the arithmetic mean of the Jacobian matrices can be expressed as 𝐴₁𝑥^(𝑘)=12𝐹^′𝑥^(𝑘)+𝐹^′𝑦^(𝑘)=𝐹^′(𝛼4)𝐼+3𝐶₂𝑒^(𝑘)+23𝐶²₂+53𝐶₃𝑒^2(𝑘)+−43𝐶³₂+43𝐶₂𝐶₃+23𝐶₃𝐶₂+56𝐶27⁴𝑒^3(𝑘)+𝑂𝑒^4(𝑘). (3.20) From an analogous reasoning as in (3.8), we obtain 𝐴₁𝑥^(𝑘)−1=𝐼+𝑌₁𝑒^(𝑘)+𝑌₂𝑒^2(𝑘)+𝑌₃𝑒^3(𝑘)𝐹^′(𝛼)⁻¹+𝑂𝑒^4(𝑘), (3.21) where 𝑌₁4=−3𝐶₂,𝑌₂=109𝐶²₂−53𝐶₃,𝑌₃=−56𝐶27⁴+20𝐶27³²+89𝐶₂𝐶₃+149𝐶₃𝐶₂. (3.22) So, ℎ𝑥^(𝑘)=𝐴₁𝑥^(𝑘)−1𝐹𝑥^(𝑘)=𝑒^(𝑘)−13𝐶₂𝑒^2(𝑘)+−29𝐶²₂−23𝐶₃𝑒^3(𝑘)+50𝐶27³²−49𝐶₂𝐶₃−19𝐶₃𝐶₂−29𝐶27⁴𝑒^4(𝑘)+𝑂𝑒^5(𝑘), (3.23) and finally, by using (3.19) and (3.23), the error equation can be expressed as 𝑒^(𝑘+1)=𝑥^(𝑘)−𝐻₁𝑥^(𝑘)𝐴₁𝑥^(𝑘)−1𝐹𝑥^(𝑘)=1−𝛼9𝐶₄+113𝐶³₂−𝐶₃𝐶₂𝑒^4(𝑘)+𝑂𝑒^5(𝑘). (3.24) So, taking into account (3.24), it can be concluded that the order of convergence of the proposed method is four.

The efficiency index of the fourth-order method free from second derivative (4th  AM) is given by IE_4thAM=ln4𝑛+2𝑛²,IE_4thAMIE_2ndNM=2𝑛+22𝑛+1>1,𝑛≥1. (3.25) This shows that the 4th  AM method is more efficient than 2nd  NM and 3rd  AM methods. We next conduct some numerical experiments to compare the methods.

4. Numerical Examples

In this section, we compare the performance of contributed method with Newton’s scheme and Arithmetic Mean Newton method (2.7). The algorithms have been written in MATLAB 7.6 and tested for the examples given below. For the following test problems, the approximate solutions are calculated correct to 500 digits by using variable arithmetic precision. We use the following stopping criterium: err=𝐹𝑥^(𝑘)+𝑥^(𝑘+1)−𝑥^(𝑘)<1⋅𝑒−100. (4.1) We have used the approximated computational order of convergence 𝑝 given by (see [15]) 𝑝_𝑐≈log𝑥^(𝑘+1)−𝑥^(𝑘)/𝑥^(𝑘)−𝑥^(𝑘−1)log𝑥^(𝑘)−𝑥^(𝑘−1)/𝑥^(𝑘−1)−𝑥^(𝑘−2). (4.2) Let Iter_max be the number of iterations required before convergence is reached and err_min be the minimum residual.

We check the mentioned method by solving the following test problems.

Test Problem 1 (TP1). Consider 𝐹(𝑥₁,𝑥₂)=0, where 𝐹∶(4,6)×(5,7)→ℝ², and 𝐹𝑥₁,𝑥₂=𝑥²₁−𝑥₂𝑥−19,³²6−𝑥²₁+𝑥₂−17. (4.3) The Jacobian matrix is given by 𝐹^′(𝑥)=^2𝑥1−1_{−2𝑥1(1/2)𝑥22+1}. The starting vector is 𝑥⁽⁰⁾=(5.1,6.1)^𝑇 and the exact solution is 𝛼=(5,6)^𝑇. In this case, the comparison of efficiencies is IE_2ndNM=ln22+4=0.1155=IE_4thNR,IE_3rdAM=ln32+2(4)=0.1099,IE_4thAM=ln42+2(4)=0.1386. (4.4)

Test Problem 2 (TP2). Consider cos𝑥₂−sin𝑥₁𝑥=0,_𝑥13−1𝑥₂=0,exp𝑥₁−𝑥²₃=0. (4.5)

The solution is 𝛼≈(0.909569,0.661227,1.575834)^𝑇. (4.6)

We choose the starting vector 𝑥⁽⁰⁾=(1,0.5,1.5)^𝑇. The Jacobian matrix is given by 𝐹^′(𝑥)=−cos𝑥₁−sin𝑥₂0𝑥^𝑥1₃ln𝑥₃1𝑥²₂𝑥^𝑥1₃𝑥₁𝑥₃exp𝑥₁0−2𝑥₃, (4.7) and has 7 nonzero elements. In this case, the comparison of efficiencies is: IE_2ndNM=ln2/(3+7) = 0.0693 = IE_4thNR,IE_3rdAM = ln3/(3+2(7)) = 0.0646,IE_4thAM = ln4/(3+2(7)) = 0.0815.

Test Problem 3 (TP3). Consider 𝑥₂𝑥₃+𝑥₄𝑥₂+𝑥₃𝑥=0,₁𝑥₃+𝑥₄𝑥₁+𝑥₃𝑥=0,₁𝑥₂+𝑥₄𝑥₁+𝑥₂𝑥=0,₁𝑥₂+𝑥₁𝑥₃+𝑥₂𝑥₃=1. (4.8) We solve this system using the initial approximation 𝑥⁽⁰⁾=(0.5,0.5,0.5,−0.2)^𝑇. The solution is 𝛼≈(0.577350,0.577350,0.577350,−0.288675)^𝑇. (4.9) The Jacobian is a matrix given by 𝐹^′(𝑥)=0𝑥₃+𝑥₄𝑥₂+𝑥₄𝑥₂+𝑥₃𝑥₃+𝑥₄0𝑥₁+𝑥₄𝑥₁+𝑥₃𝑥₂+𝑥₄𝑥₁+𝑥₄0𝑥₁+𝑥₂𝑥₂+𝑥₃𝑥₁+𝑥₃𝑥₁+𝑥₂0, (4.10) with 12 nonzero elements. In this case, the comparison of efficiencies is: IE_2ndNM=ln2/(4+12) = 0.0433 = IE_4thNR,IE_3rdAM = ln3/(4+2(12)) = 0.0392,IE_4thAM = ln4/(2+2(12)) = 0.0495.

Table 1 gives the results for Test Problems 1, 2, and 3. It is observed for all problems that the fourth-order method converges in the least number of iterations. The computational order of convergence agrees with the theory. 4th  AM gives the best results in terms of least residual, and it is the most efficient method compared to 2nd  NM and 3rd  AM. It can compete with the 4th  NR method and has the advantage of having a higher efficiency index.

Table 1 

Comparison of different methods for systems of nonlinear equations.

5. Concluding Remarks

The efficiency of the quadratically and cubically multidimensional methods is not satisfactory in most practical problems. Thus in this paper, we have extended a fourth-order method from a third-order method for systems of nonlinear equations. We have shown that the fourth-order iterative method is more efficient than the second order Newton and third-order methods. Numerical experiments have also shown that the fourth-order method is efficient.

Acknowledgments

The authors would like to thank the referees for the valuable comments and for the suggestions to improve the readability of the paper. This research was supported by Ministerio de Ciencia y Tecnología MTM2011-28636-C02-02.