|
|
A.7.1 Solving systems of polynomial equations
Here we turn our attention to the probably most popular aspect of
the solving problem: given a system of complex polynomial equations with only
finitely many solutions, compute floating point approximations for these
solutions. This is widely considered as a task for
numerical analysis. However, due to rounding errors, purely numerical
methods are often unstable in an unpredictable way.
Therefore, in many cases, it is worth investing more computing power to derive
additional knowledge on the geometric structure of the set of solutions (not
to mention the question of how to decide whether the set of solutions is
finite or not). The symbolic-numerical approach to the solving
problem combines numerical methods with a symbolic preprocessing.
Depending on whether we want to preserve the multiplicities of the solutions
or not, possible goals for a symbolic preprocessing are
- to find another system of generators (for instance, a reduced Groebner
basis) for the ideal I generated by the polynomial equations. Alternatively,
find a system of polynomials defining an ideal which has the same radical
as I (see Computing Groebner and Standard Bases,
resp. radical).
In any case, the goal should be to find a system for which a
numerical solution can be found more easily and in a more stable way.
For systems with a large number of generators, the first step in a SINGULAR
computation could be to reduce the number of generators by applying
the interred command (see interred). Another goal might be
- to decompose the system into several smaller (or, at least, more
accessible) systems of polynomial equations. Then, the set of solutions of
the original system is obtained by taking the union of the sets of solutions
of the new systems.
Such a decomposition can be obtained in several ways: for instance, by
computing a triangular decomposition (see triang_lib) for the ideal
I, or by applying the factorizing Buchberger algorithm (see
facstd), or by computing a primary decomposition of I
(see primdec_lib).
Moreover, the equational modelling of a problem frequently causes unwanted
solutions, for instance, zero as a multiple solution. Not only for stability
reasons, one is frequently interested to get rid of those.
This can be done by computing the saturation of I with respect to
an ideal having the excess components as set of solutions (see sat).
The SINGULAR libraries solve.lib and triang.lib provide
several commands for solving systems of polynomial equations (based on a
symbolic-numerical approach via Groebner bases, resp. resultants). In
the example below, we show some of these commands at work.
| | LIB "solve.lib";
ring r=0,x(1..5),dp;
poly f0= x(1)^3+x(2)^2+x(3)^2+x(4)^2-x(5)^2;
poly f1= x(2)^3+x(1)^2+x(3)^2+x(4)^2-x(5)^2;
poly f2=x(3)^3+x(1)^2+x(2)^2+x(4)^2-x(5)^2;
poly f3=x(4)^2+x(1)^2+x(2)^2+x(3)^2-x(5)^2;
poly f4=x(5)^2+x(1)^2+x(2)^2+x(3)^2;
ideal i=f0,f1,f2,f3,f4;
ideal si=std(i);
//
// dimension of a solution set (here: 0) can be read from a Groebner bases
// (with respect to any global monomial ordering)
dim(si);
==> 0
//
// the number of complex solutions (counted with multiplicities) is:
vdim(si);
==> 108
//
// The given system has a multiple solution at the origin. We use facstd
// to compute equations for the non-zero solutions:
option(redSB);
ideal maxI=maxideal(1);
ideal j=sat(si,maxI)[1]; // output is Groebner basis
vdim(j); // number of non-zero solutions (with mult's)
==> 76
//
// We compute a triangular decomposition for the ideal I. This requires first
// the computation of a lexicographic Groebner basis (we use the FGLM
// conversion algorithm):
ring R=0,x(1..5),lp;
ideal j=fglm(r,j);
list L=triangMH(j);
size(L); // number of triangular components
==> 7
L[1]; // the first component
==> _[1]=x(5)^2+1
==> _[2]=x(4)^2+2
==> _[3]=x(3)-1
==> _[4]=x(2)^2
==> _[5]=x(1)^2
//
// We compute floating point approximations for the solutions (with 30 digits)
def S=triang_solve(L,30);
==>
==> // 'triang_solve' created a ring, in which a list rlist of numbers (the
==> // complex solutions) is stored.
==> // To access the list of complex solutions, type (if the name R was assig\
ned
==> // to the return value):
==> setring R; rlist;
setring S;
size(rlist); // number of different non-zero solutions
==> 28
rlist[1]; // the first solution
==> [1]:
==> 0
==> [2]:
==> 0
==> [3]:
==> 1
==> [4]:
==> (-I*1.41421356237309504880168872421)
==> [5]:
==> -I
//
// Alternatively, we could have applied directly the solve command:
setring r;
def T=solve(i,30,1,"nodisplay"); // compute all solutions with mult's
==>
==> // 'solve' created a ring, in which a list SOL of numbers (the complex so\
lutions)
==> // is stored.
==> // To access the list of complex solutions, type (if the name R was assig\
ned
==> // to the return value):
==> setring R; SOL;
setring T;
size(SOL); // number of different solutions
==> 4
SOL[1][1]; SOL[1][2]; // first solution and its multiplicity
==> [1]:
==> [1]:
==> 1
==> [2]:
==> 1
==> [3]:
==> 1
==> [4]:
==> (i*2.449489742783178098197284074706)
==> [5]:
==> (i*1.732050807568877293527446341506)
==> [2]:
==> [1]:
==> 1
==> [2]:
==> 1
==> [3]:
==> 1
==> [4]:
==> (-i*2.449489742783178098197284074706)
==> [5]:
==> (i*1.732050807568877293527446341506)
==> [3]:
==> [1]:
==> 1
==> [2]:
==> 1
==> [3]:
==> 1
==> [4]:
==> (i*2.449489742783178098197284074706)
==> [5]:
==> (-i*1.732050807568877293527446341506)
==> [4]:
==> [1]:
==> 1
==> [2]:
==> 1
==> [3]:
==> 1
==> [4]:
==> (-i*2.449489742783178098197284074706)
==> [5]:
==> (-i*1.732050807568877293527446341506)
==> 1
SOL[size(SOL)]; // solutions of highest multiplicity
==> [1]:
==> [1]:
==> [1]:
==> 0
==> [2]:
==> 0
==> [3]:
==> 0
==> [4]:
==> 0
==> [5]:
==> 0
==> [2]:
==> 32
//
// Or, we could remove the multiplicities first, by computing the
// radical:
setring r;
ideal k=std(radical(i));
vdim(k); // number of different complex solutions
==> 29
def T1=solve(k,30,"nodisplay"); // compute all solutions with mult's
==>
==> // 'solve' created a ring, in which a list SOL of numbers (the complex so\
lutions)
==> // is stored.
==> // To access the list of complex solutions, type (if the name R was assig\
ned
==> // to the return value):
==> setring R; SOL;
setring T1;
size(SOL); // number of different solutions
==> 29
SOL[1];
==> [1]:
==> 1
==> [2]:
==> 1
==> [3]:
==> 1
==> [4]:
==> (-i*2.449489742783178098197284074706)
==> [5]:
==> (-i*1.732050807568877293527446341506)
|
|
Online Manual
