Read-only forum archive

Help solving a more difficult multivariate polynomial system

Help solving a more difficult multivariate polynomial system

mcody · Tue Oct 29, 2013 6:30 am

Greetings,
I have been using Singular over the last month or so to try to solve some systems of multivariate polynomial equations. I have been successful, to some extent, solving some systems, such as the following:

Code:
LIB "solve.lib";
// ** loaded /usr/share/Singular/LIB/solve.lib (13733,2010-12-06)
// ** loaded /usr/share/Singular/LIB/triang.lib (13499,2010-10-15)
// ** loaded /usr/share/Singular/LIB/elim.lib (14661,2012-03-05)
// ** loaded /usr/share/Singular/LIB/ring.lib (15322,2012-10-12)
// ** loaded /usr/share/Singular/LIB/primdec.lib (14732,2012-03-30)
// ** loaded /usr/share/Singular/LIB/absfact.lib (14191,2011-05-04)
// ** loaded /usr/share/Singular/LIB/matrix.lib (13658,2010-11-16)
// ** loaded /usr/share/Singular/LIB/nctools.lib (14246,2011-05-26)
// ** loaded /usr/share/Singular/LIB/random.lib (14661,2012-03-05)
// ** loaded /usr/share/Singular/LIB/poly.lib (14852,2012-04-30)
// ** loaded /usr/share/Singular/LIB/inout.lib (13499,2010-10-15)
// ** loaded /usr/share/Singular/LIB/general.lib (14191,2011-05-04)
ring R = real,(x1,y1,x2,y2),dp;
poly eq1 = 16.00000*x1^2+25.00000*y1^2-400.00000;
poly eq2 = 16.00000*x2^2-960.00000*x2+25.00000*y2^2+14000.00000;
poly eq3 = 32.00000*x1^2-32.00000*x1*x2+50.00000*y1^2-50.00000*y1*y2;
poly eq4 = 32.00000*x1*x2-960.00000*x1+50.00000*y1*y2-32.00000*x2^2+960.00000*x2-50.00000*y2^2;
ideal I = eq1,eq2,eq3,eq4;
def AC=solve(I,6,0,"nodisplay");
// 'solve' created a ring, in which a list SOL of numbers (the complex solutions)
// is stored.
// To access the list of complex solutions, type (if the name R was assigned
// to the return value):setring AC;
SOL;
[1]:
   [1]:
      -0.0000805975
   [2]:
      -4.0005
   [3]:
      30.000081
   [4]:
      -4.000009
[2]:
   [1]:
      1.666729
   [2]:
      3.7714
   [3]:
      28.333272
   [4]:
      -3.771226
[3]:
   [1]:
      1.666729
   [2]:
      -3.7714
   [3]:
      28.333272
   [4]:
      3.771226
[4]:
   [1]:
      -0.0000805975
   [2]:
      4.0005
   [3]:
      30.000081
   [4]:
      4.000009
[5]:
   [1]:
      15
   [2]:
      (i*11.313848)
   [3]:
      15
   [4]:
      (i*11.313709)
[6]:
   [1]:
      15
   [2]:
      (-i*11.313848)
   [3]:
      15
   [4]:
      (-i*11.313709)

The results are about as I expected, although I am surprised to see the two complex results that are actually points. An explanation for those results would be appreciated, but that is not my main questions.

What I am trying to do now is solve for a more difficult system of multivariate polynomials, which I provide the following example:
Code:
> LIB "solve.lib";
// ** loaded /usr/share/Singular/LIB/solve.lib (13733,2010-12-06)
// ** loaded /usr/share/Singular/LIB/triang.lib (13499,2010-10-15)
// ** loaded /usr/share/Singular/LIB/elim.lib (14661,2012-03-05)
// ** loaded /usr/share/Singular/LIB/ring.lib (15322,2012-10-12)
// ** loaded /usr/share/Singular/LIB/primdec.lib (14732,2012-03-30)
// ** loaded /usr/share/Singular/LIB/absfact.lib (14191,2011-05-04)
// ** loaded /usr/share/Singular/LIB/matrix.lib (13658,2010-11-16)
// ** loaded /usr/share/Singular/LIB/nctools.lib (14246,2011-05-26)
// ** loaded /usr/share/Singular/LIB/random.lib (14661,2012-03-05)
// ** loaded /usr/share/Singular/LIB/poly.lib (14852,2012-04-30)
// ** loaded /usr/share/Singular/LIB/inout.lib (13499,2010-10-15)
// ** loaded /usr/share/Singular/LIB/general.lib (14191,2011-05-04)
> ring R = real,(x1,y1,x2,y2),dp;
> poly eqn1 = 20.50000*x1^2-9.00000*x1*y1+20.50000*y1^2-400.00000;
> poly eqn2 = 22.75000*x2^2+7.79423*x2*y2-1365.00000*x2+18.25000*y2^2-233.82686*y2+20075.00000;
> poly eqn3 = 41.00000*x1^2-18.00000*x1*y1-41.00000*x1*x2+9.00000*x1*y2+41.00000*y1^2+9.00000*y1*x2-41.00000*y1*y2;
> poly eqn4 = 45.50000*x1*x2+7.79423*x1*y2-1365.00000*x1+7.79423*y1*x2+36.50000*y1*y2-233.82686*y1-45.50000*x2^2-15.58846*x2*y2+1365.00000*x2-36.50000*y2^2+233.82686*y2;
> ideal I = eqn1,eqn2,eqn3,eqn4;
> def AC=solve(I,8,0,60, "nodisplay");
   ? ideal not zero-dimensional
   ? leaving solve.lib::solve
