On Some Improved Harmonic Mean Newton-Like Methods for Solving Systems of Nonlinear Equations

In this work, we have developed a fourth order Newton-like method based on harmonic mean and its multi-step version for solving system of nonlinear equations. The new fourth order method requires evaluation of one function and two first order Fréchet derivatives for each iteration. The multi-step version requires one more function evaluation for each iteration. The proposed new scheme does not require the evaluation of second or higher order Fréchet derivatives and still reaches fourth order convergence. The multi-step version converges with order 2r+4, where r is a positive integer and r ≥ 1. We have proved that the root α is a point of attraction for a general iterative function, whereas the proposed new schemes also satisfy this result. Numerical experiments including an application to 1-D Bratu problem are given to illustrate the efficiency of the new methods. Also, the new methods are compared with some existing methods.


Introduction
An often discussed problem in many applications of science and technology is to find a real zero of a system of nonlinear equations F (x) = 0, where F (x) = (f 1 (x), f 2 (x), ..., f n (x)) T , x = (x 1 , x 2 , ..., x n ) T , f i : R n → R, ∀i = 1, 2, . . ., n and F : D ⊂ R n → R n is a smooth map and D is an open and convex set, where we assume that α = (α 1 , α 2 , ..., α n ) T is a zero of the system and x (0) = x (0) 1 , x (0) 2 , ..., x (0) n T is an initial guess sufficiently close to α.For example, problems of the above type arise while solving boundary value problems for differential equations.The differential equations are reduced to system of nonlinear equations, which are in turn solved by the familiar Newton's iteration method having convergence order two [1].The Newton method (2 nd NM) is given by Homeier [2] has proposed a third order iterative method called Harmonic Mean Newton's method for solving a single nonlinear equation.Analogous to this method [2], we consider the following extension to solve a system of nonlinear equation F (x) = 0, henceforth called as 3 rd HM : We note that ))] −1 is the average of the inverses of two Jacobians.In general, such third order methods free of second derivatives like Equation ( 2) can be used for solving system of nonlinear equations.These methods require one function evaluation and two first order Fréchet derivative evaluations.The convergence analysis of a few such methods using point of attraction theory can be found in [3].This 3 rd HM method is more efficient than Halley's method because it does not require the evaluation of a third order tensor of n 3 values while finding the number of function evaluations.
Furthermore, the 3 rd HM methods are less efficient than two-step fourth order Newton's method (4 th N R) which was recently rediscovered by Noor et al. [4] using the variational iteration technique.Recently Sharma et al. [5] developed the fourth order method, which is given by Cordero et al. [6] presented a sixth order method, which is given by Recently, an improved fourth order version from a third order method for solving a single nonlinear equation is found in [7].In the current paper, similar to the method found in [7], a multivariate version having fourth order convergence is proposed.The rest of this paper is organized as follows.In Section 2, we present a new algorithm (optimal) that has fourth order convergence by using only three function evaluations and a multi-step version with order 2r + 4, where r is a positive integer and r ≥ 1 for solving systems of nonlinear equations.
In Section 3, we study the convergence analysis of the new methods using the point of attraction theory.Section 4 presents numerical examples and comparison with some existing methods.Furthermore, we also study an application problem, i.e., the 1-D Bratu problem [8].A brief conclusion is given in Section 5.

Development of the Methods
Babajee [7] has recently improved the 3 rd HM method to get a fourth order method for single equation This method is one of the member in the family of higher order multi-point iterative methods based on power mean for solving single nonlinear equation by Babajee et al. [9].
We next extend the above idea to the multivariate case.For the method given in Equation ( 2), we propose an improved fourth order Harmonic Mean Newton's method (4 th HM ) for solving systems of nonlinear equations as follows: where I is the n × n identity matrix.We further improve the 4 th HM method by additional function evaluations to get a multi-step version called (2r + 4) th HM method given by Note that this multi-step version has order 2r + 4, where r is a positive integer and r ≥ 0. The case r = 0 is the 4 th HM method.

