Choosing Improved Initial Values for Polynomial Zerofinding in Extended Newbery Method to Obtain Convergence

Abstract

In all polynomial zerofinding algorithms, a good convergence requires a very good initial approximation of the exact roots. The objective of the work is to study the conditions for determining the initial approximations for an iterative matrix zerofinding method. The investigation is based on the Newbery's matrix construction which is similar to Fiedler's construction associated with a characteristic polynomial. To ensure that convergence to both the real and complex roots of polynomials can be attained, three methods are employed. It is found that the initial values for the Fiedler's companion matrix which is supplied by the Schmeisser's method give a better approximation to the solution in comparison to when working on these values using the Schmeisser's construction towards finding the solutions. In addition, empirical results suggest that a good convergence can still be attained when an initial approximation for the polynomial root is selected away from its real value while other approximations should be sufficiently close to their real values. Tables and figures on the errors that resulted from the implementation of the method are also given.

1. Introduction

In recent years, various researches have been studied on the zerofinding algorithms. For the first time, Galois established that a general direct method for calculating zeroes in terms of explicit formulas exists only for general polynomials of degree less than five. Thus finding the polynomial roots with higher degree needs numerical methods and each algorithm possesses its own advantages and disadvantages. Wilkinson [1, 2] pointed out that there is no general zerofinding algorithm that can suit any polynomial with arbitrary degree. In this paper, the zerofinding technique is considered for the class of unitary polynomials. Zerofinding unitary polynomials have been based to determine companion matrix eigenvalues. Let 𝑢(𝑧) be a unitary polynomial of degree 𝑛 as follows: 𝑢(𝑧)=𝑧^𝑛+𝑎_𝑛−1𝑧^𝑛−1+⋯+𝑎₀. (1.1) If 𝐴is its companion matrix associated with 𝑢, then det(𝐴−𝜆𝐼)=(−1)^𝑛𝑝(𝜆). (1.2) Conventional methods for numerically solving polynomials, and contemporary numerical methods from linear algebra, linear programming, and Fourier analysis, have been developed for the solution of (1.1). Most of these methods rely on a good initial approximation of the roots to ensure convergence besides stability considerations. It becomes the aim of this work to seek for an effective resolution that avoids the inaccuracy of root finding, in particular for the case of ill-conditioned algebraic or polynomial equations as in the case of higher degree polynomials and polynomials with closed or multiple roots.

The paper is organized as follows.

In Section 2, we have reviewed the iterative methods which have been used for finding roots of polynomials. In Section 3, the basis of the Fiedler’s theorems is reviewed. In Section 4, we have introduced Fiedler’s method by considering the initial values of Schmeisser’s method. In Sections 5 and 6, we have illustrated the solutions of polynomials by considering the initial values from a section of the complex plane and initial values from the circle with a certain radius, 𝑅. In Section 7, we have presented the results of choosing initial values for arbitrary degree polynomial in the Fiedler’s method to attain the convergence of the roots.

It is to be noted that in Sections 4, 5, 6, and 7 the tables given indicate the accuracy of our results. Moreover, the errors of the methods are shown by the figures. Importantly, in order to implement our methods and to obtain the results as illustrated by the figures and tables, we have utilized Matlab and Maple software. In Section 8, the analysis of the results is discussed. Finally, in Section 9, the conclusion of this research is given.

2. Review on Existing Methods

Graeffe’s root-squaring method replaces the given polynomial by another polynomial whose roots are the squares of the original polynomial. Newton’s method is an iterative procedure based on a Taylor series of the polynomial about the approximate root.

As for the study by Foster [3]: “Convergence requires a very good initial approximation of the exact root.” The algorithm of Jenkins and Traub involves three stages and the roots have to be computed in an approximately increasing order of magnitude in order to avoid instability that arises when deflating with a large root [4, 5]. The Laguerre’s method has cubic convergence for simple roots and also has linear convergence for multiple roots but each iteration requires that the first and second derivatives be evaluated at the estimated root, which makes the method computationally expensive [3, 6]. Trefethen and Toh [5, 7] studied on the convergence between roots of a given polynomial and eigenvalues of the Frobenius companion matrix [8] and also Traub and Reid have shown that these two sets are comparable.