> dim(I);
// ** I is no standard basis
2
>

This is where, with my current understanding of the mathematical concepts and of Singular, I am lost. I don't really quite understand what is meant by "zero-dimensional" versus non-zero-dimensional". Nor do I really understand how to determine the dimension of my systems of equations. Is this the Krull dimension or the Gelfand-Kirillov dimension? For that matter what is the difference between the two? Why is the ideal I "no standard basis", but dim(I) equals 2?

I would really appreciate some assistance in trying to solve this system of equations. I have searched through previous postings on this forum and read hints that perhaps only a reordering of the polynomials may be needed. I'm not exactly sure.

Regards,

Mac Cody

Re: Help solving a more difficult multivariate polynomial system

Guest · Tue Oct 29, 2013 4:06 pm

Dear mcody,

for dim you have to respect the documentaton;
how dim works it is somewhat unexpected,
but in summary, you have to compute a standard basis before calling dim to compute the ideal dimension.
Regarding to "ideal is not zerodimensional", you need a higher floating point precision (e.g. 50):
Code:
LIB "solve.lib";

ring R = (real,50),(x1,y1,x2,y2),dp;
poly eqn1 = 20.50000*x1^2-9.00000*x1*y1+20.50000*y1^2-400.00000;
poly eqn2 = 22.75000*x2^2+7.79423*x2*y2-1365.00000*x2+18.25000*y2^2-233.82686*y2+20075.00000;
poly eqn3 = 41.00000*x1^2-18.00000*x1*y1-41.00000*x1*x2+9.00000*x1*y2+41.00000*y1^2+9.00000*y1*x2-41.00000*y1*y2;
poly eqn4 = 45.50000*x1*x2+7.79423*x1*y2-1365.00000*x1+7.79423*y1*x2+36.50000*y1*y2-233.82686*y1-45.50000*x2^2-15.58846*x2*y2+1365.00000*x2-36.50000*y2^2+233.82686*y2;
ideal I = eqn1,eqn2,eqn3,eqn4;
dim( std(I) ); // ==0

ideal I1 = eqn1;

dim(std(I1));
ideal I2 = I1+eqn2;
dim(std(I2));
ideal I3 = I2+eqn3;
dim(std(I3));
ideal I4 = I3+eqn4;
dim(std(I4)); // ==0; I4==I

def AC=solve( I , 8 , 0, 60, "nodisplay");
setring AC;
SOL;


Best,

Jack

Re: Help solving a more difficult multivariate polynomial system

mcody · Wed Oct 30, 2013 3:22 am

Jack,

Thanks for your quick reply. I tried exactly as you suggested and this is what I got:
Code:
> LIB "solve.lib";
// ** loaded /usr/share/Singular/LIB/solve.lib (13733,2010-12-06)
// ** loaded /usr/share/Singular/LIB/triang.lib (13499,2010-10-15)
// ** loaded /usr/share/Singular/LIB/elim.lib (14661,2012-03-05)
// ** loaded /usr/share/Singular/LIB/ring.lib (15322,2012-10-12)
// ** loaded /usr/share/Singular/LIB/primdec.lib (14732,2012-03-30)
// ** loaded /usr/share/Singular/LIB/absfact.lib (14191,2011-05-04)
// ** loaded /usr/share/Singular/LIB/matrix.lib (13658,2010-11-16)
// ** loaded /usr/share/Singular/LIB/nctools.lib (14246,2011-05-26)
// ** loaded /usr/share/Singular/LIB/random.lib (14661,2012-03-05)
// ** loaded /usr/share/Singular/LIB/poly.lib (14852,2012-04-30)
// ** loaded /usr/share/Singular/LIB/inout.lib (13499,2010-10-15)
// ** loaded /usr/share/Singular/LIB/general.lib (14191,2011-05-04)
> ring R = (real,50),(x1,y1,x2,y2),dp;
> poly eqn1 = 20.50000*x1^2-9.00000*x1*y1+20.50000*y1^2-400.00000;
> poly eqn2 = 22.75000*x2^2+7.79423*x2*y2-1365.00000*x2+18.25000*y2^2-233.82686*y2+20075.00000;
> poly eqn3 = 41.00000*x1^2-18.00000*x1*y1-41.00000*x1*x2+9.00000*x1*y2+41.00000*y1^2+9.00000*y1*x2-41.00000*y1*y2;
> poly eqn4 = 45.50000*x1*x2+7.79423*x1*y2-1365.00000*x1+7.79423*y1*x2+36.50000*y1*y2-233.82686*y1-45.50000*x2^2-15.58846*x2*y2+1365.00000*x2-36.50000*y2^2+233.82686*y2;
> ideal I = eqn1,eqn2,eqn3,eqn4;
> dim( std(I) );
0
> ideal I1 = eqn1;
> dim(std(I1));
3
> ideal I2 = I1+eqn2;
> dim(std(I2));
2
> ideal I3 = I2+eqn3;
> dim(std(I3));
1
> ideal I4 = I3+eqn4;
> dim(std(I4));
0
> def AC=solve( I , 8 , 0, 60, "nodisplay");
   ? ideal not zero-dimensional
   ? leaving solve.lib::solve
