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

Authors: Diyashvir K. R. Babajee; Alicia Cordero; Fazlollah Soleymani; Juan R. Torregrosa
Journal: Journal of Applied Mathematics
 (2012)

---

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


---

## Body

## 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 F(x) at an iterate point is singular or almost singular, namely |F′(x)|≈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:
(1.1)y(k)=x(k)-[F′(x(k))-diag(β1f1,β2f2,…,βnfn)(x(k))]-1F(x(k)),x(k+1)=y(k)-[F′(y(k))-diag(β1f1,β2f2,…,βnfn)(y(k))]-1F(y(k)),
denoted by NR, where fi,i=1,…,n are the coordinate functions of F.However, the condition of|F′(x)|≠0 cannot be always avoided, though as discussed in [10], the generalized inverse could be used instead of F′(x)-1, 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 Section3 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 functionF(x)=(f1(x),f2(x),…,fn(x))T,x=(x1,x2,…,xn)T,fi:ℝn→ℝ,foralli=1,2,…,n and F:D⊂ℝn→ℝn is a smooth map and D is an open and convex set. We assume that α=(α1,α2,…,αn)T is a zero of the system, and y=(y1,y2,…,yn)T is an initial guess sufficiently close to α.Consider Taylor expansion ofF(x) about y as follows:
(2.1)F(x)=F(y)+F′(y)(x-y)+G(x),
where G(x)=O((x-y)2). Replace y by x(k) into (2.1) to get
(2.2)F(x)≈F(x(k))+F′(x(k))(x-x(k)).
Solve (2.2) for x to get
(2.3)x≈x(k)-[F′(x(k))]-1F(x(k))+[F′(x(k))]-1F(x).
If x is an approximation of α, we denote it by x(k+1) and get
(2.4)x(k+1)≈x(k)-[F′(x(k))]-1F(x(k))+[F′(x(k))]-1F(x(k+1)),
that is,
(2.5)x(k+1)-x(k)≈[F′(x(k))]-1(F(x(k+1))-F(x(k))).Waziri et al. [11] proposed a technique to approximate the Jacobian inverse [F′(x(k))]-1 by a diagonal matrix Dk=diag(di(k)), which is updated at each iteration, and where the ith diagonal entry of the matrix Dk is defined to be
(2.6)di(k)=xi(k)-xi(k-1)fi(x(k))-fi(x(k-1)),
where xi(k) is the ith component of the vector x(k), for all i=1,2,…,n and k=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 fi(x(k))-fi(x(k-1)) is not equal to zero or |fi(x(k))-fi(x(k-1))|>δ, where δ is a positive real number for all i=1,2,…,n. If else, then set di(k)=di(k-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]
(2.7)x(k+1)=G3rdAM(x(k))=x(k)-A1(x(k))-1F(x(k)),
where
(2.8)A1(x)=F′(x)+F′(y)2,
wherein y 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 n3 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

