Changeset a5f564 in git
 Timestamp:
 Dec 16, 2000, 2:14:43 AM (22 years ago)
 Branches:
 (u'jengelhdatetime', 'ceac47cbc86fe4a15902392bdbb9bd2ae0ea02c6')(u'spielwiese', '00e2e9c41af3fde1273eb3633f4c0c7c3db2579d')
 Children:
 e35f0bd56e7ac915c46b0f6708641f7d5d3e6706
 Parents:
 4508ce5e48c7f51f11a43264563eaabbb7992df1
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

Singular/LIB/ntsolve.lib
r4508ce5 ra5f564 1 /////////////////////////////////////////////////////////////////////////////// 2 version="$Id: ntsolve.lib,v 1.7 20001215 12:08:46 Singular Exp $"; 1 //(GMG, last modified 16.12.00) 2 /////////////////////////////////////////////////////////////////////////////// 3 version="$Id: ntsolve.lib,v 1.8 20001216 01:14:43 greuel Exp $"; 3 4 info=" 4 LIBRARY: ntsolve.lib ONE REAL SOLUTION OF POLYNOMIAL SYSTEMS (NEWTON ITERATION)5 AUTHORS: Wilfred Pohl, 5 LIBRARY: ntsolve.lib REAL NEWTON SOLVING OF POLYNOMIAL SYSTEMS 6 AUTHORS: Wilfred Pohl, email: pohl@mathematik.unikl.de 6 7 Dietmar Hillebrand 7 8 8 9 PROCEDURES: 9 nt_solve(G,..); find one real root of 0dimensional ideal G10 triMNewton(G,..); find one real root for 0dim.triangular system G10 nt_solve(G,..); find one real root of 0dimensional ideal G 11 triMNewton(G,..); find one real root for 0dim triangular system G 11 12 "; 12 /////////////////////////////////////////////////////////////////////////////// 13 14 proc nt_solve( ideal gls, ideal ini, intvec ipar ) 15 "USAGE: nt_solve(gls,ini,ipar); 16 gls: the equations 17 ini: the ideal of initial values 18 ipar: control 19 ipar[1]  max. number of iterations 20 ipar[2]  accuracy, have the l2norm . 21 for the initial error eps0 = gls(ini) 22 accept solution sol with 23 gls(sol) < eps0*(0.1^ipar[2]) 24 ipar[3]  some output for contol if != 0 25 defaults  100, 10, 0 26 ASSUME: gls is a zerodimensional ideal with 27 nvars(basering) = size(gls) (> 1) 28 RETURN: ideal of one solution (if found) 29 0 (else) 13 14 /////////////////////////////////////////////////////////////////////////////// 15 16 proc nt_solve (ideal gls, ideal ini, list #) 17 "USAGE: nt_solve(gls,ini[,ipar]); gls ideal, ini ideal, ipar list or intvec 18 gls: contains the equations, for which a solution will be computed 19 ini: ideal of initial values (approximate solutions to start with) 20 ipar: control integers (default: ipar = 100,10) 21 ipar[1]  max. number of iterations 22 ipar[2]  accuracy (we have the l_2norm .): accept solution 23 sol if gls(sol) < eps0*(0.1^ipar[2]) where eps0 = gls(ini) 24 is the initial error 25 ASSUME: gls is a zerodimensional ideal with nvars(basering) = size(gls) (>1) 26 RETURN: ideal, coordinates of one solution (if found), 0 else 27 NOTE: printlevel >0: displays comments (default =0) 30 28 EXAMPLE: example nt_solve; shows an example 31 29 " … … 34 32 int di = size(gls); 35 33 if (nvars(basering) != di){ 36 ERROR(" wrong dimension");}34 ERROR("// wrong number of equations");} 37 35 if (size(ini) != di){ 38 ERROR(" wrong number of initial values");}36 ERROR("// wrong number of initial values");} 39 37 int prec = system("getPrecDigits"); // precision 40 38 41 39 int i1,i2,i3; 42 i1 = size(ipar); 43 int itmax, acc, prot; 40 int itmax, acc; 41 intvec ipar; 42 if ( size(#)>0 ){ 43 i1=1; 44 if (typeof(#[1])=="intvec") {ipar=#[1];} 45 if (typeof(#[1])=="int") {ipar[1]=#[1];} 46 if ( size(#)>1 ){ 47 i1=2; 48 if (typeof(#[2])=="int") {ipar[2]=#[2];} 49 } 50 } 51 52 int prot = printlevelvoice+2; // prot=printlevel (default:prot=0) 44 53 if (i1 < 1){itmax = 100;}else{itmax = ipar[1];} 45 54 if (i1 < 2){acc = prec/2;}else{acc = ipar[2];} 46 if (i1 < 3){prot = 0;}else{prot = ipar[3];}47 55 if ((acc <= 0)(acc > prec1)){acc = prec1;} 48 56 49 57 int dpl = di+1; 50 string out; // for prot != 0 and more58 string out; 51 59 out = "ring rnewton=(real,prec),("+varstr(basering)+"),(c,dp);"; 52 60 execute(out); 61 53 62 ideal gls1=imap(rn,gls); 54 63 module nt,sub; … … 58 67 if (i1>0){ 59 68 setring rn; kill rnewton; 60 ERROR("one var not in equation");} 69 ERROR("// one var not in equation");} 70 61 71 list direction; 62 72 ideal ini1; … … 66 76 dum = genau; 67 77 genau = genau^acc; 78 68 79 for (i1=di;i1>0;i1){ 69 80 sub[i1]=sub[i1]+gls1[i1]*gen(dpl);} … … 71 82 for (i1=di;i1>0;i1){ 72 83 nt = subst(nt,var(i1),ini1[i1]);} 84 73 85 // now we have in sub the general structure 74 86 // and in nt the structure with subst. vars … … 76 88 // compute initial error 77 89 y1 = ml2norm(nt,genau); 78 if(prot){out=" initial error = "+string(y1);out;}90 dbprint(prot,"// initial error = "+string(y1)); 79 91 y2 = genau*y1; 80 92 81 93 // begin of iteration 82 94 for(i3=1;i3<=itmax;i3++){ 83 if(prot){out=" Nr. "+string(i3);out;}95 dbprint(prot,"// iteration: "+string(i3)); 84 96 85 97 // find newton direction … … 95 107 } 96 108 } 97 if(prot){out=" dumping = "+string(dum);out;}109 dbprint(prot,"// dumping = "+string(dum)); 98 110 99 111 // new value … … 104 116 nt = subst(nt,var(i1),ini1[i1]);} 105 117 y1 = ml2norm(nt,genau); 106 if(prot){out=" error = "+string(y1);out;}118 dbprint(prot,"// error = "+string(y1)); 107 119 if(y1<y2){break;} // we are ready 108 120 } 109 121 110 122 if (y1>y2){ 111 " WARNING: no convergence";}123 "// ** WARNING: iteration bound reached with error > error bound!";} 112 124 setring rn; 113 125 ini = imap(rnewton,ini1); … … 121 133 ideal gls = x2+y2+z210, y2+z3+w8, xy+yz+xz+w5  1,w3+y; 122 134 ideal ini = 3.1,2.9,1.1,0.5; 123 intvec ipar = 200,0 ,1;135 intvec ipar = 200,0; 124 136 ideal sol = nt_solve(gls,ini,ipar); 125 137 sol; … … 234 246 } 235 247 /////////////////////////////////////////////////////////////////////////////// 236 237 ///////////////////////////////////////////////////////////////////////////////238 // solver.lib //239 // algorithms for solving algebraic system of dimension zero //240 // written by Dietmar Hillebrand and ... //241 // //242 ///////////////////////////////////////////////////////////////////////////////243 244 245 246 ///////////////////////////////////////////////////////////////////////////////247 248 // 248 249 // Multivariate Newton for triangular systems 249 // 250 // 251 /////////////////////////////////////////////////////////////////////////////// 252 253 proc triMNewton (ideal G, list a, number err, int itb) 254 "USAGE: triMNewtion(G,a,err,itb); 255 ideal G=g1,..,gn, a triangular system 256 in n vars, i.e gi=gi(var(ni+1),..,var(n)); 257 list of numbers a, an approximation of a common zero of G, 258 (a[i] to be substituted in var(i)); 259 number err, an error bound; 260 int itb, the maximal number of iterations performed. 261 RETURN: an improved approximation for a common zero of G; 250 // algorithms for solving algebraic system of dimension zero 251 // written by Dietmar Hillebrand 252 /////////////////////////////////////////////////////////////////////////////// 253 254 proc triMNewton (ideal G, ideal a, list #) 255 "USAGE: triMNewton(G,a[,ipar]); G ideal, a ideal, ipar list or intvec 256 ASSUME: G: g1,..,gn, a triangular system of n equations in n vars, 257 i.e gi=gi(var(ni+1),..,var(n)) 258 a: ideal of numbers, coordinates of an approximation of a common 259 zero of G to start with (with a[i] to be substituted in var(i)) 260 ipar: control integer vector (default: ipar = 100,10) 261 ipar[1]  max. number of iterations 262 ipar[2]  accuracy (we have as norm . absolute value ): 263 accept solution sol if G(sol) < G(a)*(0.1^ipar[2]) 264 itb: int, the maximal number of iterations performed 265 err: number, an error bound 266 RETURN: an ideal, coordinates of a better approximation of a zero of G 262 267 EXAMPLE: example triMNewton; shows an example 263 268 " 264 269 { 270 int prot = printlevel; 271 int i1,i2,i3; 272 intvec ipar; 273 if ( size(#)>0 ){ 274 i1=1; 275 if (typeof(#[1])=="intvec") {ipar=#[1];} 276 if (typeof(#[1])=="int") {ipar[1]=#[1];} 277 if ( size(#)>1 ){ 278 i1=2; 279 if (typeof(#[2])=="int") {ipar[2]=#[2];} 280 } 281 } 282 int itb, err; 283 if (i1 < 1) {itb = 100;} else {itb = ipar[1];} 284 if (i1 < 2) {err = 10;} else {err = ipar[2];} 285 265 286 if (itb == 0) 266 287 { 267 dbprint("iteration bound performed!");268 288 dbprint(prot,"// ** iteration bound reached with error > error bound!"); 289 return(a); 269 290 } 270 291 … … 289 310 if (J[i,ni+1] == 0) 290 311 { 291 print("Error: ideal not radical!");312 ERROR("// ideal not radical"); 292 313 return(); 293 314 } … … 304 325 for (i = 1; i <= n; i++) 305 326 { 306 if ( abs(h[i]) > err)327 if ( abs(h[i]) > (1/10)^err) 307 328 { 308 329 fertig = 0; … … 312 333 if ( not fertig ) 313 334 { 335 if (prot > 0) 336 { 337 "// error:"; print(abs(h[i])); 338 "// iterations to be performed: "+string(itb); 339 } 314 340 for (i = 1; i <= n; i++) 315 341 { 316 342 a[i] = a[i] + h[i]; 317 343 } 318 return(triMNewton(G,a,err,itb1)); 344 ipar = itb1,err; 345 return(triMNewton(G,a,ipar)); 319 346 } 320 347 else … … 325 352 example 326 353 { "EXAMPLE:"; echo = 2; 327 ring r=real,(z,y,x),(lp); 328 ideal i=x^21,y^2+x43,z2y4+x1; 329 list a=2,3,4; 330 number e=1.0e10; 331 332 list l = triMNewton(i,a,e,20); 354 ring r = (real,30),(z,y,x),(lp); 355 ideal i = x^21,y^2+x43,z2y4+x1; 356 ideal a = 2,3,4; 357 intvec e = 20,10; 358 ideal l = triMNewton(i,a,e); 333 359 l; 334 360 } 335 336 337 //////////////////////////////////////////////////////////////////////////////// 361 /////////////////////////////////////////////////////////////////////////////// 338 362 339 363 static proc abs( number r) … … 348 372 } 349 373 } 350 374 /////////////////////////////////////////////////////////////////////////////// 375
Note: See TracChangeset
for help on using the changeset viewer.