>

I also tried solving I4 with the following results:
Code:
> def AC=solve( I4 , 8 , 0, 60, "nodisplay");
// ** redefining AC **
   ? ideal not zero-dimensional
   ? leaving solve.lib::solve
>

Did I do something wrong? Am I missing something here? Any suggestions?

Thanks again,

Mac

Re: Help solving a more difficult multivariate polynomial system

Jack · Wed Oct 30, 2013 7:41 pm

uhh-ohhh. Your example works only in the current development version, all previous Singular versions are buggy here.
Even a call
Code:
def AC=solve( std(I) , 8 , 0, 60, "nodisplay");

will not help (crash)

You have two options:
    - either to wait for the new 4.0 release (several weeks or months)
    - or you have to build singular by hand :

For experienced users/developers: you can find here
Code:
github.com/Singular/Sources/wiki/Building-Singular-from-source (it is a URL. I removed the http prefix since it is not allowed for guests to post urls )

a recipe how to do it.

If you are not a developer , the easiest way would be to ask an administrator or a friend.
Otherwise ask for furhter support.

Some hints:

install Linux packages 'git', 'automake', 'autoconf', 'libtool', 'g++', 'gmp-devel', 'ntl-devel', 'mpfr-devel', and probably others

then make a local project copy and try to build the linux binary:
Code:
git clone -b spielwiese https://github.com/Singular/Sources.git Singular-spielwiese
cd Singular-spielwiese/; ./autogen.sh
mkdir BUILD; cd BUILD
../configure
make

if something fails in the "git clone" step, google for it. usually it is a config issue
if something failes in the configure step, you probably need to install another linux package (with development stuff, package is usually called $libname-devel )

then if 'make' succeeded, look for a file "Singular/Singular"

Re: Help solving a more difficult multivariate polynomial system

Jack · Wed Oct 30, 2013 8:21 pm

if you just need the result output, here it is
(I did not check for correctness)

Code:
> SOL;
[1]:
   [1]:
      -1.0265836
   [2]:
      -4.5275683
   [3]:
      30.78896161
   [4]:
      -4.76959024
[2]:
   [1]:
      2.20823573
   [2]:
      4.34100979
   [3]:
      29.56570656
   [4]:
      -4.56465285
[3]:
   [1]:
      0.47176359
   [2]:
      -4.28965946
   [3]:
      28.12512359
   [4]:
      4.60704423
[4]:
   [1]:
      0.95678666
   [2]:
      4.5275332
   [3]:
      29.15137
   [4]:
      4.76955949
[5]:
   [1]:
      (15.79355826+i*0.19988323)
   [2]:
      (3.26332842+i*14.80559908)
   [3]:
      (15.79355826+i*0.19988323)
   [4]:
      (3.26332842+i*14.80559908)
[6]:
   [1]:
      (15.79355826-i*0.19988323)
   [2]:
      (3.26332842-i*14.80559908)
   [3]:
      (15.79355826-i*0.19988323)
   [4]:
      (3.26332842-i*14.80559908)
[7]:
   [1]:
      (33.3421346-i*66.53630777)
   [2]:
      (72.35258221+i*17.86330541)
   [3]:
      (33.3421346-i*66.53630777)
   [4]:
      (72.35258221+i*17.86330541)
[8]:
   [1]:
      (33.3421346+i*66.53630777)
   [2]:
      (72.35258221-i*17.86330541)
   [3]:
      (33.3421346+i*66.53630777)
   [4]:
      (72.35258221-i*17.86330541)

Re: Help solving a more difficult multivariate polynomial system

mcody · Thu Oct 31, 2013 6:47 am

Jack,

I downloaded, compiled, and installed the current development version, as per your directions (or so I thought - see below). I was able to successfully solve my system of equations, obtaining the same answers you listed in your follow-up post. I still find it odd (given my limited understanding) that there are four complex solutions that are effectively points. What is the cause of this or what does it mean? I wasn't expecting any complex solutions. Was was just expecting the four real-valued solutions.

The odd thing is that the configure script did not fail when it didn't find NTL and FLINT. There were tell-tale warnings in the config.log file, but I didn't check there until later. Given that configure did not explicitly fail, I blindly went ahead with compiling and installing Singular. Apparently solving my problem didn't depend on either of the missing libraries or some internal work-arounds were performed. The only complaint that occurred was when ESingular started, producing the following messge:
Code:
// ** Could not get 'InfoFile'.
// ** Either set environment variable 'SINGULAR_INFO_FILE' to 'InfoFile',
// ** or make sure that 'InfoFile' is at "/usr/local/bin/../info/singular.hlp"
// ** InfoFile:

I will go back and recompile Singular once I get the proper libraries in place. This will require upgrading my GCC, due to FLINT's dependency on newer GMP and MPFR libraries.