Letp be the order of a method and d be defined as
(3.1)d=d0n+∑j=1qdjnj+1,
where d0 and dj represent the number of times F and F(j) are to be evaluated, respectively. The definition of the Logarithms Informational Efficiency or Efficiency Index for nonlinear systems (IE) [1] is given by
(3.2)IE=lnpd0n+∑j=1qdjnj+1.
The efficiency indices of the Newton method (2nd  NM) and the third-order methods free from second derivatives (3rd  AM) are given by
(3.3)IE2ndNM=ln2n+n2,IE3rdAM=ln3n+2n2,
respectively. We observe that IE3rdAM>IE2ndNM, if n=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 (IE4thNR = ln4/(2n+2n2) = IE2ndNR). 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 f(x)=0. A simplified example is
(3.4)yk=xk-23f(xk)f′(xk),xk+1=xk-2f(xk)f′(xk)+f′(yk)[1-14(f′(yk)f′(xk)-1)+34(f′(yk)f′(xk)-1)2].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:(3.5)x(k+1)=G4thAM(x(k))=x(k)-H1(x(k))[A1(x(k))]-1F(x),
where H1(x)=(I-(1/4)(τ(x)-I)+(3/4)(τ(x)-I)2),τ(x)=[F′(x)]-1F′(x-(2/3)u(x)), and I is the n×n identity matrix and u(x)=[F′(x)]-1F(x). Thus, by simplifying, we can produce the following fourth-order convergent algorithm.Proposed AlgorithmStep 1.
Solve the linear systemF′(x(k))u(x(k))=F(x(k)),Step 2.
Calculate(i)
y(k)=x(k)-(2/3)u(x(k)),(ii)
τ(x(k))=[F′(x(k))]-1F′(y(k)),Step 3.
Solve the linear systemA1(x(k))h(x(k))=F(x(k)), where A1(x(k))=(1/2)(F′(x(k))+F′(y(k))),Step 4.
Calculatex(k+1)=x(k)-(I-(1/4)(τ(x(k))-I)+(3/4)(τ(x(k))-I)2)h(x(k)).The main theorem is going to be demonstrated by means of then-dimensional Taylor expansion of the functions involved. Let F:D⊆ℝn→ℝn be sufficiently Frechet differentiable in D. By using the notation introduced in [14], the qth derivative of F at u∈ℝn,q≥1, is the q-linear function F(q)(u):ℝn×⋯×ℝn→ℝn such that F(q)(u)(v1,…,vq)∈ℝn. It is easy to observe that(1)
F(q)(u)(v1,…,vq-1,·)∈ℒ(ℝn)(2)
F(q)(u)(vσ(1),…,vσ(q))=F(q)(u)(v1,…,vq), for all permutation σ of {1,2,…,q}.So, in the following, we will denote:(a)
F(q)(u)(v1,…,vq)=F(q)(u)v1⋯vq,(b)
F(q)(u)vq-1F(p)vp=F(q)(u)F(p)(u)vq+p-1.It is well known that forα+h∈ℝn lying in a neighborhood of a solution α of the nonlinear system F(x)=0, Taylor’s expansion can be applied (assuming that the Jacobian matrix F′(α) is nonsingular), and
(3.6)F(α+h)=F′(α)[h+∑q=2p-1Cqhq]+O(hp),
where Cq=(1/q!)[F′(α)]-1F(q)(α),q≥2. We observe that Cqhq∈ℝn since F(q)(α)∈ℒ(ℝn×⋯×ℝn,ℝn) and [F′(α)]-1∈ℒ(ℝn).In addition, we can expressF′ as
(3.7)F′(α+h)=F′(α)[I+∑q=2p-1qCqhq-1]+O(hp),
where I is the identity matrix. Therefore, qCqhq-1∈ℒ(ℝn). From (3.7), we obtain
(3.8)[F′(α+h)]-1=[I+X1h+X2h2+X3h3+⋯][F′(α)]-1+O(hp),
where
(3.9)X1=-2C2,X2=4C22-3C3,X3=-8C23+6C2C3+6C3C2-4C4,⋮
We denote e(k)=x(k)-α as the error in the kth iteration. The equation
(3.10)e(k+1)=Le(k)p+O(e(k)p+1),
where L is a p-linear function L∈ℒ(ℝn×⋯×ℝn,ℝn), is called the error equation, and p is the order of convergence. Observe that e(k)p is (e(k),e(k),…,e(k)).Theorem 3.1.
LetF:D⊆ℝn→ℝn be sufficiently Frechet differentiable at each point of an open convex neighborhood D of α∈ℝn, that is, a solution of the system F(x)=0. Let us suppose that F′(x) is continuous and nonsingular in α, and x(0) close enough to α. Then the sequence {x(k)}k≥0 obtained using the iterative expression (3.5) converges to α with order 4, being the error equation
(3.11)e(k+1)=(19C4+113C23-C3C2)e(k)4+O(e(k)5).Proof.
From (3.6) and (3.7), we obtain
(3.12)F(x(k))=F′(α)[e(k)+C2e(k)2+C3e(k)3+C4e(k)4]+O(e(k)5),F′(x(k))=F′(α)[I+2C2e(k)+3C3e(k)2+4C4e(k)3+5C5e(k)4]+O(e(k)5),
where Ck=(1/k!)[F′(α)]-1F(k)(α),k=2,3,…, and e(k)=x(k)-α.
From (3.8), we have
(3.13)[F′(x(k))]-1=[I+X1e(k)+X2e(k)2+X3e(k)3][F′(α)]-1+O(e(k)4),
where X1=-2C2,X2=4C22-3C3 and X3=-8C23+6C2C3+6C3C2-4C4.
Then,(3.14)[F′(x(k))]-1F(x(k))=e(k)-C2e(k)2+2(C22-C3)e(k)3+(-4C23+4C2C3+3C3C2-3C4)e(k)4+O(e(k)5),
and the expression for y(k) is
(3.15)y(k)=α+13e(k)+23C2e(k)2-43(C22-C3)e(k)3+(2C4-83C2C3-2C3C2+8C23)e(k)4+O(e(k)5).
The Taylor expansion of Jacobian matrixF′(y(k)) is
(3.16)F′(y(k))=F′(α)[I+2C2(y(k)-α)+3C3(y(k)-α)2+4C4(y(k)-α)3F′(α)k+5C5(y(k)-α)4]+O(e(k)5)=F′(α)[I+N1e(k)+N2e(k)2+N3e(k)3]+O(e(k)4),
where
(3.17)N1=23C2,N2=43C22+13C3,N3=-83C2(C22-C3)+43C3C2+427C4.
Therefore,
(3.18)τ(x(k))=[F′(x(k))]-1F′(y(k))=I+(N1+X1)e(k)+(N2+X1N1+X2)e(k)2+(N3+X1N2+X2N1+X3)e(k)3+O(e(k)4)=I-43C2e(k)+(4C22-83C3)e(k)2+(-323C23+8C2C3+163C3C2-10427C4)e(k)3+O(e(k)4),
and then,
(3.19)H1(x(k))-I=-14(τ(x(k))-I)+34(τ(x(k))-I)2=13C2e(k)+(13C22+23C3)e(k)2+(-163C23+23C2C3+43C3C2+2627C4)e(k)3+O(e(k)4).
On the other hand, the arithmetic mean of the Jacobian matrices can be expressed as(3.20)A1(x(k))=12(F′(x(k))+F′(y(k)))=F′(α)[I+43C2e(k)+(23C22+53C3)e(k)2F′(α)a+(-43C23+43C2C3+23C3C2+5627C4)e(k)3]+O(e(k)4).
From an analogous reasoning as in (3.8), we obtain
(3.21)[A1(x(k))]-1=[I+Y1e(k)+Y2e(k)2+Y3e(k)3][F′(α)]-1+O(e(k)4),
where
(3.22)Y1=-43C2,Y2=109C22-53C3,Y3=-5627C4+2027C23+89C2C3+149C3C2.
So,
(3.23)h(x(k))=A1(x(k))-1F(x(k))=e(k)-13C2e(k)2+(-29C22-23C3)e(k)3+(5027C23-49C2C3-19C3C2-2927C4)e(k)4+O(e(k)5),
and finally, by using (3.19) and (3.23), the error equation can be expressed as
(3.24)e(k+1)=x(k)-H1(x(k))A1(x(k))-1F(x(k))-α=(19C4+113C23-C3C2)e(k)4+O(e(k)5).
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(3.25)IE4thAM=ln4n+2n2,IE4thAMIE2ndNM=2n+22n+1>1,n≥1.
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:
(4.1)err=∥F(x(k))∥+∥x(k+1)-x(k)∥<1·e-100.
We have used the approximated computational order of convergence p given by (see [15])
(4.2)pc≈log(∥x(k+1)-x(k)∥)/(∥x(k)-x(k-1)∥)log(∥x(k)-x(k-1)∥)/(∥x(k-1)-x(k-2)∥).
Let Itermax be the number of iterations required before convergence is reached and errmin be the minimum residual.We check the mentioned method by solving the following test problems.Test Problem 1 (TP1). Consider F(x1,x2)=0, where F:(4,6)×(5,7)→ℝ2, and
(4.3)F(x1,x2)=(x12-x2-19,x236-x12+x2-17).
The Jacobian matrix is given by F′(x)=(2x1-1-2x1(1/2)x22+1). The starting vector is x(0)=(5.1,6.1)T and the exact solution is α=(5,6)T. In this case, the comparison of efficiencies is
(4.4)IE2ndNM=ln22+4=0.1155=IE4thNR,IE3rdAM=ln32+2(4)=0.1099,IE4thAM=ln42+2(4)=0.1386.Test Problem 2 (TP2). Consider
(4.5)cosx2-sinx1=0,x3x1-1x2=0,expx1-x32=0.The solution is(4.6)α≈(0.909569,0.661227,1.575834)T.We choose the starting vectorx(0)=(1,0.5,1.5)T. The Jacobian matrix is given by
(4.7)F′(x)=(-cosx1-sinx20x3x1lnx31x22x3x1x1x3expx10-2x3),
and has 7 nonzero elements. In this case, the comparison of efficiencies is: IE2ndNM=ln2/(3+7) = 0.0693 = IE4thNR,IE3rdAM = ln3/(3+2(7)) = 0.0646,IE4thAM = ln4/(3+2(7)) = 0.0815.Test Problem 3 (TP3). Consider
(4.8)x2x3+x4(x2+x3)=0,x1x3+x4(x1+x3)=0,x1x2+x4(x1+x2)=0,x1x2+x1x3+x2x3=1.
We solve this system using the initial approximation x(0)=(0.5,0.5,0.5,-0.2)T. The solution is
(4.9)α≈(0.577350,0.577350,0.577350,-0.288675)T.
The Jacobian is a matrix given by
(4.10)F′(x)=(0x3+x4x2+x4x2+x3x3+x40x1+x4x1+x3x2+x4x1+x40x1+x2x2+x3x1+x3x1+x20),
with 12 nonzero elements. In this case, the comparison of efficiencies is: IE2ndNM=ln2/(4+12) = 0.0433 = IE4thNR,IE3rdAM = ln3/(4+2(12)) = 0.0392,IE4thAM = ln4/(2+2(12)) = 0.0495.Table1 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.
Methods
TP1
TP2
TP3
Iter
max
err
min
p
c
Iter
max
err
min
p
c
Iter
max
err
min
p
c
2nd NM
7
8.8e–113
2
9
1.7e–107
2
8
1.7e–144
2.02
3rd AM
5
1.9e–140
3.02
7
5.9e–255
3
6
6.3e–276
3.02
4th NR
4
8.8e–113
4.02
5
1.7e–107
4
5
8.9e–291
4.11
4th AM
4
4.9e–102
4.02
6
2.7e–299
4
5
1.0e–237
4.14


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



---

*Source: 165452-2012-11-06.xml*