source: git/Singular/LIB/solve.lib @ 0d69875

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