Numerical methods for finding the roots of polynomials with real and complex coefficients

Sklyar Alexander Yakovlevich PhD in Technical Science Associate Professor; Department of Applied Mathematics; Russian Technological University (MIREA) 119602, Russia, Moscow, Vernadsky , 78 askliar@mail.ru Other publications by this author

Introduction

There are many tasks of a very different nature in which the definition of the roots of polynomials is required. Algebraic equations arise in the study of equilibrium states of complex thermodynamic and mechanical systems, they often appear in aerodynamics, in flight mechanics. For example, the speed of the fastest climb of an airplane is determined from an algebraic equation of the eighth degree. Algebraic equations also arise when performing various geometric calculations – determining the points of intersection and conjugation of curved contours, when designing smooth surfaces, well-streamlined bodies, and many other tasks.

Methods for solving equations to the second degree have been known since the time of ancient Greece. In the XVI century, analytical expressions for polynomials of degree 3 (Cordano formula) and degree 4 were obtained by Scipio del Ferro, Tartaglia, Ferrari ^[3].

In the early 20s, Abel and then Galois in the 30s of the XIX century proved that such formulas for equations of the nth degree in the general case for any n ≥ 5 obviously cannot be found.

Various methods are currently used for the numerical approximate solution of equations of higher degrees, such as the Lobachevsky method, the Hitchcock method, the Gorner scheme for dividing a polynomial into a binomial and a square trinomial ^{[1, 4]} and others.

Another type of algorithms based on obtaining iterative formulas, by extracting simple and quadratic multipliers from polynomials, followed by comparing the entries of polynomials with the remainder when the roots are approximate, and without remainder when the values of the roots are exact, is considered in articles by Chye Yong Un and A.B. Shein ^[5,6,7].

A number of algorithms for solving polynomial problems and an extensive bibliography are given by G. P. Kutischev ^[8]. It is also worth noting the approaches for solving them in ^{[9, 10, 11]}.

Note that the very existence of a large number of different methods generally indicates that there is not a single "completely satisfactory" one.

1. Removing multiple roots

A polynomial of degree n has exactly n roots and is represented as

In general, not all roots of x _k are different.

Finding the roots of a polynomial encounters certain difficulties in cases where it has multiple roots.

Let there be a multiple root x* with multiplicity m, then P _n(x) can be represented as

Therefore, P _n(x) and P_'n(x) will have a common divisor (x-x*)^m-1.

Then, using Euclid's algorithm, we can find the largest divisor Q(x) of the polynomials P _n(x) and P_'n(x).

And the original polynomial P _n(x) can be represented as

The polynomial H(x) will have the same roots as P _n(x), but will not have multiple roots. Thus, the initial task is reduced to finding the roots of a polynomial whose roots are all different.

2. Determining the range of roots

Consider the equation (a ₀ will be assumed to be equal to 1)

According to the well-known theorem ^[1,2], all roots x _k (k = 1, 2, ..., n) of a polynomial in the complex plane lie in a ring

r <| x_k |< R, where

At the same time, we should immediately note that this statement can be somewhat strengthened.

Let's replace the variable x=y/t. In this case, the original equation with respect to y will take the form

Given the nature of the substitution, you can write

By changing the value of t, it is possible to reduce the estimate of the upper bound.

Let the maximum value take the modulus of the coefficient at y _m. Then the score for this coefficient will be no less than

t* is located from

Thus, it can be argued that the estimate lies between the values corresponding to the current value of t and t*.

At the same time, when t changes, other coefficients also change. Thus, we get restrictions on the intersection of the lines of their change.

If t _r does not lie within the interval t and t*, then it can be ignored. From the rest, select the value closest to t.

If there are none, then the optimal t is _r=t*.

If the signs of derivatives

if they match, then further improvement is impossible. Optimal t _r=t*.

If the signs of the derivatives match, then we take m equal to k and repeat the procedure.

Consider the derivative of the k coefficient at the point t

Consider as an initial polynomial of the fourth degree

t⁴-10t³+35t²-50^t+24=(t-1)(t-2)(t-3)(t-4)

The initial upper bound is 51, m=1, t=1,

Points of intersection

And finally we get

3. Search for valid roots

Let's consider separately the search for positive and negative roots.