Many thanks for your help!

Mac

Re: Help solving a more difficult multivariate polynomial system

Jack · Thu Oct 31, 2013 11:35 am

Dear Mac,


Quote:
The odd thing is that the configure script did not fail when it didn't find NTL and FLINT.

Yes, Singular is working around missing NTL. FLINT is optional.

Regarding the output of complex numbers,
you could ask the "solve.lib" authors (see 'Singular/LIB/solve.lib' for email address )
or report the issue here
Code:
http://www.singular.uni-kl.de/index.php/singular-report-bugs.html
(clear milestone entry, set version to singular-spielwiese , priority: minor)

When looking at the documentation in the 'solve.lib' file, you will read, that it is about complex solving
of polynomial systems. I can imagine that it is simpler to do it over complex numbers than over reals.
Cumbersome fact is, that you can't pass an ideal with complex coefficients to 'solve',
due to Singular's internal implementation of polynomial rings over complex numbers...
Whatever.

I would appreciate, if you could look for some other bugs and issues in Singular you can think of
and report them. E.g. test routines you regularily use with known examples, special cases or random input
and verify the output by using other CAS(e.g. free: Sagemath, Macaulay2 , commertial: Magma, Maple, Matlab, Mathematica, ),
computation by hand or reference results for known examples.

Thanks,


Jack

Re: Help solving a more difficult multivariate polynomial system

gorzel · Thu Oct 31, 2013 7:16 pm

Quoting hannes from
viewtopic.php__q__f-10__and__t-1937__and__p-2881#p2881

Quote:
As soon as one calls groebner or std in a ring with inexact coefficients
(like real or complex), most probably the result will be unusable.


What does this mean?
Keep in mind that Singular is mainly a system for symbolic computation with polynomial expressions
over the rational numbers (or derivatives from this). With your definition solve performs a std-comoutation
over the reals which is not numerically good.

You will obtain better results when your input is an ideal with precise rational coefficients,
in fact the equations of the first example have integer coefficients,
and do a Groebner basis computation over Q[x1,y1,x2,y2] to simplify the ideal.

A call with solve to obtain numerical values should be the last step.

Code:
LIB "solve.lib";

ring R0 = 0,(x1,y1,x2,y2),dp;
ideal I =
16*x1^2+25*y1^2-400,
16*x2^2+25*y2^2-960*x2+14000,
32*x1^2+50*y1^2-32*x1*x2-50*y1*y2,
32*x1*x2-32*x2^2+50*y1*y2-50*y2^2-960*x1+960*x2;
proc displaynumsolve(ideal I)
{
option(redSB);
ideal Istd = std(I);
if (dim(Istd)>0) {ERROR("non-isolated solution");}
"There are ",vdim(std(I))," isolated solutions";

list L = facstd(std(I));

for (int n=1;n<=size(L);n++)
{
   "n = ",n;
   L[n];
   "numerical solutions: (" + string(maxideal(1)) + ") =";
  solve(L[n]);
}
}
displaynumsolve(I);


Then the result are even better in the first example:

Code:
There are  6  isolated solutions:

n =  1
_[1]=y2+4
_[2]=x2-30
_[3]=y1+4
_[4]=x1
numerical solutions: (x1,y1,x2,y2) =
[1]:
   [1]:
      0
   [2]:
      -4
   [3]:
      30
   [4]:
      -4

n =  2
_[1]=y2-4
_[2]=x2-30
_[3]=y1-4
_[4]=x1
numerical solutions: (x1,y1,x2,y2) =
[1]:
   [1]:
      0
   [2]:
      4
   [3]:
      30
   [4]:
      4

n =  3
_[1]=x2-15
_[2]=y1-y2
_[3]=x1-15
_[4]=y2^2+128
numerical solutions: (x1,y1,x2,y2) =
[1]:
   [1]:
      15
   [2]:
      (-i*11.3137085)
   [3]:
      15
   [4]:
      (-i*11.3137085)
[2]:
   [1]:
      15
   [2]:
      (i*11.3137085)
   [3]:
      15
   [4]:
      (i*11.3137085)

n =  4
_[1]=3*x2-85
_[2]=y1+y2
_[3]=3*x1-5
_[4]=9*y2^2-128
numerical solutions: (x1,y1,x2,y2) =
[1]:
   [1]:
      1.66666667
   [2]:
      3.77123617
   [3]:
      28.33333333
   [4]:
      -3.77123617
[2]:
   [1]:
      1.66666667
   [2]:
      -3.77123617
   [3]:
      28.33333333
   [4]:
      3.77123617


Try the same with your second example with exact rational coefficients.

Re: Help solving a more difficult multivariate polynomial system

mcody · Sun Nov 03, 2013 6:58 am

gorzel,

Thanks for your insightful reply!

I had originally used your approach successfully on my first example because that example indeed has integer coefficients. That encouraged me that I was headed in the "correct" direction. I was stuck, though, with my second example because it possesses real coefficients. In the end, I actually do want to obtain a symbolic solution to a specific set of multivariate polynomial equations. I just thought, due to my lack of understanding of Singular and this branch of mathematics, that it would be "simpler" to approach the problem with a numeric solution first. That appears to be my mistaken assumption. My studying of Singular indicated to me that symbolic solutions would take a long time to solve, may not be solvable due to hardware computational and storage limitations, and a solution would probably have a large number of terms.

