This section describes an algorithm, which uses generalized eigenvalues and eigenvectors, for solving a system of two bivariate polynomials , with
. It treats multiple roots, which are usually considered as obstacles for numerical linear algebra techniques. The numerical difficulties are handled with the help of singular value decomposition (SVD) of the matrix of eigenvectors associated to each eigenvalue. The extension to
equations could be treated similarly.
A vector is called a {generalized eigenvector} associated to the eigenvalue
if
. The matrices
and
have
generalized eigenvalues if and only if
. If
, then
can be finite, empty, or infinite. Note that if
then
. Moreover, if
is invertible then
, which is the ordinary spectrum of
.\
Recall that an matrix
with polynomial coefficients can be equivalently written as a polynomial with
matrix coefficients. If
, we obtain
, where
are
matrices.
Definition: The companion matrices of ![]() ![]() ![]()
|
And we have the following interesting property :
Proposition: With the above notation, the following equivalence holds:
|
Hereafter we will consider systems having only a {finite} number of roots. This condition is not quite restrictive since it is sufficient (and necessary) to require that polynomials
and
are coprime in
(otherwise one can divide them by their
).
In the case where one of the curve is a line, i.e. one of the polynomial has degree 1, the problem is reduced to solving a univariate polynomial. Indeed, one may assume e.g. that and thus we are looking for the roots of
. Let us denote them by
. Then we know that
where the 's are assumed to be distinct and
is a non-zero constant in
. The integer
, for all
, is called the {multiplicity} of the root
, and it turns out that
if
do not vanish at infinity in the direction of the
-axis (in other words, if the homogeneous part of
of highest degree do not vanish when
). This later condition can be avoid in the projective setting: let
be the homogeneous polynomial obtained by homogenizing
with the new variable
, then we have
where is an integer corresponding to the multiplicity of the root at infinity, and
. Moreover, it turns out that the roots
and their corresponding multiplicities can be computed by eigenvalues and eigenvectors computation.
In the sequel we will generalize this approach to the case where and
are bivariate polynomials of arbitrary degree. For this we first need to recall the notion of multiplicity in this context. Then we will show how to recover the roots from multiplication maps.
A modern definition of the intersection multiplicity of and
at a point
is (see [Ful84])
where the point is here abusively (but usually) identified with its corresponding prime ideal in
,
denoting the ordinary localization of the ring
by this prime ideal. As a result, the finite
-algebra
(which is actually finite if and only if
and
are coprime in
) can be decomposed as the direct sum
and consequently
The intersection multiplicities can be computed using a resultant. The main idea is to project ``algebraically'' the intersection points on the -axis. To do this let us see both polynomials
and
in
where
, that is to say as univariate polynomials in
with coefficients in the ring
which is a domain; we can rewrite
Their resultant (with respect to ) is an element in
which is non-zero, since
and
are assumed to be coprime in
, and which can be factorized, assuming that
and
do not vanish simultaneously, as
where the 's are distinct elements in
and
(as sets). For instance, if all the
are distinct then we have
and
for all
. Moreover, we have (see e.g. [Ful84] (S)):
Proposition: For all ![]() ![]() ![]() ![]() ![]() ![]() |
As a corollary, if all the 's are distinct (this can be easily obtained by a linear change of coordinates
) then
is nothing but the valuation of
at
. Another corollary of this result is the well-know Bezout theorem for algebraic plane curves, which is better stated in the projective context. Let us denote by
and
the homogeneous polynomials in
obtained from
and
by homogenization with the new variable
.
Proposition: [Bezout theorem] If ![]() ![]() ![]() ![]() ![]() ![]() |
{proof} It follows straightforwardly from the previous proposition and the well-known fact that is a homogeneous polynomial in
of degree
(see e.g. [Lang02], [VWe48II]). {proof}
Suppose given such a point and assume that
and
are not both zero. If
is of dimension one, then it is easy to check that there is only one point
such that
and moreover any element in
is a multiple of the vector
. Therefore to compute
, we only need to compute a basis of
and compute the quotient of its second coordinate by its first one. In case of multiple roots, this construction can be generalized as follows:
In this theorem, observe that the construction of and
is always possible because we assumed that both polynomials
and
depend on the variable
, which implies that the number of rows in the matrix
is always strictly greater than the
.
The previous theorem shows how to recover the -coordinates of the roots of the system
from the Sylvester matrix. However, instead of the Sylvester matrix we can use the Bezout matrix to compute
(see section [*] ). This matrix have all the required properties to replace the Sylvester matrix in all the previous results, except in proposition [*] . Indeed, the Bezout matrix being smaller than the Sylvester matrix, the condition given in remark [*] is not always fulfilled; using Bezout matrix in proposition [*] requires that for all roots
of
,
This condition may failed: take for instance in the system
Note that the use of the Bezout matrix gives, in practice, a faster method because it is a smaller matrix than the Sylvester matrix, even if its computation takes more time (the save in the eigen-computations is greater).
For difficulties ``c'', it is necessary to use the singular value decomposition: If a generalized eigenvalue of
and
, and
the matrix of the associated eigenvectors, then the singular value decomposition of
is
and is the rank of
. Let
The matrix is also a matrix of eigenvectors, because the changes are done on the columns, and V is invertible. Hence
and
have the same rank
, and
is a matrix of linearly independent eigenvectors. Consequently
can be used to build
and
.
This leads to the following complete algorithm:
{Algorithm { {Numerical approximation of the root of a square bivariate system} {{-0.3cm}}
{Input:} Two bivariate polynomials
template<class C> Seq<VectDse<C> > solve(p,q,EigenRes<C>());
synaps/resultant/bivariate_res.h
.#include <synaps/msolve/projection/EigenRes.h> #ifdef SYNAPS_WITH_LAPACK using std::cout; using std::cin; using std::endl; typedef MPol<double> mpol_t; int main(int argc, char **argv) { mpol_t p1("16*x0^4-160*x0^2*x1^2+400*x1^4-32*x0^3+160*x0*x1^2+6*x0^2-50*x1^2+10*x0+1.5625"), p2("x0^2-x0*x1+x1^2-1.2*x0-0.0625"); cout<<"p1="<<p1<<endl; cout<<"p2="<<p2<<endl; Seq<VectDse<double> > res; Seq<int> mux,muy; res= solve(p1,p2, EigenRes<double>(1e-3)); res= solve(p1,p2,mux, EigenRes<double>(1e-6)); res= solve(p1,p2,mux,muy, EigenRes<double>(1e-7)); cout<<"Intersection points:"<<endl; for(unsigned i=0;i<res.size();++i) std::cout<<"("<<res[i][0]<<","<<res[i][1]<<") of multiplicity "<<mux[i]<<" "<<muy[i]<<std::endl; } #else #include <iostream> int main() { std::cout <<" LAPACK is not connected"<<std::endl; return(77); } #endif //SYNAPS_WITH_LAPACK
![]() |