Convergence Analysis
The main theorem is going to be demonstrated by means of the n-dimensional Taylor expansion of the functions involved.In the following, we use certain notations and results found in [10]: which lies in a neighborhood of a solution α of the nonlinear system F (x) = 0, Taylor's expansion can be applied (assuming Jacobian matrix F (α) is nonsingular) to obtain where ).Also, we can expand F (α + h) in Taylor series where I is the identity matrix.It is also noted that qC q h q−1 ∈ L(R n ).Denote e (k) = x (k) − α, so the error at the In order to prove the convergence order for the Equation ( 6), we need to recall some important definitions and results from the theory of point of attraction.
if there is an open neighborhood S of α defined by such that S ⊂ D and, for any x (0) ∈ S, the iterating {x (k) } defined by Equation (10) all lie in D and converge to α.
Theorem 1 (Ostrowski Theorem).[11] Assume that G : then α is a point of attraction for We now prove a general result that shows α is a point of attraction of a general iteration function which is a solution of the system F (x) = 0. Suppose that P, Q, R : D ⊂ R n → R n are sufficiently Fréchet differentiable functionals (depending on F ) at each point in D with P (α) = α, Q(α) = 0 and R(α) = 0.Then, there exists a ball on which the mapping Since P (x) is differentiable in α and P (α) = α, we can assume that δ was chosen sufficiently small such that for all x ∈ S with > 0 depending on δ and = 0 in case P (x) = x.
Since P , Q and R are continuously differentiable functions, then Q , R and R are bounded: Now by mean value theorem for integrals and Combining, we have which shows that G(x) is differentiable in α since δ and are arbitrary and is a solution of the system F (x) = 0. Let us suppose that x ∈ S = S(α, δ) and F (x) is continuous and nonsingular in α , and x (0) is close enough to α.Then α is a point of attraction of the sequence {x (k) } obtained using the iterative expression Equation ( 6).Furthermore, the sequence converges to α with order 4, where the error equation obtained is Proof: We first show that α is a point of attraction using Theorem 2. In this case, Now, since F (α) = 0, we have so that ρ(G (α)) = 0 < 1 and by Ostrowski's theorem, α is a point of attraction of Equation ( 6).
We next establish the fourth order convergence of this method.From Equations ( 8) and ( 9) we obtain and where e (k) = x (k) − α.
We have where and the expression for y(x (k) ) is 2 )e (k) 4 + O(e (k) 5 ).
The Taylor expansion of the Jacobian matrix F (y(x (k) )) is Therefore, +O(e (k) 4 ) and then Also, where (16) On the other hand, using Equations ( 14) and ( 16), the harmonic mean can be expressed as Using Equations ( 15) and ( 17), we have + O(e (k) 4 ) (18) Finally, by using Equations ( 13) and (18) in Equation ( 6) with some simplifications, the error equation can be expressed as: Thus from Equation (19), it can be concluded that the order of convergence of the 4 th HM method is four.
For the case r ≥ 1 we state and prove the following theorem.
Theorem 4. Let F : D ⊆ R n −→ R n be sufficiently Fréchet differentiable at each point of an open convex neighborhood D of α ∈ R n that is a solution of the system F (x) = 0. Let us suppose that x ∈ S = S(α, δ) and F (x) is continuous and nonsingular in α, and x (0) is close enough to α.Then α is a point of attraction of the sequence {x (k) } obtained using the iterative expression Equation (7).Furthermore the sequence converges to α with order 2r + 4, where r is a positive integer and r ≥ 1.

Numerical Examples
In this section, we compare the performance of the contributed Equations ( 6) and ( 7) with different methods given in Equations ( 1)-( 5).The numerical experiments have been carried out using MATLAB 7.6 software for the test problems given below.The approximate solutions are calculated correct to 1000 digits by using variable precision arithmetic.We use the following stopping criterion for the iterations: We have used the approximated computational order of convergence p c given by (see [12]) Let M be the number of iterations required for reaching the minimum residual err min .
Test Problem 2 (TP2) We consider the following system given in [3]: The solution is α ≈ (0.909569, 0.661227, 1.575834) T .We choose the starting vector x (0) = (1, 0.5, 1.5) T .The Jacobian matrix has 7 non-zero elements and it is given by Test Problem 3 (TP3) We consider the following system given in [3]: We solve this system using the initial approximation x (0) = (0.5, 0.5, 0.5, −0.2) T .The solution of this system is α ≈ (0.577350, 0.577350, 0.577350, −0.288675) T .The Jacobian matrix that has 12 non-zero elements is given by Table 1 shows the results for the test problems (TP1, TP2, TP3), from which we conclude that the 10 th HM method is the most efficient method with least number of iterations and residual error.In Table 2, we have given CPU time for the proposed methods and some existing methods.Next, we consider the (2r + 4) th HM family of methods for finding the least value of r and thus the value of p in order to get the number of iteration M = 2 and err min = 0. To achieve this, TP1 requires r = 6 (p = 16), TP2 requires r = 18 (p = 40) and TP3 requires r = 8 (p = 20).Furthermore, it is observed that the order of convergence p depends on the test problem and its starting vector.