Given the nature of what I'm trying to solve, it appears impractical, or at least difficult, to always restate the polynomial coefficients as exact rational coefficients. Perhaps it might be easier to symbolically solve my specific set of multivariate polynomial equations. I am just experiencing a difficult time doing so. I just tried (again) to reduce the polynomial set to a Groebner basis, using std(), and after tens of minutes of processing, my computer crashed. Perhaps you can enlighten me as to what to do.

Here, in all its "glory" is the system of multivariate equations that I am trying to solve. This is the code I just tried that resulted in my system crash.
Code:
LIB "solve.lib";
ring R = (0,a1,b1,c1,d1,e1,f1,a2,b2,c2,d2,e2,f2),(x1,y1,x2,y2),dp;
poly eq1 = a1*x1^2 + b1*x1*y1 + d1*x1 + c1*y1^2 + e1*y1 + f1;
poly eq2 = a2*x2^2 + b2*x2*y2 + d2*x2 + c2*y2^2 + e2*y2 + f2;
poly eq3 = 2*a1*x1^2 + 2*b1*x1*y1 - 2*a1*x1*x2 - b1*x1*y2 + d1*x1 + 2*c1*y1^2 - b1*y1*x2 - 2*c1*y1*y2 + e1*y1 - d1*x2 - e1*y2;
poly eq4 = 2*a2*x1*x2 + b2*x1*y2 + d2*x1 + b2*y1*x2 + 2*c2*y1*y2 + e2*y1 - 2*a2*x2^2 - 2*b2*x2*y2 - d2*x2 - 2*c2*y2^2 - e2*y2;
ideal I = eq1,eq2,eq3,eq4;
Istd = std(I);

The coefficients for two conic sections (actually ellipses in my need) are represented by a1,b1,c1,d1,e1, and f1 for ellipse E1, and a2,b2,c2,d2,e2, and f2 for ellipse E2. The variables x1,y1,x2 and y2 represent points on the two ellipses. I believe that this is solvable, but my lack of understanding and experience with Singular is preventing me from properly setting things up to make it happen.

Thanks for your help in giving me greater understanding of Singular!

Regards,

Mac

Re: Help solving a more difficult multivariate polynomial system

gorzel · Wed Nov 06, 2013 2:56 pm

There are two questions to be discussed but I will not give a complete answer. These comments are mostly
sketches, which will prove that just giving an equation system is not sufficient to solve a geometric problem.

1.) your old second example:
Quote:
I was stuck, though, with my second example because it possesses real coefficients.
In the end, I actually do want to obtain a symbolic solution to a specific set of multivariate polynomial equations.

Most likely the coefficients of your second example are still algebraic numbers.
Perhaps the numbers 7.79423, 15.58846, 233.82686 are the square roots of
60.75 == 243/4, 243 == 3^5, 5468 == 2^2*1367, respectively?

Then define the algebraic number field Q(\sqrt{3},\sqrt{1367})
by a common primitive element, called c. This can be by found hand,
in the general case you may use primitive.lib which also gives terms to reexpress the
original algebraic elements in terms of the new primitive element.
Then do the computation in
Code:
ring ralg =(0,c),(x1,y1,x2,y2),dp;
minpoly = ...;

The command solve will not work over these rings but you are interested in a
symbolic solution, so facstd will work and give you what you are looking for.

2. On your general form:
First I doubt that it made your system crash. I guess the computation only got stuck. Or how long let you Singular run?

Enable option(prot); After approx 2 min one only sees:
Code:
> ideal Istd = std(I);
[65535:2]2(3)s(2)sss3(3)s(4)s(5)s(7)s4(9)--s(7)s(8)
but there is no memory swell.

In general, groebner (std) computation with (a large number of) parameters will not work in practice.
The reason is that in reduction step
multivariate gcd computations for the polynomial coefficients have to be executed. This takes much time.

