![]() |
Eigen-unsupported
3.4.0
|
This module provides a QR based polynomial solver.
To use this module, add
at the start of your source file.
Classes | |
class | Eigen::PolynomialSolver< _Scalar, _Deg > |
A polynomial solver. More... | |
class | Eigen::PolynomialSolverBase< _Scalar, _Deg > |
Defined to be inherited by polynomial solvers: it provides convenient methods such as. More... | |
Functions | |
template<typename Polynomial > | |
NumTraits< typename Polynomial::Scalar >::Real | Eigen::cauchy_max_bound (const Polynomial &poly) |
template<typename Polynomial > | |
NumTraits< typename Polynomial::Scalar >::Real | Eigen::cauchy_min_bound (const Polynomial &poly) |
template<typename Polynomials , typename T > | |
T | Eigen::poly_eval (const Polynomials &poly, const T &x) |
template<typename Polynomials , typename T > | |
T | Eigen::poly_eval_horner (const Polynomials &poly, const T &x) |
template<typename RootVector , typename Polynomial > | |
void | Eigen::roots_to_monicPolynomial (const RootVector &rv, Polynomial &poly) |
The remainder of the page documents first the functions for evaluating, computing polynomials, computing estimates about polynomials and next the QR based polynomial solver.
The function
computes the coefficients \( a_i \) of
\( p(x) = a_0 + a_{1}x + ... + a_{n-1}x^{n-1} + x^n \)
where \form#44 is known through its roots i.e. \form#45.
The function
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.
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
The function
provides a maximum bound (the Cauchy one: \(C(p)\)) for the absolute value of a root of the given polynomial i.e. \( \forall r_i \) root of \( p(x) = \sum_{k=0}^d a_k x^k \), \( |r_i| \le C(p) = \sum_{k=0}^{d} \left | \frac{a_k}{a_d} \right | \) The leading coefficient \( p \): should be non zero \(a_d \neq 0\).
The function
provides a minimum bound (the Cauchy one: \(c(p)\)) for the absolute value of a non zero root of the given polynomial i.e. \( \forall r_i \neq 0 \) root of \( p(x) = \sum_{k=0}^d a_k x^k \), \( |r_i| \ge c(p) = \left( \sum_{k=0}^{d} \left | \frac{a_k}{a_0} \right | \right)^{-1} \)
Computes the complex roots of a polynomial by computing the eigenvalues of the associated companion matrix with the QR algorithm.
The roots of \( p(x) = a_0 + a_1 x + a_2 x^2 + a_{3} x^3 + x^4 \) are the eigenvalues of \( \left [ \begin{array}{cccc} 0 & 0 & 0 & a_0 \\ 1 & 0 & 0 & a_1 \\ 0 & 1 & 0 & a_2 \\ 0 & 0 & 1 & a_3 \end{array} \right ] \)
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#56 have distinct moduli i.e.
\( \forall i,j \in [1;d],~ \| r_i \| \neq \| r_j \| \).
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.70514,0) (-1.9834,0) (-0.396563,0.966801) (-0.396563,-0.966801) (-16.7513,0) Norms of the evaluations of the polynomial at the roots: 3.08019e-06 2.98023e-07 2.10915e-05 5.35758e-07 5.35758e-07 0 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.75128174,0) (-16.75128099,0) Norm of the difference: 0
|
inline |
[in] | poly | : the vector of coefficients of the polynomial ordered by degrees i.e. poly[i] is the coefficient of degree i of the polynomial e.g. \( 1 + 3x^2 \) is stored as a vector \( [ 1, 0, 3 ] \). |
|
inline |
[in] | poly | : the vector of coefficients of the polynomial ordered by degrees i.e. poly[i] is the coefficient of degree i of the polynomial e.g. \( 1 + 3x^2 \) is stored as a vector \( [ 1, 0, 3 ] \). |
|
inline |
[in] | poly | : the vector of coefficients of the polynomial ordered by degrees i.e. poly[i] is the coefficient of degree i of the polynomial e.g. \( 1 + 3x^2 \) is stored as a vector \( [ 1, 0, 3 ] \). |
[in] | x | : the value to evaluate the polynomial at. |
|
inline |
[in] | poly | : the vector of coefficients of the polynomial ordered by degrees i.e. poly[i] is the coefficient of degree i of the polynomial e.g. \( 1 + 3x^2 \) is stored as a vector \( [ 1, 0, 3 ] \). |
[in] | x | : the value to evaluate the polynomial at. |
void Eigen::roots_to_monicPolynomial | ( | const RootVector & | rv, |
Polynomial & | poly | ||
) |
Given the roots of a polynomial compute the coefficients in the monomial basis of the monic polynomial with same roots and minimal degree. If RootVector is a vector of complexes, Polynomial should also be a vector of complexes.
[in] | rv | : a vector containing the roots of a polynomial. |
[out] | poly | : the vector of coefficients of the polynomial ordered by degrees i.e. poly[i] is the coefficient of degree i of the polynomial e.g. \( 3 + x^2 \) is stored as a vector \( [ 3, 0, 1 ] \). |