source:git/Singular/LIB/solve.lib@d866a4

spielwiese
Last change on this file since d866a4 was d866a4, checked in by Gert-Martin Greuel <greuel@…>, 24 years ago
* Greuel/Wenk: example of laguerre_solve with complex coefficients git-svn-id: file:///usr/local/Singular/svn/trunk@3219 2c84dea3-7e68-4137-9b89-c4e89433aadc
• Property mode set to `100644`
File size: 7.4 KB
Line
1///////////////////////////////////////////////////////////////////////////////
2
3version="\$Id: solve.lib,v 1.5 1999-07-02 12:47:47 greuel 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;
97         of type number if the ground field is the complex numbers,
98         of type string if the ground field is the rational or real numbers
99EXAMPLE: example laguerre_solve; shows an example
100"
101{
102  int polish=2;
103  int prec=30;
104
105  if ( size(#) >= 1 )
106    {
107      prec= #[1];
108      if ( prec == 0 ) { prec = 30; }
109      if ( prec < 0 )
110        {
111          ERROR("Fisrt parameter must be positive!");
112        }
113    }
114  if ( size(#) >= 2 )
115    {
116      polish= #[2];
117      if ( polish < 0 || polish > 2 )
118        {
119          ERROR("Valid values for third parameter are:
120          0,1,2: number of iterations for approximation of roots");
121        }
122    }
123  if ( size(#) > 2 )
124    {
125      ERROR("only two parameters allowed!");
126    }
127
128  int digits= system("setFloatDigits",prec);
129
130  return(laguerre(f,polish));
131
132}
133example
134{
135  "EXAMPLE:";echo=2;
136  // Find all roots of an univariate polynomial using Laguerre's method:
137  ring rs1= 0,(x,y),lp;
138  poly f = 15x5 + x3 + x2 - 10;
139  laguerre_solve(f);
140
141  pause;
142  // Now with 10 digits precision:
143  laguerre_solve(f,10);
144
145  pause;
146  // Now with complex coefficients, precision is 20 digits:
147  ring rsc= (real,20,I),x,lp;
148  poly f = (15+I*5)*x^5 + (0.25+I*2)*x^3 + x2 - 10*I;
149  list l = laguerre_solve(f);
150  l;
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  pause;
196  // compute resultant matrix (Macaulay resultant matrix)
197  ring rdq= (0,u0,u1,u2),(x0,x1,x2),lp;
198  ideal h=  homog(imap(rsq,i),x0);
199  hgls;
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  pause;
206  // compute numerical sparse resultant matrix
207  setring rsq;
208  ideal ir= 15+2*x1+5*x2,x1^2 + x2^2 - 10,x1^2 + x1*x2 + 2*x2^2 - 16;
209  module mn = mp_res_mat(ir);
210  print(mn);
211  // computing sparse resultant
212  det(mn);
213}
214///////////////////////////////////////////////////////////////////////////////
215
216proc interpolate( ideal p, ideal w, int d )
217"USAGE:   interpolate(p,v,d); p,v ideals, d int
218ASSUME:  ground field K is the rational numbers,
219         p and v consist of numbers of the ground filed K, p must have
220         n elements, v must have N=(d+1)^n elements where n=nvars(basering)
221         and d=deg(f) (f is the unknown polynomial in K[x1,...,xn])
222COMPUTE: polynomial f with values given by v at points p1,..,pN derived from p;
223         more precisely: consider p as point in K^n and v as N elements in K,
224         let p1,..,pN be the points in K^n obtained by evaluating all monomials
225         of degree 0,1,...,N at p in lexicographical order,
226         then the procedure computes the polynomial f satisfying f(pi) = v[i]
227RETURN:  polynomial f of degree d
228NOTE:    mainly useful for n=1, with f satisfying f(p^(i-1)) = v[i], i=1..d+1
229EXAMPLE: example interpolate; shows an example
230"
231{
232  return(vandermonde(p,w,d));
233}
234example
235{
236  "EXAMPLE:";  echo=2;
237  ring r1 = 0,(x),lp;
238  // determine f with deg(f) = 4 and
239  // v = values of f at points 3^0, 3^1, 3^2, 3^3, 3^4
240  ideal v=16,0,11376,1046880,85949136;
241  interpolate( 3, v, 4 );
242
243  ring r2 = 0,(x,y),dp;
244  // determine f with deg(f) = 3 and
245  // v = values of f at 16 points (2,3)^0=(1,1),...,(2,3)^15=(2^15,3^15)
246  // valuation point (2,3)
247  ideal p = 2,3;
248  ideal v= 1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16;
249  interpolate( p,v,3 );
250}
251///////////////////////////////////////////////////////////////////////////////
252
253
254
255
256
257
Note: See TracBrowser for help on using the repository browser.