To get an impression what can be done (but this may not give the same) is to
compute with entire variables and then map back to parameters,
Do the computation first in prime characteric:
(If the solution is zero-dimensional one can lift it to char 0 by multiple parallel computations).
Code:
> ring R = 32003,(a1,b1,c1,d1,e1,f1,a2,b2,c2,d2,e2,f2,x1,y1,x2,y2),dp;
> poly eq1 = a1*x1^2 + b1*x1*y1 + d1*x1 + c1*y1^2 + e1*y1 + f1;
> poly eq2 = a2*x2^2 + b2*x2*y2 + d2*x2 + c2*y2^2 + e2*y2 + f2;
> poly eq3 = 2*a1*x1^2 + 2*b1*x1*y1 - 2*a1*x1*x2 - b1*x1*y2 + d1*x1 + 2*c1*y1^2 - b1*y1*x2 - 2*c1*y1*y2 + e1*y1 - d1*x2 - e1*y2;
> poly eq4 = 2*a2*x1*x2 + b2*x1*y2 + d2*x1 + b2*y1*x2 + 2*c2*y1*y2 + e2*y1 - 2*a2*x2^2 - 2*b2*x2*y2 - d2*x2 - 2*c2*y2^2 - e2*y2;
> ideal I = eq1,eq2,eq3,eq4;
> option(prot);
> ideal Istd = std(I);
[15:2]3(3)s(2)sss4s(3)s(5)s5(6)s(10)s(13)s(16)s(19)s6(20)s(24)--s(26)-s(28)s(29)s(31)----s(29)7----ss(33)s(35)s(39)s(43)-ss(46)--s(47)-----s(45)s(50)s(52)s(55)--s(57)8---------s(50)s(55)s(58)s(61)s(65)s(69)---s(68)--------------------s(52)----s9(53)s(56)s(60)s(62)----s(63)-s(67)s(71)s(75)-------------s(67)s(70)
---s-s(73)---s(71)--s10--s(73)-s(77)s(82)------s(78)----ss(83)--s(84)------s(87)s(90)-s(94)------s(92)-s(95)------s(92)s(97)s(101)s(105)s(108)s(110)s(113)-s(115)-s(116)s(122)-s(126)----11-------s(119)s(123)s(129)s(132)s(137)---ss(141)---------------s(130)s(134)-s(137)----------s(133)s(138)s(143)s(147)s(154)
s(158)------------s(151)-s(154)----------------s(140)s(142)--s(143)------s(141)--s(143)s(146)s(150)-s(154)--s(155)s(160)s(163)----s(164)---s(166)12------------------s(154)-s(158)s(163)s(166)---s(168)s(174)---------------s(162)-----s(160)s(163)-----------s(157)s(161)-s(165)-------s(161)-s(164)s(170)s(174)-----------
---------s(157)-s(159)------------------------s(139)s(143)------------s(134)-------------------13--------------s(104)----(100)-s(101)s(104)s(107)s(111)s(117)s(122)s(126)-s(130)s(133)s(137)-s(140)-s(146)s(151)-s(158)s(165)-----------s(157)---------------------------------------------------------(100)------14---------------------------
s(70)----------------------s(53)---s(54)s(57)-------s(56)----s(59)---s(62)s(67)-s(72)-s(78)--s(82)---s(86)----s15(84)--------------s(73)----------------s(61)----------s(57)s(61)s(66)----s(68)s(75)s(80)s(86)s(91)---s(92)s(100)--s(105)---s(110)--s(116)--s(120)-s(126)--s(130)--s(135)--16----------------------s(114)-------------
-(100)------s(99)s(104)s(110)--s(111)-----------(100)--------s(97)-s(104)s(110)--s(113)------s--s(117)-----------------(100)-s(106)s(114)---------s(112)s(120)------ss(126)-----17--------------------s(104)----(100)------------------s(86)s(91)----s(92)---------
----------s(81)--s(86)--s(90)--------s(89)----s(93)s(101)-(100)---s(103)---(100)---s(104)----(100)--s(106)------(100)---s(104)18----(100)-------------s(90)-----------------s(79)-----s(78)-----------------s(68)s(75)-s(81)--------s(78)--------s-s(83)s(90)-s(97)--
---s--s(100)----s(104)---s19(107)-------(100)----------s(92)------------------s(79)s(83)-----------------------s(66)----s(68)s(75)---------------s(67)------ss(75)--------------20----------------------------------s(34)------------s(27)-----s(28)------s(29)21------------------
-----------
product criterion:503 chain criterion:25916
Even
Code:
option(redSB);list L = facstd(Istd);
will work
(it takes 15 min on Intel(R) Core(TM) i5-3470 CPU @ 3.20GHz, 6MB cache)
and give a list with 6 ideals,
Anyhow, the result will be large polynomials with huge coefficients in terms of the parameters.
But will this be useful?
The first component mapped back to
Code:
ring Rp = (32003,a1,b1,c1,d1,e1,f1,a2,b2,c2,d2,e2,f2),(x1,y1,x2,y2),dp;
is:
Code:
_[1]=y1-y2
_[2]=x1-x2
_[3]=(a1*b2-b1*a2)*x2*y2+(a1*c2-c1*a2)*y2^2+(a1*d2-d1*a2)*x2+(a1*e2-e1*a2)*y2+(a1*f2-f1*a2)
_[4]=(a1*b2-b1*a2)*x2^2+(-b1*c2+c1*b2)*y2^2+(-b1*d2+d1*b2)*x2+(-b1*e2+e1*b2)*y2+(-b1*f2+f1*b2)
_[5]=(a1^3*b2*c2^2-a1^2*b1*a2*c2^2-a1^2*b1*b2^2*c2-2*a1^2*c1*a2*b2*c2+a1^2*c1*b2^3+2*a1*b1^2*a2*b2*c2+2*a1*b1*c1*a2^2*c2-2*a1*b1*c1*a2*b2^2+a1*c1^2*a2^2*b2-b1^3*a2^2*c2+b1^2*c1*a2^2*b2-b1*c1^2*a2^3)*y2^3+(2*a1^3*b2*c2*e2-a1^3*c2^2*d2-2*a1^2*b1*a2*c2*e2-a1^2*b1*b2^2*e2-2*a1^2*c1*a2*b2*e2+2*a1^2*c1*a2*c2*d2+a1^2*c1*b2^2*d2+a1^2*d1*a2*c2^2-a1^2*d1*b2^2*c2-2*a1^2*e1*a2*b2*c2+a1^2*e1*b2^3+2*a1*b1^2*a2*b2*e2+2*a1*b1*c1*a2^2*e2-2*a1*b1*c1*a2*b2*d2+2*a1*b1*d1*a2*b2*c2+2*a1*b1*e1*a2^2*c2-2*a1*b1*e1*a2*b2^2-a1*c1^2*a2^2*d2-2*a1*c1*d1*a2^2*c2+2*a1*c1*e1*a2^2*b2-b1^3*a2^2*e2+b1^2*c1*a2^2*d2-b1^2*d1*a2^2*c2+b1^2*e1*a2^2*b2-2*b1*c1*e1*a2^3+c1^2*d1*a2^3)*y2^2+(-a1^3*b2^2*f2+a1^3*b2*d2*e2-a1^3*c2*d2^2+2*a1^2*b1*a2*b2*f2-a1^2*b1*a2*d2*e2+a1^2*c1*a2*d2^2-a1^2*d1*a2*b2*e2+2*a1^2*d1*a2*c2*d2-a1^2*e1*a2*b2*d2+a1^2*f1*a2*b2^2-a1*b1^2*a2^2*f2+a1*b1*d1*a2^2*e2+a1*b1*e1*a2^2*d2-2*a1*b1*f1*a2^2*b2-2*a1*c1*d1*a2^2*d2-a1*d1^2*a2^2*c2+a1*d1*e1*a2^2*b2+b1^2*f1*a2^3-b1*d1*e1*a2^3+c1*d1^2*a2^3)*x2+(a1^3*b2*c2*f2+a1^3*b2*e2^2
-a1^3*c2*d2*e2-a1^2*b1*a2*c2*f2-a1^2*b1*a2*e2^2-a1^2*b1*b2^2*f2-a1^2*c1*a2*b2*f2+a1^2*c1*a2*d2*e2+a1^2*d1*a2*c2*e2-a1^2*d1*b2^2*e2-2*a1^2*e1*a2*b2*e2+a1^2*e1*a2*c2*d2+a1^2*e1*b2^2*d2-a1^2*f1*a2*b2*c2+a1^2*f1*b2^3+2*a1*b1^2*a2*b2*f2+a1*b1*c1*a2^2*f2+2*a1*b1*d1*a2*b2*e2+2*a1*b1*e1*a2^2*e2
-2*a1*b1*e1*a2*b2*d2+a1*b1*f1*a2^2*c2-2*a1*b1*f1*a2*b2^2-a1*c1*d1*a2^2*e2-a1*c1*e1*a2^2*d2+a1*c1*f1*a2^2*b2-a1*d1*e1*a2^2*c2+a1*e1^2*a2^2*b2-b1^3*a2^2*f2-b1^2*d1*a2^2*e2+b1^2*e1*a2^2*d2+b1^2*f1*a2^2*b2-b1*c1*f1*a2^3-b1*e1^2*a2^3+c1*d1*e1*a2^3)*y2+(a1^3*b2*e2*f2-a1^3*c2*d2*f2-a1^2*b1*a2*e2*f2+a1^2*c1*a2*d2*f2+a1^2*d1*a2*c2*f2-a1^2*d1*b2^2*f2-a1^2*e1*a2*b2*f2-a1^2*f1*a2*b2*e2+a1^2*f1*a2*c2*d2+a1^2*f1*b2^2*d2+2*a1*b1*d1*a2*b2*f2+a1*b1*e1*a2^2*f2+a1*b1*f1*a2^2*e2
-2*a1*b1*f1*a2*b2*d2-a1*c1*d1*a2^2*f2-a1*c1*f1*a2^2*d2-a1*d1*f1*a2^2*c2+a1*e1*f1*a2^2*b2-b1^2*d1*a2^2*f2+b1^2*f1*a2^2*d2-b1*e1*f1*a2^3+c1*d1*f1*a2^3)
components 2,...,5 reduce with std(L[i]) to 1. The sixth component is especially large.

