Zeros of Polynomial Equations in the Complex Plane

Abstract

In this article we will discuss Laguerre's algorithm for finding zeros of polynomial equations of degree n in the complex plane. We will derive a particular case and analyze convergence for simple roots. We will discuss and verify some examples using the implemented algorithm, showing some of its virtues and limitations.

1. Introduction

The problem of finding zeros of a polynomial of degree n has substantially influenced the development of mathematics through the centuries and is still quite important in practical and theoretical applications.

Currently, the study of polynomials of the form

$$p(z) = \sum_{i=0}^n a_i z^i \quad (1)$$

does not play a central role in mathematics or applied mathematics. In particular, many of the problems that arise are linearized and solved using linear algebra tools and fast Fourier transforms, which may involve roots of polynomial equations.

However, equations of this type are still of interest in some areas of applied mathematics and algorithms continue to be developed and studied.

As was proven in 1827 by Abel, for polynomials of degree equal to or greater than 4 it is not possible to find a formula to obtain the exact solution, and therefore, any problem that requires the solution of such a polynomial requires us to use numerical methods.

In the present article we discuss Laguerre's method in the complex plane and an implementation of it with some different types of polynomials.

2. Theorical Results

2.1 Initial Definitions and Theorems

Definition 1: We say that an algorithm has convergence order

$$\alpha$$

at root

$$\eta$$

if,

$$\lim_{n \rightarrow \infty} \frac{|z_{n+1} - \eta|}{|z_n - \eta|^\alpha} = C \quad (2)$$

Definition 2: A numerical method is said to be global if it does not require an initial value sufficiently close to the zero of the polynomial to converge, otherwise we say it is local.

Theorem 1: Let p(z) be a polynomial of degree n, n > 1, with complex coefficients. Then p(z) has n complex roots.

Theorem 2: Let p(z) be a normalized polynomial with degree m > 4 and let

$$P = {p_i : i \in {1,...,l}}$$

be the set of roots. For

$$u_0 \in \mathbb{C} \setminus P$$

with

$$p'(u_0) \neq 0$$

or

$$p''(u_0) \neq 0$$

the sequence

$$L_p(u_0) = (u_n)_{n \in \mathbb{N}}$$

defined by:

$$u_{n + 1} = u_{n} - \frac{m}{q(u_{n}) + s(u_{n})r(u_{n})} \quad (3)$$

with

$$n \in \mathbb{N}$$

and

