You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.

133 lines
4.4 KiB

  1. #ifndef EIGEN_POLYNOMIALS_MODULE_H
  2. #define EIGEN_POLYNOMIALS_MODULE_H
  3. #include <Eigen/Core>
  4. #include <Eigen/src/Core/util/DisableStupidWarnings.h>
  5. #include <Eigen/Eigenvalues>
  6. // Note that EIGEN_HIDE_HEAVY_CODE has to be defined per module
  7. #if (defined EIGEN_EXTERN_INSTANTIATIONS) && (EIGEN_EXTERN_INSTANTIATIONS>=2)
  8. #ifndef EIGEN_HIDE_HEAVY_CODE
  9. #define EIGEN_HIDE_HEAVY_CODE
  10. #endif
  11. #elif defined EIGEN_HIDE_HEAVY_CODE
  12. #undef EIGEN_HIDE_HEAVY_CODE
  13. #endif
  14. /** \ingroup Unsupported_modules
  15. * \defgroup Polynomials_Module Polynomials module
  16. *
  17. *
  18. *
  19. * \brief This module provides a QR based polynomial solver.
  20. *
  21. * To use this module, add
  22. * \code
  23. * #include <unsupported/Eigen/Polynomials>
  24. * \endcode
  25. * at the start of your source file.
  26. */
  27. #include "src/Polynomials/PolynomialUtils.h"
  28. #include "src/Polynomials/Companion.h"
  29. #include "src/Polynomials/PolynomialSolver.h"
  30. /**
  31. \page polynomials Polynomials defines functions for dealing with polynomials
  32. and a QR based polynomial solver.
  33. \ingroup Polynomials_Module
  34. The remainder of the page documents first the functions for evaluating, computing
  35. polynomials, computing estimates about polynomials and next the QR based polynomial
  36. solver.
  37. \section polynomialUtils convenient functions to deal with polynomials
  38. \subsection roots_to_monicPolynomial
  39. The function
  40. \code
  41. void roots_to_monicPolynomial( const RootVector& rv, Polynomial& poly )
  42. \endcode
  43. computes the coefficients \f$ a_i \f$ of
  44. \f$ p(x) = a_0 + a_{1}x + ... + a_{n-1}x^{n-1} + x^n \f$
  45. where \f$ p \f$ is known through its roots i.e. \f$ p(x) = (x-r_1)(x-r_2)...(x-r_n) \f$.
  46. \subsection poly_eval
  47. The function
  48. \code
  49. T poly_eval( const Polynomials& poly, const T& x )
  50. \endcode
  51. evaluates a polynomial at a given point using stabilized H&ouml;rner method.
  52. The following code: first computes the coefficients in the monomial basis of the monic polynomial that has the provided roots;
  53. then, it evaluates the computed polynomial, using a stabilized H&ouml;rner method.
  54. \include PolynomialUtils1.cpp
  55. Output: \verbinclude PolynomialUtils1.out
  56. \subsection Cauchy bounds
  57. The function
  58. \code
  59. Real cauchy_max_bound( const Polynomial& poly )
  60. \endcode
  61. provides a maximum bound (the Cauchy one: \f$C(p)\f$) for the absolute value of a root of the given polynomial i.e.
  62. \f$ \forall r_i \f$ root of \f$ p(x) = \sum_{k=0}^d a_k x^k \f$,
  63. \f$ |r_i| \le C(p) = \sum_{k=0}^{d} \left | \frac{a_k}{a_d} \right | \f$
  64. The leading coefficient \f$ p \f$: should be non zero \f$a_d \neq 0\f$.
  65. The function
  66. \code
  67. Real cauchy_min_bound( const Polynomial& poly )
  68. \endcode
  69. provides a minimum bound (the Cauchy one: \f$c(p)\f$) for the absolute value of a non zero root of the given polynomial i.e.
  70. \f$ \forall r_i \neq 0 \f$ root of \f$ p(x) = \sum_{k=0}^d a_k x^k \f$,
  71. \f$ |r_i| \ge c(p) = \left( \sum_{k=0}^{d} \left | \frac{a_k}{a_0} \right | \right)^{-1} \f$
  72. \section QR polynomial solver class
  73. Computes the complex roots of a polynomial by computing the eigenvalues of the associated companion matrix with the QR algorithm.
  74. The roots of \f$ p(x) = a_0 + a_1 x + a_2 x^2 + a_{3} x^3 + x^4 \f$ are the eigenvalues of
  75. \f$
  76. \left [
  77. \begin{array}{cccc}
  78. 0 & 0 & 0 & a_0 \\
  79. 1 & 0 & 0 & a_1 \\
  80. 0 & 1 & 0 & a_2 \\
  81. 0 & 0 & 1 & a_3
  82. \end{array} \right ]
  83. \f$
  84. However, the QR algorithm is not guaranteed to converge when there are several eigenvalues with same modulus.
  85. Therefore the current polynomial solver is guaranteed to provide a correct result only when the complex roots \f$r_1,r_2,...,r_d\f$ have distinct moduli i.e.
  86. \f$ \forall i,j \in [1;d],~ \| r_i \| \neq \| r_j \| \f$.
  87. With 32bit (float) floating types this problem shows up frequently.
  88. However, almost always, correct accuracy is reached even in these cases for 64bit
  89. (double) floating types and small polynomial degree (<20).
  90. \include PolynomialSolver1.cpp
  91. In the above example:
  92. -# a simple use of the polynomial solver is shown;
  93. -# the accuracy problem with the QR algorithm is presented: a polynomial with almost conjugate roots is provided to the solver.
  94. Those roots have almost same module therefore the QR algorithm failed to converge: the accuracy
  95. of the last root is bad;
  96. -# a simple way to circumvent the problem is shown: use doubles instead of floats.
  97. Output: \verbinclude PolynomialSolver1.out
  98. */
  99. #include <Eigen/src/Core/util/ReenableStupidWarnings.h>
  100. #endif // EIGEN_POLYNOMIALS_MODULE_H
  101. /* vim: set filetype=cpp et sw=2 ts=2 ai: */