1-D Bratu Problem
The 1-D Bratu problem [8] is given by with the boundary conditions U (0) = U (1) = 0.The 1-D planar Bratu problem has two known, bifurcated, exact solutions for values of λ < λ c , one solution for λ = λ c and no solution for λ > λ c .
The critical value of λ c is simply 8(η 2 − 1), where η is the fixed point of the hyperbolic cotangent function coth (x).The exact solution to Equation ( 26) is known and can be presented here as where θ is a constant to be determined, which satisfies the boundary conditions and is carefully chosen and assumed to be the solution of the differential Equation (26).Using a similar procedure as in [14], we show how to obtain the critical value of λ.Substitute Equation (27) in Equation ( 26), simplify and collocate at the point x = 1 2 because it is the midpoint of the interval.Another point could be chosen, but low order approximations are likely to be better if the collocation points are distributed somewhat evenly throughout the region.Then, we have Differentiating Equation ( 28) with respect to θ and setting dλ dθ = 0, the critical value λ c satisfies By eliminating λ from Equations ( 28) and (29), we have the value of θ c for the critical λ c satisfying for which θ c = 4.798714560 can be obtained using an iterative method.We then get λ c = 3.513830720 from Equation (28). Figure 1 illustrates this critical value of λ.The finite dimensional problem using standard finite difference scheme is given by with discrete boundary conditions U 0 = U N = 0 and the step size h = 1/N .There are N − 1 unknowns (n = N − 1).The Jacobian is a sparse matrix and its typical number of nonzero per row is three.It is known that the finite difference scheme converges to the lower solution of the 1-D Bratu using the starting vector U (0) = (0, 0, .., 0) T .We use N = 101 (n = 100) and test for 350 λ's in the interval (0, 3.5] (interval width = 0.01).For each λ, we let M λ be the minimum number of iterations for which U (k+1) j − U (k) j 2 < 1e − 13, where the approximation U (k) j is calculated correct to 14 decimal places.Let M λ be the mean of iteration number for the 350 λ's. Figure 2 and Table 3 give the results for the 1-D Bratu problem, where M represents number of iterations for convergence.It can be observed from the six methods considered in Table 3 that as λ increases to its critical value, the number of iterations required for convergence increase.However, as the order of method increases, the mean of iteration number decreases.The 6 th HM is the most efficient method among the six methods because it has the lowest mean iteration number and the highest number of λ converging in 2 iterations.For each λ, we find the minimum order of the (2r + 4) th HM family so that we reach convergence in 2 iterations and the results are shown in Figure 3.It can be observed that as the value of λ increases, the value of p required for convergence in 2 iterations also increases.For λ ∈ [0.01, 0.04], we require p = 4 (4 th HM ).For λ ∈ [0.05, 0.35], we require p = 6 (6 th HM ).For λ ∈ [0.36, 0.83], we require p = 8 (8 th HM ).For λ ∈ [0.84, 1.29], we require p = 10 (10 th HM ).For λ ∈ [1.30, 1.66], we require p = 12 (12 th HM ).For λ ∈ [1.66, 1.95], we require p = 14 (14 th HM ).For λ ∈ [1.96, 2.19], we require p = 16 (16 th HM ).For λ ∈ [2.20, 2.37], we require p = 18 (18 th HM ).For λ ∈ [2.38, 2.52], we require p = 20 (20 th HM ).For λ ∈ [2.53, 2.64], we require p = 22 (22 th HM ) and so on.We notice that the width of the interval decrease and the order of the family is very high as λ tends to its critical value.Finally, for λ = 3.5, we require p = 260 to reach convergence in 2 iterations.

Conclusion
In this work, we have proposed a fourth order method and its multi-step version having higher order convergence using weight functions to solve systems of nonlinear equations.The proposed schemes do not require the evaluation of second or higher order Fréchet derivatives to reach fourth order or higher order of convergence.We have tested a few examples using the proposed schemes and compared them with some known schemes, which illustrate the superiority of the new schemes.Finally, the proposed new methods have been applied on a practical problem called the 1-D Bratu problem.The results obtained are interesting and encouraging for the new methods.Hence, the proposed methods can be considered competent enough to some of the existing methods.

Figure 1 .
Figure 1.Variation of θ for different values of λ.

Figure 2 .
Figure 2. Variation of number of iteration with λ for the 2 nd N M , 3 rd HM , 4 th HM and 6 th HM methods.

Table 1 .
Comparison of different methods for system of nonlinear equations.

Table 2 .
Comparison of CPU time (s).

Table 3 .
Comparison of number of λ's in different methods for 1-D Bratu problem.