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.

138 lines
4.7 KiB

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