Quadrature based Broyden-like method for systems of nonlinear equations Idowu

A new iterative method based on the quasi-Newton approach for solving systems of nonlinear equations, especially large scale is proposed. We used the weighted combination of the Trapezoidal and Simpson quadrature rules. Our goal is to enhance the efficiency of the well known Broyden method by reducing the number of iterations it takes to reach a solution. Local convergence analysis and computational results are given.


Introduction
In vector notation form, consider the system of nonlinear equations where the function F : R n → R n , is continuously differentiable in an open neighbourhood Ω ⊂ R n .One of the prominent methods used to solve (1) is Newton's method as a result of its convergence property.However, the method has some challenges such as computation and storage of Jacobian matrix at every iteration as well as inefficiency in handling large scale systems.In an attempt to overcome these challenges, quasi-Newton methods have been introduced.Among the successful variants of quasi-Newton method is the Broyden's method that generate a reasonable approximation of the Jacobian matrix with no additional evaluation of the function, (see [15] and references therein).Broyden's method is given by where matrix B p is the approximation of F (x p ), such that the quasi-Newton equation is satisfied for each p.On the other hand, Broyden method needs n 2 (n is the length of vector x) storage locations, therefore, for large scale systems, this might lead to severe memory constraints [14].
Over the years, development of algorithms cheaper than Newton's method and some of its variants has been an area of intense research.Among others, several modifications of Newton's method and its variants were proposed to reduce computational cost, accelerate the convergence and reduce evaluations of functions in each step of the 131 iterations (see [3][4][5]7,17,22] and references therein).Just like quadrature rules have been used to propose Newtonlike methods [17,[19][20][21][22], some efforts at using quadrature formulas to propose Broyden-like methods includes those of [13,14], all of which are directed towards the development of cost effective approaches.Based on this fact and frequent occurrence of systems of nonlinear equations in science and engineering, it is interesting to present an approach which will improve further the Jacobian approximation using weighted combination of Trapezoidal and Simpson (TS) quadrature formulas.
The paper is organized as follows: Section 2 is for the development of the proposed method; in Section 3, the convergence result is discussed; the numerical results are shown in Section 4; and Section 5 is for the conclusion.

Development of the Method
Let x * be a root of the nonlinear equation F (x) = 0, where F is sufficiently differentiable.Newton's method originates from the Taylor's series expansion of the function (of a single variable) f (x) about the point x 1 : where f and its first and second derivatives, f ′ and f ′′ , are calculated at x 1 .For a multiple variable function F , from the above equation, it is obvious that The matrix of partial derivatives appearing in (4) is the Jacobian J, where ∫ x xp F ′ (t)dt is multiple integrals as follows: The most obvious approach is to treat the multiple integral as a nested sequence of one-dimensional integrals, and to use one-dimensional quadrature with respect to each argument in turn [11].So we can approximate ∫ x xp F ′ (t)dt with the weight combination of the Trapezoidal and Simpson quadrature rules.The authors ( [9], [10], [14], etc.) and references therein have proposed various methods by the approximation of the indefinite integral in (4) using Newton Cotes formulae of order zero and one.Approximating the integral in (4) by the average of Trapezoidal and Simpson (TS) quadrature rules yields: where ) With the above formulation, we have the following two-step iterative schemes for solving the nonlinear system (1).

TS METHOD
For a given x 0 using initial matrix B 0 = I, compute the approximate solution x p+1 by the iterative schemes ) 132 In the following, we recall the most common notions and results about the convergence of an iterative method needed for subsequent development.

Local Convergence Result
Our aim is to prove that the proposed TS method converges superlinearly.To achieve this, we consider the following standard assumptions on the nonlinear function (ii) There exists x * ∈ Ω such that F (x * ) = 0 and F (x * ) is nonsingular and continuous for every x; (iii) F ′ (x) is Lipschitz continuous and hence satisfies the Lipschitz condition of order 1 such that there exists a positive constant λ such that The following is the main result, which is a modified result of [14].

