# source:git/Singular/LIB/solve.lib@da408f

spielwiese
Last change on this file since da408f was da408f, checked in by Moritz Wenk <wenk@…>, 24 years ago
*wenk: changed uressolve CMD_3 -> CMD_M (4) laguerre CMD2_ -> CMD_3 removed "setFloatDigits" in extra.cc fixed output (2.2e33 -> 2.2e+33) adapted solve.lib to uressolve, laguerre, extended examples git-svn-id: file:///usr/local/Singular/svn/trunk@3247 2c84dea3-7e68-4137-9b89-c4e89433aadc
• Property mode set to `100644`
File size: 7.6 KB
Line
1///////////////////////////////////////////////////////////////////////////////
2
3version="\$Id: solve.lib,v 1.12 1999-07-08 10:18:13 wenk Exp \$";
4info="
5LIBRARY: solve.lib     PROCEDURES TO SOLVE POLYNOMIAL SYSTEMS
6AUTHOR:  Moritz Wenk,  email: wenk@mathematik.uni-kl.de
7
8PROCEDURES:
9 ures_solve(i,..);      find all roots of 0-dimensional ideal i with resultants
10 mp_res_mat(i,..);      multipolynomial resultant matrix of ideal i
11 laguerre_solve(p,..);  find all roots of univariate polynom p
12 interpolate(i,j,d);    interpolate poly from evaluation points i and results j
13";
14
15///////////////////////////////////////////////////////////////////////////////
16
17proc ures_solve( ideal gls, list # )
18"USAGE:   ures_solve(i[,k,l,m]); i ideal, k,l,m integers
19         k=0: use sparse resultant matrix of Gelfand, Kapranov and Zelevinsky
20         k=1: use resultant matrix of Macaulay (k=0 is default)
21         l>0: defines precision of fractional part if ground field is Q
22         m=0,1,2: number of iterations for approximation of roots (default=2)
23ASSUME:  i is a zerodimensional ideal with
24         nvars(basering) = ncols(i) = number of vars actually occuring in i
25RETURN:  list of all (complex) roots of the polynomial system i = 0,
26         of type number if the ground field is the complex numbers,
27         of type string if the ground field is the rational or real numbers
28EXAMPLE: example ures_solve; shows an example
29"
30{
31  int typ=0;
32  int polish=2;
33  int prec=30;
34
35  if ( size(#) >= 1 )
36  {
37    typ= #[1];
38    if ( typ < 0 || typ > 1 )
39    {
40      ERROR("Valid values for second parameter k are:
41      0: use sparse Resultant (default)
42      1: use Macaulay Resultant");
43    }
44  }
45  if ( size(#) >= 2 )
46  {
47    prec= #[2];
48    if ( prec == 0 ) { prec = 30; }
49    if ( prec < 0 )
50    {
51      ERROR("Third parameter l must be positive!");
52    }
53  }
54  if ( size(#) >= 3 )
55  {
56    polish= #[3];
57    if ( polish < 0 || polish > 2 )
58    {
59      ERROR("Valid values for fourth parameter m are:
60      0,1,2: number of iterations for approximation of roots");
61    }
62  }
63
64  if ( size(#) > 3 )
65  {
66    ERROR("only three parameters allowed!");
67  }
68
69  return(uressolve(gls,typ,prec,polish));
70
71}
72example
73{
74  "EXAMPLE:";echo=2;
75  // compute the intersection points of two curves
76  ring rsq = 0,(x,y),lp;
77  ideal gls=  x2 + y2 - 10, x2 + xy + 2y2 - 16;
78  ures_solve(gls);
79  // result is a list (x,y)-coordinates as strings
80
81  // now with complex coefficient field, precision is 20 digits
82  ring rsc= (real,20,I),(x,y),lp;
83  ideal i = (2+3*I)*x2 + (0.35+I*45.0e-2)*y2 - 8, x2 + xy + (42.7)*y2;
84  list l= ures_solve(i);
85  // result is a list of (x,y)-coordinates of complex numbers
86  l;
87  // check the result
88  subst(subst(i[1],x,l[1][1]),y,l[1][2]);
89}
90///////////////////////////////////////////////////////////////////////////////
91
92proc laguerre_solve( poly f, list # )
93"USAGE:   laguerre_solve( p[,l,m]); f poly, l,m integers
94         l>0: defines precision of fractional part if ground field is Q
95         m=0,1,2: number of iterations for approximation of roots (default=2)
96ASSUME:  p is an univariate polynom
97RETURN:  list of all (complex) roots of the polynomial p;
98         of type number if the ground field is the complex numbers,
99         of type string if the ground field is the rational or real numbers
100EXAMPLE: example laguerre_solve; shows an example
101"
102{
103  int polish=2;
104  int prec=30;
105
106  if ( size(#) >= 1 )
107  {
108    prec= #[1];
109    if ( prec == 0 ) { prec = 30; }
110    if ( prec < 0 )
111    {
112      ERROR("Fisrt parameter must be positive!");
113    }
114  }
115  if ( size(#) >= 2 )
116  {
117    polish= #[2];
118    if ( polish < 0 || polish > 2 )
119    {
120      ERROR("Valid values for third parameter are:
121      0,1,2: number of iterations for approximation of roots");
122    }
123  }
124  if ( size(#) > 2 )
125  {
126    ERROR("only two parameters allowed!");
127  }
128
129  return(laguerre(f,prec,polish));
130
131}
132example
133{
134  "EXAMPLE:";echo=2;
135  // Find all roots of an univariate polynomial using Laguerre's method:
136  ring rs1= 0,(x,y),lp;
137  poly f = 15x5 + x3 + x2 - 10;
138  laguerre_solve(f);
139
140  // Now with 10 digits precision:
141  laguerre_solve(f,10);
142
143  // Now with complex coefficients, precision is 20 digits:
144  ring rsc= (real,20,I),x,lp;
145  poly f = (15.4+I*5)*x^5 + (25.0e-2+I*2)*x^3 + x2 - 10*I;
146  list l = laguerre_solve(f);
147  l;
148  // check result, value of substituted poly should be near to zero
149  subst(f,x,l[1]);
150  subst(f,x,l[2]);
151}
152///////////////////////////////////////////////////////////////////////////////
153
154proc mp_res_mat( ideal i, list # )
155"USAGE:   mp_res_mat(i[,k]); i ideal, k integer
156         k=0: sparse resultant matrix of Gelfand, Kapranov and Zelevinsky
157         k=1: resultant matrix of Macaulay (k=0 is default)
158ASSUME:  nvars(basering) = ncols(i)-1 = number of vars actually occuring in i,
159         for k=1 i must be homogeneous
160RETURN:  module representing the multipolynomial resultant matrix
161EXAMPLE: example mp_res_mat; shows an example
162"
163{
164  int typ=2;
165
166  if ( size(#) == 1 )
167  {
168    typ= #[1];
169    if ( typ < 0 || typ > 1 )
170    {
171      ERROR("Valid values for third parameter are:
172      0: sparse resultant (default)
173      1: Macaulay resultant");
174    }
175  }
176  if ( size(#) > 1 )
177  {
178    ERROR("only two parameters allowed!");
179  }
180
181  return(mpresmat(i,typ));
182
183}
184example
185{
186  "EXAMPLE:";echo=2;
187  // compute resultant matrix in ring with parameters (sparse resultant matrix)
188  ring rsq= (0,u0,u1,u2),(x1,x2),lp;
189  ideal i= u0+u1*x1+u2*x2,x1^2 + x2^2 - 10,x1^2 + x1*x2 + 2*x2^2 - 16;
190  module m = mp_res_mat(i);
191  print(m);
192  // computing sparse resultant
193  det(m);
194
195  // compute resultant matrix (Macaulay resultant matrix)
196  ring rdq= (0,u0,u1,u2),(x0,x1,x2),lp;
197  ideal h=  homog(imap(rsq,i),x0);
198  h;
199
200  module m = mp_res_mat(h,1);
201  print(m);
202  // computing Macaulay resultant (should be the same as above!)
203  det(m);
204
205  // compute numerical sparse resultant matrix
206  setring rsq;
207  ideal ir= 15+2*x1+5*x2,x1^2 + x2^2 - 10,x1^2 + x1*x2 + 2*x2^2 - 16;
208  module mn = mp_res_mat(ir);
209  print(mn);
210  // computing sparse resultant
211  det(mn);
212}
213///////////////////////////////////////////////////////////////////////////////
214
215proc interpolate( ideal p, ideal w, int d )
216"USAGE:   interpolate(p,v,d); p,v ideals, d int
217ASSUME:  ground field K is the rational numbers,
218         p and v consist of numbers of the ground filed K, p must have
219         n elements, v must have N=(d+1)^n elements where n=nvars(basering)
220         and d=deg(f) (f is the unknown polynomial in K[x1,...,xn])
221COMPUTE: polynomial f with values given by v at points p1,..,pN derived from p;
222         more precisely: consider p as point in K^n and v as N elements in K,
223         let p1,..,pN be the points in K^n obtained by evaluating all monomials
224         of degree 0,1,...,N at p in lexicographical order,
225         then the procedure computes the polynomial f satisfying f(pi) = v[i]
226RETURN:  polynomial f of degree d
227NOTE:    mainly useful for n=1, with f satisfying f(p^(i-1)) = v[i], i=1..d+1
228EXAMPLE: example interpolate; shows an example
229"
230{
231  return(vandermonde(p,w,d));
232}
233example
234{
235  "EXAMPLE:";  echo=2;
236  ring r1 = 0,(x),lp;
237  // determine f with deg(f) = 4 and
238  // v = values of f at points 3^0, 3^1, 3^2, 3^3, 3^4
239  ideal v=16,0,11376,1046880,85949136;
240  interpolate( 3, v, 4 );
241
242  ring r2 = 0,(x,y),dp;
243  // determine f with deg(f) = 3 and
244  // v = values of f at 16 points (2,3)^0=(1,1),...,(2,3)^15=(2^15,3^15)
245  // valuation point (2,3)
246  ideal p = 2,3;
247  ideal v= 1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16;
248  poly ip= interpolate( p,v,3 );
249  ip;
250  // compute poly at point 2,3, result must be 2
251  subst(subst(ip,x,2),y,3);
252  // compute poly at point 2^15,3^15, result must be 16
253  subst(subst(ip,x,2^15),y,3^15);
254}
255///////////////////////////////////////////////////////////////////////////////
Note: See TracBrowser for help on using the repository browser.