For the case of polynomials with repeated roots, Hull and Mathon [9] presented an iterative polynomial zerofinding algorithm such that the iterations not only converge to simple roots but also converge to multiple roots. In 2005 Yan and Chieng [10] introduced a method that theoretically resolves the multiple-root issue. The proposed method adopts the Euclidean algorithm to obtain the greatest common divisor (GCD) of a polynomial and its first derivative. The multiple roots are then defaulted into simple ones and then the multiplicities of the roots are determined and calculated accordingly by applying conventional root-finding methods. In 2007, Winkler [11] denoted that GCD computations by Uspensky’s algorithm enable the multiplicity of each root to be calculated, and the initial estimates of the roots of a polynomial are obtained by solving several lower degree polynomials, all of whose roots are simple.

In some work, pejorative manifold have been applied. For example, Zeng [12] presented an algorithm which transforms the singular root-finding problem into a regular nonlinear least squares problem on a pejorative manifold and calculates multiple roots simultaneously from a given multiplicity structure and initial root approximations.

Besides stability considerations in most of the conventional zerofinding methods, convergence requires a good initial approximation of the exact roots. In this study, we consider the importance of choosing good initial approximation of the roots to ensure that convergence is attained. We present generally how to choose initial values by applying Fiedler’s theorems and remarks, and the hybrid between Schmeisser’s and Fiedler’s methods. The work partly focuses on the comparison of errors between the Schmeisser’s method and the Schmeisser-Fiedler’s method when the initial values for the Fiedler’s method are generated from the Schmeisser’s method, for solving the same polynomial. Moreover, this study also discusses the error of finding roots of a polynomial by using the Fiedler’s method, choosing initial values on a complex plane and on a circle. However, Malek and Vaillancourt [13] has similarly investigated on the finding of the roots of polynomials by choosing the initial values through the mentioned ways without paying attention to the comparison and condition of choosing desired initial values. In this study, we have especially investigated on the effects of attaining convergence, despite choosing only one initial value that is not sufficiently close to its exact value. The upcoming tables and figures show the associated error of the corresponding computations. What is more, the polynomials used in this research are not restricted to only a particular class of polynomials. It is also highlighted that one of the main tasks of this research is the implementation of all the methods that we have described here for solving polynomials and drawing related figures by Matlab and Maple software.

3. Fiedler’s Method

The basis of Fiedler’s method is a reflection of an important theorem in linear algebra: all roots of the characteristic polynomial of a real symmetric matrix are real. In fact, Fiedler’s method is Newbery’s expanded method [14] and it determines real symmetric matrix for polynomial with real roots. Required initial values in Fiedler’s method are chosen by some different ways: from the initial values supplied by Schmeisser’s method, randomly taken from a region in the complex plane, or from a circle with a large radius.

In the method of Fiedler, there are some important theorems for obtaining the companion matrix which are given as Theorems 3.1, 3.2, and 3.5 below. In fact, Fiedler’s Theorem is an advantage of general theorem described below.