3. Some general remarks:

Switching to bigger machines and buying more memory is not a resort to make a problem computational feasible.

Successful strategies are the choice of good ordering (whereas this has to be found out by experiments)
and mod p techniques (modstd.lib)

But the main questions are:

What is the problem to be solved, and
what kind of result do one wants to get?

My experience is that rethinking the problem will save time and energy.
In such way you will outperform a slow computation.
-----
Christian Gorzel

Re: Help solving a more difficult multivariate polynomial system

Jack · Wed Nov 06, 2013 6:24 pm

Regarding the example
Code:
LIB "solve.lib";
ring R = (0,a1,b1,c1,d1,e1,f1,a2,b2,c2,d2,e2,f2),(x1,y1,x2,y2),dp;
poly eq1 = a1*x1^2 + b1*x1*y1 + d1*x1 + c1*y1^2 + e1*y1 + f1;
poly eq2 = a2*x2^2 + b2*x2*y2 + d2*x2 + c2*y2^2 + e2*y2 + f2;
poly eq3 = 2*a1*x1^2 + 2*b1*x1*y1 - 2*a1*x1*x2 - b1*x1*y2 + d1*x1 + 2*c1*y1^2 - b1*y1*x2 - 2*c1*y1*y2 + e1*y1 - d1*x2 - e1*y2;
poly eq4 = 2*a2*x1*x2 + b2*x1*y2 + d2*x1 + b2*y1*x2 + 2*c2*y1*y2 + e2*y1 - 2*a2*x2^2 - 2*b2*x2*y2 - d2*x2 - 2*c2*y2^2 - e2*y2;
ideal I = eq1,eq2,eq3,eq4;
def Istd = std(I);


