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