Theorem 3.1 (see [15]). For𝑛≥1,assume that 𝑏₁,𝑏₂,…,𝑏_𝑛 are 𝑛 distinct numbers and 𝑣(𝑥)=^𝑛_𝑖=1𝑥−𝑏_𝑘. (3.1) Let 𝑢(𝑥) and 𝑤(𝑥) be polynomials of degree 𝑛 such that 𝑢(𝑏_𝑘)≠0 for each 𝑘=1,2,…,𝑛. Define 𝑛×𝑛 matrices A, C, and 𝐷 as follows.
For 𝑖,𝑘=1,2,…,𝑛, let 𝐴=diag(−𝑤(𝑏_𝑘)/𝑢(𝑏_𝑘)), and let 𝐶=(𝑐_𝑖𝑘) such that if𝑖≠𝑘𝑐_𝑖𝑘=𝑤𝑏_𝑖/𝑢𝑏_𝑖−𝑤𝑏_𝑘/𝑢𝑏_𝑘𝑏_𝑖−𝑏_𝑘,else𝑐_𝑖𝑘=𝑤(𝑡)𝑢(𝑡)^′_𝑡=𝑏. (3.2) Let𝑑=diag(𝑑_k) such that for a fixed constant δ≠ 0, 𝛿𝑣^′(𝑏_𝑘)𝑑²_𝑘−𝑢(𝑏_𝑘)=0is satisfied.
Then for each 𝑡 with  𝑣(𝑡)≠0, the number −𝑤(𝑡)/𝑉(𝑡)is an eigenvalue of (𝑢(𝑡)/𝑉(𝑡))𝐴−𝛿𝐷𝐶𝐷and 𝑑_𝑘/(𝑡−𝑏_𝑘)is the corresponding eigenvector.
We present the important result of the above theorem as follows.

Result 1. It is seen that by the selection of 𝑡₀as a root of 𝑢(𝑡) in the above theorem, the matrix 𝛿𝐷𝐶𝐷 will have an eigenvalue given by −𝑤(𝑡₀)/𝑉(𝑡₀) and this number will be equal to 𝑡₀ if unitary polynomials 𝑢(𝑥)and 𝑤(𝑥) are assumed such that 𝑤(𝑥)=𝑥[𝑢(𝑥)−𝑣(𝑥)] since we can write −𝑤(𝑡₀)/𝑉(𝑡₀)=−𝑡₀[𝑢(𝑡₀)−𝑉(𝑡₀)]/𝑉(𝑡₀)=𝑡₀.

Theorem 3.2 (see [16, Fiedler’s Theorem]). Assume that 𝑢(𝑥) is a unitary polynomial of degree   𝑛>1, and 𝑏₁,…,𝑏_𝑛∈𝐂 are 𝑛 distinct numbers such that 𝑢(𝑏_𝑘)≠0 for 𝑘=1,2…,𝑛. Consider 𝑣(𝑥)=^𝑛_𝑘=1𝑥−𝑏_𝑘,𝐵=diag𝑑_𝑘, (3.3) and define the matrix 𝐴 is a chaos of the matrix 𝐵, such that for a fixed 𝛿≠0, 𝐴=𝐵−𝛿𝑑𝑑^𝑇∈ℂ^𝑛×𝑛, (3.4) where, for 𝑖,𝑘=1,2,…,𝑛. If𝑖≠𝑘then𝑎_𝑖𝑘=−𝛿𝑑_𝑖𝑑_𝑘,  else𝑎_𝑘𝑘=𝑏_𝑘−𝛿𝑑²_𝑘such that 𝑑_𝑘 is a root of𝛿𝑣^′(𝑏_𝑘)𝑑²_𝑘−𝑢(𝑏_𝑘)=0, then det(𝐴−𝜆𝐼)=(−1)^𝑛𝑢(𝜆).
If roots of 𝑢 are distinct and real and 𝑏₁,…,𝑏_𝑛are approximations of the roots, then 𝛿 can be chosen as +1 or −1 in such a way that 𝑑_𝑘∈ℝ, thus 𝐴 is real symmetric [16].

Remark 3.3. If 𝑢(𝑥) and 𝑏_𝑘are all real then each 𝑑_𝑘is real or imaginary.

Remark 3.4. Schmeisser [16] have proved that if the roots of 𝑢 are distinct and single and the numbers 𝑏₁,…,𝑏_𝑛 are approximations of these roots, then the matrix 𝐴 in Theorem 3.2 is a unitary matrix.