Let the initial polynomial have the form (the highest coefficient a ₀, without reducing generality, can be assumed to be equal to 1).

Let's introduce b _k=(a k_+|ak|)/2 and c _k=(a k_-|ak|)/2. For a _k>0 b _k=a _k, c _k=0 for a _k≤0 b _k=0. c _k=-a _k. In these designations

The polynomials P₊(x) and P_-(x) for x>0 are continuous non-negative monotonically increasing functions.

Using the introduced polynomials P₊(x) and P_-(x) and the upper bound of the values of the roots of R, we will look for roots in the range [A, B], where A=0, B=R. The lower bound can be refined to the value r* by obtaining the range [r*, R*] instead of the range [r, R], but this does not matter in principle.

Let's take a closer look at the algorithm for finding roots.

Calculate the values of P₊(x) and P_-(x) at the ends of the interval and call the Root function to calculate the roots with the parameters A, B, P₊(A), P_-(A) and P₊(B), P_-(B). The Root function returns either the found root (value greater than 0), or -1 if the polynomial P(x) has no roots in this interval.

The algorithm of the Root function has the following form.

1. Calculate the values of the initial polynomial at the ends of the interval P(A)=P₊(A)-P_-(A) and P(B)= P₊(B)-P-(B)

2. If P(A)=0, we return the found root A, if P(B)=0, we return the found root B, if P(A)P(B), we call the standard STROOT procedure for finding roots in a given interval [A, B] (for example, by dichotomy or chords) and return, the value found for it.

3. Calculate Q ₁=P₊(A)-P_-(B) – the lower bound of the values of P(x) in the interval [A, B] and Q ₂=P₊(B)-P_-(A) – the upper bound of the values of P(x) in the interval [A, B].

4. If Q ₁ Q ₂≥0, then there are no roots in this interval and return the value -1.

5. Calculate C=(A+B)/2.

Calling the Root function to calculate the roots with the parameters A, C, P₊(A), P_-(A) and P₊(C), P_-(C).

6. If the calculated value of x is >0, then we return the found root of x.

7. Calling the Root function to calculate the roots with the parameters C, B, P₊(C), P_-(C) and P₊(B), P_-(B). We return the found root x (if there are no roots, then -1).

The above recursive algorithm allows you to find the root of a polynomial in the specified range or make sure that there are no valid roots inside it.

After finding the root u, we proceed to search for the next root, replacing the original polynomial P _n(x) with a polynomial of lesser degree P _n-1(x), P _n(x)=(x-u)P _n-1(x). The polynomial P _n-1(x) has all remaining roots P _n(x). Repeating the above algorithm, we find all the positive roots.

To find negative roots, replace the original polynomial P _n(x) with the polynomial (-1)ⁿ P _n(-x), the roots of which have the same values, but with the opposite sign. Having found all the positive roots of the second polynomial, we will thereby find all the real roots of the original polynomial.

Since there are no other real roots, the remaining roots represent pairs of complex conjugate numbers and the remaining polynomial has an even degree.

4. Search for complex roots

Let the polynomial P _{2 n}(x) with real coefficients have no real roots, then it is represented as

Consider the representation of a polynomial of degree 2 n in the form

To standardize calculations for all m, we introduce q_-2=q_-1=q _n-1=q _n≡0, q _n-2=1.

We will find the roots based on optimizing the approximation of the original polynomial Pn(x) by the polynomial Rn(x).

Let's introduce the function

For given values p ₀, p ₁, the coefficients q _i are found from the minimization requirement F(p,q).

Or

For given values of q _i, the coefficients p ₀, p ₁ are found from the minimization requirement F(p,q).

Or

To find the polynomial x ²+p ₁ x+p ₀, you can use the following iterative algorithm.

1. Set the initial value of the coefficients p ₀, p ₁.

2. Based on (4.5), we calculate the coefficients q _i. The problem is reduced to solving a system of linear equations (SLA) for a 5-diagonal matrix. Given the type of matrix, the problem is reduced to a number of preparatory operations of difficulty O(n ²) and solving a SLOUGH of dimension 3×3.

3. The resulting set of coefficients q _i is used to obtain adjusted values p _0, p 1 in accordance with (4.6).

