Dear list,
Polynomial solver currently finds roots only for real polynomial
coefficients. Having Eigen::ComplexEigenSolver lets Polynomial module to
compute roots for complex polynomial coefficients. I have made necessary
modifications. I have attached corresponding patch files. Unfortunately
module cannot deduce type of the solver from partial specialization of
corresponding matrix type, such that
template<typename RealScalar> class EigenSolver {};
template<typenameRealScalar> class
EigenSolver<Eigen::Matrix<RealScalar,_Deg,_Deg> > {// EigenSolver };
template<typename RealScalar> class
EigenSolver<Eigen::Matrix<std::complex<RealScalar>,_Deg,_Deg > > {
//ComplexEigenSolver};
So I have to add additional template parameter EigenSolverType with
default value for real coefficients:
template< typename _Scalar, int _Deg , typename EigenSolverType =
EigenSolver<Matrix<_Scalar,_Deg,_Deg> > >
class PolynomialSolver : public PolynomialSolverBase<_Scalar,_Deg>{};
I have to replace in number of places in the file Companion.h (mainly
in functions balanced() and balance()) Scalar with RealScalar, where
variables are really real.
With this code I have run test against Matlab and it gave similar
results.
Test
VectorXcd roots = VectorXcd::Random(4);
VectorXcd polynomialCoefficients;
roots_to_monicPolynomial(roots, polynomialCoefficients);
cout << "roots : " << endl;
cout << setprecision(14) << roots << endl;
cout << "polynomialCoefficients : " << endl;
cout << setprecision(14) << polynomialCoefficients << endl;
PolynomialSolver<std::complex<double>, Dynamic,
ComplexEigenSolver<MatrixXcd> > psolve(polynomialCoefficients);
VectorXcd computedRoots = psolve.roots();
cout << "computedRoots : " << endl;
cout << setprecision(14) << computedRoots << endl;
for(int index = 0; index < computedRoots.size(); ++index)
{
cout << setprecision(14) << "polynom at computedRoots * 1.0 e-11:
" << poly_eval(polynomialCoefficients, computedRoots(index))*1.0e11 << endl;
}
Output:
roots :
(0.12717062898648,-0.99749748222297)
(0.61748100222785,-0.61339152195807)
(-0.04025391399884,0.17001861629078)
(0.79192480239265,-0.29941709646901)
polynomialCoefficients :
(0.091649023983597,-0.091441416775918)
(0.24020598856634,0.37401934653925)
(-0.16301627948124,-1.8544616197629)
(-1.4963225196081,1.7402874843593)
(1,0)
computedRoots :
(-0.04025391399884,0.17001861629078)
(0.79192480239265,-0.29941709646901)
(0.61748100222785,-0.61339152195807)
(0.12717062898648,-0.99749748222297)
polynom at computedRoots * 1.0 e-11:
(8.3266726846887e-006,-1.2490009027033e-005)
polynom at computedRoots * 1.0 e-11:
(4.7184478546569e-005,-1.6653345369377e-005)
polynom at computedRoots * 1.0 e-11:
(2.2204460492503e-005,-1.3877787807814e-005)
polynom at computedRoots * 1.0 e-11:
(1.5163809063462e-005,-2.7286889655471e-005)
Matlab for the same coefficients gives following roots (remember that
matlab array of coefficients and Eigen one are reversed to each other)
computedRoots :
0.127170628986480 - 0.997497482222969i
0.617481002227849 - 0.613391521958071i
0.791924802392650 - 0.299417096469009i
-0.040253913998840 + 0.170018616290780i
Best regards,
Oleg Shirokobrod