A.8.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).

• 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.

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); // 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.44948974278317809819728407471) ==> [5]: ==> (i*1.73205080756887729352744634151) ==> [2]: ==> [1]: ==> 1 ==> [2]: ==> 1 ==> [3]: ==> 1 ==> [4]: ==> (-i*2.44948974278317809819728407471) ==> [5]: ==> (i*1.73205080756887729352744634151) ==> [3]: ==> [1]: ==> 1 ==> [2]: ==> 1 ==> [3]: ==> 1 ==> [4]: ==> (i*2.44948974278317809819728407471) ==> [5]: ==> (-i*1.73205080756887729352744634151) ==> [4]: ==> [1]: ==> 1 ==> [2]: ==> 1 ==> [3]: ==> 1 ==> [4]: ==> (-i*2.44948974278317809819728407471) ==> [5]: ==> (-i*1.73205080756887729352744634151) ==> 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.44948974278317809819728407471) ==> [5]: ==> (-i*1.73205080756887729352744634151)