4. We estimate the magnitude of the calculation error in accordance with (4.3). If the error is less than the specified threshold, then we finish the algorithm. Otherwise, we proceed to step 2 of the algorithm.

As a result, we get the polynomial x ₂ +p ₁ x+ p ₀ giving the first pair of roots. The coefficients of qi give the polynomial

which can be used by the above algorithm to obtain the next pair of roots.

Note that the convergence rate of this algorithm is low. To speed up, you can use the processing of previously obtained values.

Let's assume that at the ith and i+1st steps we got the values p _i,1, p _i+1,1 and p _i,0, p _i+1,0 with residuals F(p,q) – Fa, Fb, respectively.

Calculate the values p _{c 1}=p _i,1+2(p _i+1,1-p _i+1,1) and p _{c 0}= p _i+1,0+2(p _i+1,0-p _i,0) and the corresponding Fc value.

If F _c<F _b, replace the points and the values in them a=b, b=c and repeat the calculation. If F _c<F _b, then by quadratic interpolation we find the point c*, taking it as a suboptimal value and return to the main algorithm.

At the optimum point, both (4.5 and 4.6) must be performed simultaneously

In matrix form, the aggregate (4.5, 4.6) is represented as

At |A|≠0, the algorithm converges to the solution of the problem.

In addition, from the point of view of practical implementation, the hybrid algorithm looks much more efficient.

Several iterations of the basic algorithm lead us to the vicinity of the root. Considering that in the absence of multiple roots |P_'n(x)| is guaranteed not to turn to 0, for P _n(x)=0, a scheme using derivative values (Newton's method) or its more subtle modifications ^[12]) provides sufficiently fast convergence.

Let the approximation s=x ²+p ₁ x+p ₀ be obtained. Its root will be

The value of P _n(z) in the neighborhood of z ₀ is represented as

This allows you to calculate the approximation of the z* root.

At z ₀, close to z*, the algorithm converges quite quickly. As a result, we get for s=x ²+p ₁ x+p ₀

The resulting expression is used to obtain the polynomial Q _n-2(x) from P _n(x)=(x ²+p ₁ x+p ₀)Q _n-2(x), after which we look for the following roots, polynomials that coincide with the roots of Q _n-2(x).

5. Search for roots of polynomials with complex coefficients

Consider the equation (a ₀ will be assumed to be equal to 1). The coefficients a _k are complex numbers.

Let its roots be x ₁, x ₂, ..., x _n, then the roots of the equation

will be

Let's construct a polynomial R _{2 n}(x)=P _n(x)Q _n(x). This polynomial has roots

The quantities b _k are representable as sums and products of complex conjugate numbers and, therefore, are real numbers. Thus, the polynomial R _{2 n}(x) is a polynomial with real coefficients.

Accordingly, the search for the roots of the polynomial P _n(x) with complex coefficients is reduced to the search for the roots of the polynomial R _{2 n}(x) with real coefficients. Note that all its real roots will be multiples, and some of the obtained complex conjugate roots will be extraneous.

Conclusion

The sequence of the proposed algorithms makes it possible to find all the real and complex roots of the polynomial.

To find the roots of a polynomial of degree n with real coefficients, an algorithm is constructed that includes the following main steps:

· definition of multiple roots,

· highlighting a range of roots;

· finding intervals guaranteed to contain roots (in time O(ln(R/H)), where R is an estimate of the total range of roots, and H is the length of the interval guaranteed to contain the root);

· iterative construction of trinomials that serve as an estimate of the values of pairs of complex conjugate roots (the number of iterations does not necessarily have to meet the requirements for the accuracy of the solution, it is enough to fulfill the requirement of root isolation, after which the solution is achieved by traditional methods, for example, the Newton method).

To find the roots of a polynomial of degree n with complex coefficients, an auxiliary polynomial of degree 2 n with real coefficients is constructed, containing all the roots of the original polynomial of degree n with complex coefficients, for which the above algorithm is used. The computational complexity of the algorithm is quite acceptable for calculating the roots of polynomials of almost any degree found in real problems.

The software implementation of the proposed algorithms is performed in C++ in the Windows operating system environment. Given the computational nature of the algorithms, transferring them to a different operating environment does not cause any problems.