1. With missing NTL there is indeed a memory leak/high memory usage at least in "spielwiese" => Singular should complain and not build if NTL is missing!

2. This computation seems not to finish in "spielwiese" anyway , while it does in 3.1.6 or 3.1.1 ( don't know, if the result is reliable. )

Since the last post from Gorzel contains a long output line, it is now in this thread necessary to scroll to the right to find the input forms for a reply.

Re: Help solving a more difficult multivariate polynomial system

mcody · Thu Nov 07, 2013 7:44 am

Jack wrote:
Regarding the example
Code:
LIB "solve.lib";
ring R = (0,a1,b1,c1,d1,e1,f1,a2,b2,c2,d2,e2,f2),(x1,y1,x2,y2),dp;
poly eq1 = a1*x1^2 + b1*x1*y1 + d1*x1 + c1*y1^2 + e1*y1 + f1;
poly eq2 = a2*x2^2 + b2*x2*y2 + d2*x2 + c2*y2^2 + e2*y2 + f2;
poly eq3 = 2*a1*x1^2 + 2*b1*x1*y1 - 2*a1*x1*x2 - b1*x1*y2 + d1*x1 + 2*c1*y1^2 - b1*y1*x2 - 2*c1*y1*y2 + e1*y1 - d1*x2 - e1*y2;
poly eq4 = 2*a2*x1*x2 + b2*x1*y2 + d2*x1 + b2*y1*x2 + 2*c2*y1*y2 + e2*y1 - 2*a2*x2^2 - 2*b2*x2*y2 - d2*x2 - 2*c2*y2^2 - e2*y2;
ideal I = eq1,eq2,eq3,eq4;
def Istd = std(I);


1. With missing NTL there is indeed a memory leak/high memory usage at least in "spielwiese" => Singular should complain and not build if NTL is missing!

I have recompiled spielwiese after finding the NTL 5.4.2 package for my Kubuntu system. I have not tried the example above yet.
Quote:
2. This computation seems not to finish in "spielwiese" anyway , while it does in 3.1.6 or 3.1.1 ( don't know, if the result is reliable. )

I tried running the example above on 3.1.6 and it did not finish after over one-half hour of waiting. My computer slowed down a lot at the end and I had to kill the Singular process. My system is and AMD Athlon 64 3200+ (2200 MHz)) with 1 GB of RAM and running Kubuntu 10.04 (lucid).
Quote:
Since the last post from Gorzel contains a long output line, it is now in this thread necessary to scroll to the right to find the input forms for a reply.

Indeed! All of the information is excellent, though.

Re: Help solving a more difficult multivariate polynomial system

mcody · Thu Nov 07, 2013 7:52 am

Christian Gorzel,

I want to acknowledge and thank you for the excellent information! I am still trying to process and understand it all. I hope to continue to pursue this more tomorrow. My less capable system, which is described in my reply to Jack, is having a difficult time processing successfully. It appears to slow down significantly such that I have to kill Singular. I need a system upgrade. :oops:

Re: less capable system

Jack · Thu Nov 07, 2013 12:13 pm

Hello,


I did not observe the memory usage peak using singular 3.1.6 ,
but I guess 8GB will be sufficient.
Remark: I did the computation on a server with 24 GB RAM.

To perform the computation you could get a login at a computer center where they have machines with plenty of RAM...

By the way, did you try out other computer algebra systems?


Best,


Jack

Re: less capable system

mcody · Fri Nov 08, 2013 4:52 am

Jack wrote:
I did not observe the memory usage peak using singular 3.1.6 ,
but I guess 8GB will be sufficient.
Remark: I did the computation on a server with 24 GB RAM.

To perform the computation you could get a login at a computer center where they have machines with plenty of RAM...

I have a laptop with an AMD Turion X2 (dual core) processor with 2 GB of RAM that I could use immediately. I also have a (possibly) more powerful computer with 4 GB RAM that I have been slowly working on as a new upgrade workstation. I will have to accelerate my efforts in switching over to that machine. I have to look into a computer center otherwise.

Quote:
By the way, did you try out other computer algebra systems?

Before working with Singular, I was using a Groebner basis library I found at Matlab Central. It appears to have code errors and (likely) I wasn't using it with complete correctness. A know that the latest version of Matlab contains Groebner basis functions. It is not really possible to use it at the moment. If there are other open source algebra systems that you can recommend, I am willing to listen. From my (limited) research, I felt that Singular offered me the best shot at solving my system of equations. Of course, that is dependent upon my using Singular correctly and having sufficient computing resources.

Mac