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 |
|
---|
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 |
|
---|
27 | /**
|
---|
28 | * \defgroup Polynomials_Module Polynomials module
|
---|
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: */
|
---|