$$q(z) = \frac{p'(z)}{p(z)}$$

1.

$$r(z) = \sqrt{(m-1)(mt(z) - q^2(z))}$$

where

$$t(z) = q^2(z) - \frac{p''(z)}{p(z)}$$

and

$$s(z) = \begin{cases} 1, & \text{if } (\text{Re}\, q(z))(\text{Re}\, r(z)) + (\text{Im}\, q(z))(\text{Im}\, r(z)) \geq 0 \ -1, \end{cases}$$

otherwise, with

$$q(u_n) + s(u_n)r(u_n) \neq 0$$

If there exists a simple root

$$p \in P$$

such that

1. $$|u_0 - p| \leq \frac{1}{2m-1} \min{|p - p_i|, p \in P}$$

then

$$L_p(u_0)$$

converges to the limit p and

1. $$|u_n - p| \leq \lambda^n|u_0 - p|$$

with

$$\lambda = \frac{15} {16}$$

for all

$$n \in \mathbb{N} \setminus {0}$$

If all roots of p(z) are simple, then in 2.,

$$\mu_p = \min{|p - p_i|, p \in P \setminus {p_i}}$$

can be replaced by the lower bound

$$(1 + \frac{1}{d_0} \max_{1 \leq j \leq \binom{m}{2}} |d_j|)^{-\frac{1}{2}} < \min{\mu_p | p \in P}$$

with coefficients

$$d_j$$

of the polynomial

$$D_p(z) = \prod_{1 \leq i < k \leq m}(z - (p_i - p_k)^2) = \sum_{j=0}^{\binom{m}{2}} d_j z^j \quad (4)$$

These coefficients can be calculated using only coefficients of p(z). All roots of p(z) are simple if, and only if,

$$d_0 \neq 0$$

Theorem 3: Let

$$P(z) = \sum_{i=0}^{n} a_i z^i$$

be a polynomial. We can evaluate its value at a point using the following procedure:

$$\begin{cases} P_0(z) = a_n \ P_i(z) = a_{n-i} + zP_{i-1}(z), \quad i \in {1,...,n} \end{cases}$$

The proof is immediate, simply rewriting the polynomial.

Theorem 4 (Horner's Algorithm for Polynomial Division): Let

$$P_n(z) = (z-r)P_{n-1}(z)$$

be a polynomial.

Then, if

$$P_{n-1}(z) = b_0 + b_1z + b_2z^2 + ... + b_{n-1}z^{n-1}$$

, we have:

$$b_{n-1} = a_n, \quad b_{n-2} = a_{n-1} + rb_{n-1}, \quad b_0 = a_1 + rb_1$$

To prove this result, simply equate the coefficients on both sides.

2.2 Laguerre's Method

Let

$$P_n(z)$$

be a polynomial of degree

$$n \geq 1$$

By Theorem 1, we have

$$P_n(z) = (z - z_1)(z - z_2)...(z - z_n) \quad (5)$$

therefore,

$$\log|P_n| = \sum_1^n \log|z - z_i| \quad (6)$$

and then,

$$\frac{d\log|P_n|}{dz} = \sum_1^n \frac{1}{(z - z_i)} = \frac{P_n'}{P_n} \quad (7)$$

$$-\frac{d^2\log|P_n|}{dz^2} = \sum_1^n \frac{1}{(z - z_i)^2} = -\frac{P_n''}{P_n} \quad (8)$$

We define equations (7) and (8) as G and H, respectively, we have

$$H = G^2 - R$$

where

$$R = \frac{P_n''}{P_n}$$

.

To proceed, we will assume that the sought root

$$z_1$$

is isolated from the others by a distance, such that

$$a = z - z_1 \quad (9)$$

$$b = z - z_i, \quad i \in {2,...,n} \quad (10)$$

Thus,

$$\frac{1}{a} + \frac{n-1}{b} = G \quad (11)$$

$$\frac{1}{a^2} + \frac{n-1}{b^2} = H \quad (12)$$

And therefore,

$$a = \frac{n}{G \pm \sqrt{(n-1)(nH - G^2)}} \quad (13)$$

Such that we take

$$\max{\pm\sqrt{(n-1)(nH - G^2)}}$$

so that the step has the smallest possible size.

Since the factor inside the root can be negative, a can be a complex number.

If the root of the polynomial is simple, it is proven that the method converges cubically. If the root is multiple, we have linear convergence.

Theorem 2 determines for each simple root a disk centered on it, such that for all initial values in this neighborhood the convergence of the root is guaranteed.

Thus, this method offers the possibility of convergence to a complex root and if the polynomial has only real roots, it is guaranteed that the method converges independently of the given initial point, therefore convergence in this case is global.

3. Algorithm

The algorithm was developed in C language and designed to run on compilers like GCC that support the C99 and C11 standards. In these standards, the complex.h library is implemented with functions for working with complex numbers in a simple and elegant way.

Using some of the results from the section above, the necessary libraries were built for calculating the value of the polynomial at a point, polynomial division, and the method's algorithm.

The evalpoly.h and deflation.h libraries were developed from Theorems 3 and 4, respectively, just as the main algorithm, laguerre.h, was made according to what was developed in the second part of the previous section.

Rounding errors occur mainly when the deflate (deflation.h) and evalpoly (evalpoly.h) functions are called. Therefore, the calculations were done with long double to try to minimize rounding errors.

4. Results and Discussion

Several tests were performed, among which those listed in the table below.

For low-degree polynomials, we note that the precision of the solution is remarkable, allowing, for example, to calculate the square root approximation with up to 13 decimal places of precision.

We can note, according to the results obtained, that for roots with multiplicity less than or equal to 10, the algorithm behaves exceptionally well, since it was not developed with this type of problem in mind. From higher degrees onwards, we begin to notice problems such as divergence of the polynomial solution, which could not be circumvented by increasing the maximum number of iterations (since, as this method converges linearly for multiple roots, we expect to need a large number of iterations for this type of problem).

For polynomials with simple roots and slightly higher degrees

$$n > 10$$

we notice better behavior, but as shown in the penultimate and last polynomials in the table, the algorithm diverged in calculating the last root, finding only 13 of 14 existing roots.

We also note that as the degree of the polynomial increases considerably, the precision drops drastically.

Test Results Table

5. Conclusion

For a first implementation, the algorithm presents good performance, but with more time and in a longer project it would be possible to considerably improve its performance, since in its current state it could not be used in a production environment to solve high-performance problems (such as in software like Mathematica or MATLAB).
The complete project can be found here

References

1. Pan Y., Victor, Solving a Polynomial Equation: Some History and Recent Progress, SIAM Review, Vol. 39, No. 2 (Jun., 1997), pp. 187-220.

2. Ahlfors, Lars, Complex Analysis, McGraw-Hill, 1979.

3. Moller, Herbert, Convergence and visualization of Laguerre's rootfinding algorithm, January 9, 2015.

4. Acton S. Forman, Numerical Methods that work, The Mathematical Association of America, 1990.

5. Kiusalaas, Jaan, Numerical Methods in Engineering with Python 3, Cambridge, 2013.