This module provides a QR based polynomial solver.
at the start of your source file.
and a QR based polynomial solver.
The remainder of the page documents first the functions for evaluating, computing
polynomials, computing estimates about polynomials and next the QR based polynomial
solver.
@section polynomialUtils convenient functions to deal with polynomials
@subsection roots_to_monicPolynomial
The function
@code
void roots_to_monicPolynomial( const RootVector& rv, Polynomial& poly )
\endcode
computes the coefficients \form#34 of
where \form#36 is known through its roots i.e. \form#37.
@subsection poly_eval
The function
@code
T poly_eval( const Polynomials& poly, const T& x )
\endcode
evaluates a polynomial at a given point using stabilized Hörner method.
The following code: first computes the coefficients in the monomial basis of the monic polynomial that has the provided roots;
then, it evaluates the computed polynomial, using a stabilized Hörner method.
\include PolynomialUtils1.cpp
Output:
Roots: 0.680375 -0.211234 0.566198 0.59688
Polynomial: -0.04857.x^0+ 0.00860842.x^1+ 0.739882.x^2+ -1.63222.x^3+ 1.x^4
Evaluation of the polynomial at the roots: -2.08167e-17 0 0 2.08167e-17
@subsection Cauchy bounds
The function
@code
Real cauchy_max_bound( const Polynomial& poly )
\endcode
provides a maximum bound (the Cauchy one: \form#38) for the absolute value of a root of the given polynomial i.e.
root of
,
The leading coefficient
: should be non zero
.The function
provides a minimum bound (the Cauchy one:
) for the absolute value of a non zero root of the given polynomial i.e.
root of
,
@section QR polynomial solver class
Computes the complex roots of a polynomial by computing the eigenvalues of the associated companion matrix with the QR algorithm.
The roots of \form#46 are the eigenvalues of
However, the QR algorithm is not guaranteed to converge when there are several eigenvalues with same modulus.
Therefore the current polynomial solver is guaranteed to provide a correct result only when the complex roots \form#48 have distinct moduli i.e.
![$ \forall i,j \in [1;d],~ \| r_i \| \neq \| r_j \| $](form_49.png)
.
With 32bit (float) floating types this problem shows up frequently.
However, almost always, correct accuracy is reached even in these cases for 64bit (double) floating types and small polynomial degree (<20).
\include PolynomialSolver1.cpp
In the above example:
-# a simple use of the polynomial solver is shown;
-# the accuracy problem with the QR algorithm is presented: a polynomial with almost conjugate roots is provided to the solver.
Those roots have almost same module therefore the QR algorithm failed to converge: the accuracy
of the last root is bad;
-# a simple way to circumvent the problem is shown: use doubles instead of floats.
Output:
Roots: 0.680375 -0.211234 0.566198 0.59688 0.823295
Complex roots: (-0.211234,0) (0.566198,0) (0.59688,0) (0.680375,0) (0.823295,0)
Real roots: -0.211234 0.566198 0.59688 0.680375 0.823295
Illustration of the convergence problem with the QR algorithm:
---------------------------------------------------------------
Hard case polynomial defined by floats: -0.957 0.9219 0.3516 0.9453 -0.4023 -0.5508 -0.03125
Complex roots: (1.19707,0) (0.705138,0) (-1.9834,0) (-0.396564,0.966801) (-0.396564,-0.966801) (-16.7513,0)
Norms of the evaluations of the polynomial at the roots: 1.89088e-06 2.32458e-06 1.13394e-05 1.29018e-06 1.29018e-06 0.164618
Using double's almost always solves the problem for small degrees:
-------------------------------------------------------------------
Complex roots: (1.19707,0) (0.70514,0) (-1.9834,0) (-0.396564,0.966801) (-0.396564,-0.966801) (-16.7513,0)
Norms of the evaluations of the polynomial at the roots: 3.78175e-07 0 2.0411e-06 2.48518e-07 2.48518e-07 0
The last root in float then in double: (-16.75127792,0) (-16.75128099,0)
Norm of the difference: 3.814697266e-06