Theorem 1
Let F : R n → R n satisfy the hypothesis of Lemma 2.1 on the set Ω. Let B k be a sequence of nonsingular matrices in L(R n ), the space of real matrices of order n.Suppose for some x 0 the sequence x k generated by ( 5) remains in Ω and lim p→∞ x p = x * where for each p, x p ̸ = x * .Then {x p } converges q-superlinearly to x * and F (x * ) = 0 if and only if where Proof Suppose (6) holds.Then (5) becomes Take the norm of both sides to have: Using vector norm properties, we have Divide through by ∥s p ∥ to have Using Lemma 2.1, Since x p → x * ∀ p, then, from (6), we have Therefore, But F ′ (x * ) is nonsingular.Thus, by Lemma 2.1, ∃ ρ > 0, p 0 ≥ 0 such that we have For all p ≥ p 0 , (8) and Equation ( 9) gives Therefore, x p converges q-superlinearly to x * .Conversely, suppose that x p converges q-superlinearly to x * and F (x * ) = 0.Then, by Lemma 2.1, there exists ρ > 0 such that we have Then, Using Lemma 2.2, we have From ( 7), we have Since x p converges to x * , then lim The proof is complete.

Numerical Results
We used six test functions with eight instances of dimension n = 5 to n = 1065, which makes a total of 48 problems, in order to check the effectiveness of the proposed methods.The x 0 stands for the initial approximation to the solution in all the tested problems.

Stopping criterion
We have to stop the program if any of the following conditions are met.
(2) the number of iterations reaches 500.(We do not want the program to run indefinitely) A comparison of the numerical test results of our two-step methods is made with these well-known methods: Classical Broyden Method [2], Trapezoidal Broyden Method [13] and Midpoint-Simpson Broyden Method [14].
Table 1 shows the results on the number of iterations (NI) and the time (CPU) needed to converge to a solution.A failure is reported (denoted by '-') in the tabulated results.All the methods were implemented using MATLAB R2012b.All computations were carried out on a PC with Intel(R) Pentium(R) processor with 4GB of RAM and CPU 2.20 GHz.The comparison of the tested methods is based on the performance profile of [8] and the comparison indices of [1].

Set of Test Problems
Problem 2 [16] Problem 3 [16] The Robustness index measures the performance of the algorithm on a wide variety of problems in their class for all reasonable values of the starting point.The Efficiency index records the use of low amount of overall resources involved in the programming, while the combined Robustness and Efficiency index measures both robustness and efficiency performances.The Robustness (R), Efficiency (E) and Combined Robustness and Efficiency (C ER ) performance of each solver are shown as proposed by [1].

Discussion
The Tables show the results obtained for the benchmark problems.The Figures show the performance profile of each of the methods for small value of τ , the best possible ratio.We assume that we have n s solvers and n p problems, the performance profile P : R → [0; 1] is defined as follows: let P and S be the set of problems and set of solvers respectively.For each problem p ∈ P and for each solver s ∈ S, we define t p;s := (number of iterations required to solve problem p by solver s) and t p;s := (computing time required to solve problem p by solver s) for Figures 1 and 2 respectively.The performance ratio is given by r p;s = tp;s min tp;s,s∈S .Then the performance profile is defined by P (τ ) := 1 np sizep ∈ P : log 2 (r p,s ) ≤ τ , ∀τ ∈ R, where P (τ ) is the probability for solver s ∈ S that a performance ratio r p,s is within a factor τ ∈ R of the best possible ratio.
In both Figures, TS climb off the graph and this indicates that it solves all of the test problems successfully.MSB has the highest probability of 0.83 and 0.57 in Figures 1 and 2 respectively followed by the proposed method with 0.27 and 0.22, TB with 0.17 and 0.18 and CB with 0.02 and 0.04.In Figure 1, TS competes with MSB for τ ≥ 0.4.For factor τ ≥ 0.3, MSB and TS have the best probabilities in Figure 2. On the set of benchmark problems, TS is superior in terms of robustness over other methods.The number of iterations obtained is promising and competitive, the CPU time is also encouraging.The results of the proposed method are better than one of CB and TB.However, MSB is more efficient than the proposed method.In conclusion, the numerical tests have shown that our new method is comparable with the well known existing Broyden-like methods.

Conclusion
The Broyden method was introduced to overcome the shortcomings of the Newton's method but still requires a higher number of iterations than the Newton's method.The Broyden-like method proposed is a promising and competitive alternative to quasi-Newton and Broyden method, especially for large scale systems.

Table 1 .
Comparison of Numerical Results

Table 2 .
Summary of Robustness, Efficiency and Combined Robustness and Efficiency Measures