Theorem 3.5 (see [14]). Let u(x) be a unitary polynomial of 𝑛≥1 such that 𝑢(𝑥)=𝑥^𝑛−𝑝𝑥^𝑛−1+𝑟(𝑥),deg𝑟≤𝑛−2 (3.5) and 𝑏₁,…,𝑏_𝑛 are complex and distinct numbers that 𝑢(𝑏_𝑘)≠0 for 𝑘=1,…,𝑛−1. Let 𝑉(𝑥)=^𝑛_𝑘=1𝑥−𝑏_𝑘,𝐵=diag𝑏_𝑘∈𝐶^{(𝑛−1)×(𝑛−1)}. (3.6) Assume that 𝐶=(𝑐_𝑘)∈𝐶^𝑛−1 is a column vector such that 𝑐_𝑘 satisfies 𝑉^′𝑏_𝑘𝑐²_𝑘+𝑈𝑏_𝑘=0𝑘=1,…,𝑛−1. (3.7) Then, there exists a bounded and symmetric matrix 𝐴∈𝐶^𝑛×𝑛, 𝑐𝐴=𝐵𝐶_𝑡𝑑𝑡ℎ𝑎𝑡𝑑=−𝑝−^𝑛−1_𝑘=1𝑏_𝑘det(𝐴−𝜆𝐼)=(−1)^𝑛𝑝(𝜆). (3.8)
If all the roots of 𝑢(𝑥)are simple and real and 𝑏_𝑘is approximation of these roots, then 𝐴 is real and symmetric, that is,𝐴^𝑇=𝐴∈𝑅^𝑛×𝑛. Thus the matrix 𝐴 is similar to Newbery’s matrix.
According to the aforementioned theorems and remarks, we can find the roots of polynomials with estimating initial values 𝑏₁,…,𝑏_𝑛 by using the methods of Fiedler and Schmeisser and also by generating the companion matrix𝐴, where𝐴=𝐵−𝛿𝑑𝑑^𝑡∈𝐶^𝑛×𝑛 and by the definition 𝑑_𝑘 as the root of𝛿𝑉^′(𝑏_𝑘)−𝑢(𝑏_𝑘)=0. We present some examples of solving polynomials by applying Fiedler’s Theorem and Schmeisser’s method. Further, we will examine the condition when only one of the approximations of the roots is far from its real value. For future study, we will go through another approach for estimation of the roots without much restriction and without compromising the convergence of the method to the exact solutions with a high degree of accuracy.

4. Hybrid of Fiedler’s Method and Schmeisser’s Method

Schmeisser [16] generated a symmetric tridiagonal matrix, 𝑇, by using a modified Euclidean algorithm. According to Schmeisser’s theorem which is based on a modified Euclidean algorithm and the matrix𝑇, we implemented the related algorithm using Matlab for solving monic polynomials. Consider a monic polynomial 𝑈(𝑥)and the corresponding matrix 𝑇,after solving, det(𝑇−𝜆𝐼)=(−1)^𝑛𝑈(𝜆), we obtain the roots of 𝑈(𝑥) approximately. In this method, we consider the obtained values of Schmeisser’s method as the desired initial values for Fiedler’s method.

Example 4.1. Consider the Wilkinson polynomial as follows: 𝑢(𝑥)=(𝑥−1)(𝑥−2)(𝑥−3)⋯(𝑥−19)(𝑥−20). (4.1) Using this method after ten iterations, we find the respective root of the polynomial and the results are shown in Table 1.

Table 1 

Achieved roots of Wilkinson polynomial by considering initial values of Schmeisser’s method.

Now, the error chart for the obtained results is given in Figure 1.

Figure 1 

Error of Schmeisser-Fiedler method Wilkinson polynomial.

