$treeview $search $mathjax
Eigen-unsupported
3.2.5
$projectbrief
|
$projectbrief
|
$searchbox |
This module provides a QR based polynomial solver. More...
Classes | |
class | PolynomialSolver< _Scalar, _Deg > |
A polynomial solver. More... | |
class | PolynomialSolverBase< _Scalar, _Deg > |
Defined to be inherited by polynomial solvers: it provides convenient methods such as
| |
Functions | |
template<typename Polynomial > | |
NumTraits< typename Polynomial::Scalar >::Real | cauchy_max_bound (const Polynomial &poly) |
template<typename Polynomial > | |
NumTraits< typename Polynomial::Scalar >::Real | cauchy_min_bound (const Polynomial &poly) |
template<typename Polynomials , typename T > | |
T | poly_eval (const Polynomials &poly, const T &x) |
template<typename Polynomials , typename T > | |
T | poly_eval_horner (const Polynomials &poly, const T &x) |
template<typename RootVector , typename Polynomial > | |
void | roots_to_monicPolynomial (const RootVector &rv, Polynomial &poly) |
This module provides a QR based polynomial solver.
To use this module, add
#include <unsupported/Eigen/Polynomials>
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.
The function
void roots_to_monicPolynomial( const RootVector& rv, Polynomial& poly )
computes the coefficients of
where is known through its roots i.e.
.
The function
T poly_eval( const Polynomials& poly, const T& x )
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 <unsupported/Eigen/Polynomials> #include <iostream> using namespace Eigen; using namespace std; int main() { Vector4d roots = Vector4d::Random(); cout << "Roots: " << roots.transpose() << endl; Eigen::Matrix<double,5,1> polynomial; roots_to_monicPolynomial( roots, polynomial ); cout << "Polynomial: "; for( int i=0; i<4; ++i ){ cout << polynomial[i] << ".x^" << i << "+ "; } cout << polynomial[4] << ".x^4" << endl; Vector4d evaluation; for( int i=0; i<4; ++i ){ evaluation[i] = poly_eval( polynomial, roots[i] ); } cout << "Evaluation of the polynomial at the roots: " << evaluation.transpose(); }
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
Real cauchy_max_bound( const Polynomial& poly )
provides a maximum bound (the Cauchy one: ) for the absolute value of a root of the given polynomial i.e.
root of
,
The leading coefficient
: should be non zero
.
The function
Real cauchy_min_bound( const Polynomial& poly )
provides a minimum bound (the Cauchy one: ) for the absolute value of a non zero root of the given polynomial i.e.
root of
,
Computes the complex roots of a polynomial by computing the eigenvalues of the associated companion matrix with the QR algorithm.
The roots of 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 have distinct moduli i.e.
.
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 <unsupported/Eigen/Polynomials> #include <vector> #include <iostream> using namespace Eigen; using namespace std; int main() { typedef Matrix<double,5,1> Vector5d; Vector5d roots = Vector5d::Random(); cout << "Roots: " << roots.transpose() << endl; Eigen::Matrix<double,6,1> polynomial; roots_to_monicPolynomial( roots, polynomial ); PolynomialSolver<double,5> psolve( polynomial ); cout << "Complex roots: " << psolve.roots().transpose() << endl; std::vector<double> realRoots; psolve.realRoots( realRoots ); Map<Vector5d> mapRR( &realRoots[0] ); cout << "Real roots: " << mapRR.transpose() << endl; cout << endl; cout << "Illustration of the convergence problem with the QR algorithm: " << endl; cout << "---------------------------------------------------------------" << endl; Eigen::Matrix<float,7,1> hardCase_polynomial; hardCase_polynomial << -0.957, 0.9219, 0.3516, 0.9453, -0.4023, -0.5508, -0.03125; cout << "Hard case polynomial defined by floats: " << hardCase_polynomial.transpose() << endl; PolynomialSolver<float,6> psolvef( hardCase_polynomial ); cout << "Complex roots: " << psolvef.roots().transpose() << endl; Eigen::Matrix<float,6,1> evals; for( int i=0; i<6; ++i ){ evals[i] = std::abs( poly_eval( hardCase_polynomial, psolvef.roots()[i] ) ); } cout << "Norms of the evaluations of the polynomial at the roots: " << evals.transpose() << endl << endl; cout << "Using double's almost always solves the problem for small degrees: " << endl; cout << "-------------------------------------------------------------------" << endl; PolynomialSolver<double,6> psolve6d( hardCase_polynomial.cast<double>() ); cout << "Complex roots: " << psolve6d.roots().transpose() << endl; for( int i=0; i<6; ++i ) { std::complex<float> castedRoot( psolve6d.roots()[i].real(), psolve6d.roots()[i].imag() ); evals[i] = std::abs( poly_eval( hardCase_polynomial, castedRoot ) ); } cout << "Norms of the evaluations of the polynomial at the roots: " << evals.transpose() << endl << endl; cout.precision(10); cout << "The last root in float then in double: " << psolvef.roots()[5] << "\t" << psolve6d.roots()[5] << endl; std::complex<float> castedRoot( psolve6d.roots()[5].real(), psolve6d.roots()[5].imag() ); cout << "Norm of the difference: " << std::abs( psolvef.roots()[5] - castedRoot ) << endl; }
In the above example:
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.72694e-07 1.43051e-06 1.76896e-05 1.74676e-06 1.74676e-06 0.0823092 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.75127983,0) (-16.75128099,0) Norm of the difference: 1.907348633e-06
NumTraits<typename Polynomial::Scalar>::Real Eigen::cauchy_max_bound | ( | const Polynomial & | poly | ) | [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. ![]() ![]() |
Precondition: the leading coefficient of the input polynomial poly must be non zero
NumTraits<typename Polynomial::Scalar>::Real Eigen::cauchy_min_bound | ( | const Polynomial & | poly | ) | [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. ![]() ![]() |
T Eigen::poly_eval | ( | const Polynomials & | poly, | |
const T & | x | |||
) | [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. ![]() ![]() |
[in] | x | : the value to evaluate the polynomial at. |
References Eigen::poly_eval_horner().
T Eigen::poly_eval_horner | ( | const Polynomials & | poly, | |
const T & | x | |||
) | [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. ![]() ![]() |
[in] | x | : the value to evaluate the polynomial at. |
Note for stability:
Referenced by Eigen::poly_eval().
void Eigen::roots_to_monicPolynomial | ( | const RootVector & | rv, | |
Polynomial & | poly | |||
) | [inline] |
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. ![]() ![]() |