The second column of Table 1 gives the eigenvalues of the matrix generated by Schmeisser’s method. These values correspond to the respective roots of the given polynomial, when applying Schmeisser’s method. Subsequently, the values which are obtained by Schmeisser’s method are used in Fiedler’s method as initial approximations of the roots and the eigenvalues of the associated companion matrix 𝐴 are then obtained. From row four to the last row of Table 1, it is clearly shown that the errors of solving the polynomial by Schmeisser’s method are higher than the errors accumulated from applying Fiedler’s method in which the desired initial values are acquired from Schmeisser’s method. Likewise, Figure 1 shows that the errors of Fiedler’s method by applying Schmeisser’s method for roots greater than 5 in Wilkinson polynomial decrease.

5. Fiedler’s Method Initial Values from a Section of the Complex Plane

In this method, we choose the initial values of Fiedler’s method taken from a section of the complex plane.

Example 5.1. Consider the polynomial𝑢(𝑥)=(𝑥−𝑖)(𝑥−2𝑖)(𝑥−3𝑖)(𝑥−4𝑖)(𝑥−5𝑖).
Using this method, we obtain the roots of this polynomial and the results are shown in Table 2.
The error chart is depicted in Figure 2.
The second column of Table 2 gives the approximated values of the roots when the initial values for solving polynomial of Example 5.1 are chosen from a complex plane using Fiedler’s method. Working on the generated matrix 𝑈, it was found that its eigenvalues converge to the respective real roots of the polynomial in the third column of Table 2. Referring to the fourth column, the errors of this method are adequately small in comparison with the real roots. Figure 2 depicts the results in Table 2, as well.

Table 2 

Achieved roots of polynomial in Example 5.1 by considering initial values of a complex plane.

Figure 2 

Error of Fiedler’s method for complex polynomial.

6. Fiedler’s Method with Initial Values from a Circle with 𝑅 Radius

In this method, we choose the initial values of Fiedler’s method from a circle with radius 𝑅. It should be taken care that the approximations converge to smaller roots if 𝑅is considered to be sufficiently large  (𝑅>10), and the method converges to larger roots if 𝑅 is assumed to be adequately small 𝑅<1).

Example 6.1. Consider the polynomial 𝑢(𝑥)=(𝑥−11)(𝑥−12)(𝑥−13)(𝑥−14) using this method, 𝑅=10 is chosen and we obtain the roots.
The results are shown in Table 3.
The error chart is depicted in Figure 3.
In Table 3, the second column points out the desired initial values for solving polynomial given in Example 6.1 by applying Fiedler’s method. They were taken from the circle with radius  𝑅=10. After computing the eigenvalues of matrix𝑈, given in the third column of Table 3, each corresponding to the respective roots of the polynomial, the errors of the method were satisfactorily small in comparison with real roots. Figure 3 illustrates the results in Table 3, as well.

Table 3 

Achieved roots of polynomial by choosing initial values on the circle.

Figure 3 

Error of R-Fiedler’s method for polynomial of degree 4.

7. Approximation of Initial Values for Fiedler’s Method for Arbitrary Degree Polynomial

In this part, after a set of research about the polynomial with each degree, we obtained that if we want to choose the initial values 𝑏1,𝑏2,…,𝑏𝑛 we are allowed to choose one of the roots to be away from the real roots but the others must be close to the real ones.

Example 7.1. Consider the below polynomial: 𝑢(𝑥)=𝑥⁴−10𝑥³+35𝑥²−50𝑥+24. (7.1) By considering the initial values as the second column in the table below, we obtain the roots of the polynomial after 10 iteration of Fiedler’s method. The results are listed in Table 4.
The error chart for the results obtained is given in Figure 4.

Table 4 

Achieved roots of polynomial in Example 7.1 by considering initial values (one of the approximations is away from real value).

Figure 4 

Error of Fiedler’s method for polynomial with arbitrary degree.

In Table 4, the second column shows the desired initial values for solving the polynomial, given in Example 7.1, by applying Fiedler’s method. In the second row, the amount of 13.3 is taken away from the exact value. In the third column, the eigenvalues of the matrix 𝑈 which corresponds to the respective roots of polynomial are shown. The results are appropriately close to the real roots. Figure 4 illustrates the results in Table 4, as well.

8. Discussion

Many numerical methods, using linear algebra, linear programming, and Fourier analysis, have been developed for the solution of the polynomial (1.1). In this stage, we describe the disadvantages of the present methods and explain the findings of our results in the form of tables and figures.

Considering the disadvantages of the zerofinding methods, Winkler mentioned that the Graeffe’s root-squaring method fails when there are roots of equal magnitudes [11, p. 3]; however, by applying Fiedler’s method the algebraic equations which have roots with almost the same modulus can be solved [17]. In addition, Bairstow’s method is only valid for polynomials containing real coefficients avoiding complex arithmetic. Moreover, the algorithm of Jenkins and Traub also involves three stages and is only valid for polynomials with real coefficients. Another insufficient method like Laguerre’s technique is not completely perfect whereby each iteration requires that the first and second derivatives be evaluated at the estimated root, which makes the method computationally expensive. Muller’s method is a variant of Newton’s method and convergence in Newton’s method requires that the estimate be sufficiently near the exact root.

It can be gathered that the above methods have been facing some issues which need to be reviewed. The information in Table 1 shows that after choosing the desired initial values from the results obtained by Schmeisser’s method, the third column of Table 1, the approximate results are reasonable, having the accuracy of nearly 10⁻¹⁰ after ten iterations. By comparing the results in columns 4 and 5 of Table 1, it reveals that the Fiedler’s method, assuming the desired initial values taken from obtained values of Schmeisser’s method, is more accurate than solving the polynomial by Schmeisser’s method entirely.

It can be seen that 75 percent of the roots have accuracy up to almost 10⁻¹⁶. Similarly, Figure 1 also verifies that in case of roots which are greater than 5 the error of Fiedler’s method in which Schmeisser’s method is applied steadily decreases.

The information in Table 2 points out the estimated initial values which are chosen of a complex plane. The results obtained by Fiedler’s method in Example 4.1 are reasonable and nearly have accuracy of 10⁻¹¹. Figure 2 confirms the same results as well.

Choosing suitable initial values on the circle with 𝑅=10 in Example 6.1 along with comparison of the third column in Table 3 and the real roots of the polynomial concludes that the results which were found by using this method are reasonable. These results roughly have accuracy of 10⁻¹⁴. Likewise, Figure 3 confirms the similar findings.

In the second column of Table 4 while the real roots are {1,2,3,4}, only one of approximation of the roots is chosen away from the exact value. In Example 7.1, we have considered an initial value approximately equals 13.3 for the real root 2. According to the third column of Table 4, the eigenvalues of matrix U correspond to the roots of polynomial. the results were adequately close to the real roots with an accuracy of 10⁻¹⁰. Similarly, Figure 4 also proves this statement.

9. Conclusion

Fiedler’s different algorithms are described. As mentioned earlier, it can be seen that among existing numerical algorithms, we are not able to say that there is a special algorithm for every arbitrary polynomial that is better than other ones and also there are the zerofinding explicit formulas for maximum fifth-degree polynomial. In order to find the roots of an arbitrary polynomial, we could find the roots of polynomial with high accuracy by using one of the algorithms presented in this paper. In the case of using these algorithms for choosing the initial values, we are able to choose these values from Schmeisser’s method or by selection from a square or circle or by an arbitrary selection that all values must be closed to the real ones except for one of them. In addition, besides stability considerations, in future work we are interested to find the root-finding algorithms with less limitation of good initial approximation of the roots to ensure convergence besides stability considerations. In this case, future studies should consider whether we can find an approach of polynomial zerofinding which ensures convergence to the roots even though some of the initial values may not necessarily be closed to the real roots.

Acknowledgments

The authors would like to acknowledge UTM Research University Grant, vote no. Q.J130000.7126.04J05, Ministry of Higher Education (MOHE), Malaysia, for supporting the research. The authors are thankful to the referees for their constructive comments which improved the presentation of the paper.