//(GMG, last modified 16.12.00) /////////////////////////////////////////////////////////////////////////////// version="$Id: ntsolve.lib,v 1.16 2009-04-06 09:17:01 seelisch Exp $"; category="Symbolic-numerical solving"; info=" LIBRARY: ntsolve.lib Real Newton Solving of Polynomial Systems AUTHORS: Wilfred Pohl, email: pohl@mathematik.uni-kl.de Dietmar Hillebrand PROCEDURES: nt_solve(G,ini,[..]); find one real root of 0-dimensional ideal G triMNewton(G,a,[..]); find one real root for 0-dim triangular system G "; LIB "general.lib"; /////////////////////////////////////////////////////////////////////////////// proc nt_solve (ideal gls, ideal ini, list #) "USAGE: nt_solve(gls,ini[,ipar]); gls,ini= ideals, ipar=list/intvec,@* gls: contains the equations, for which a solution will be computed@* ini: ideal of initial values (approximate solutions to start with),@* ipar: control integers (default: ipar = [100, 10]) @format ipar[1]: max. number of iterations ipar[2]: accuracy (we have the l_2-norm ||.||): accepts solution @code{sol} if ||gls(sol)|| < eps0*(0.1^ipar[2]) where eps0 = ||gls(ini)|| is the initial error @end format ASSUME: gls is a zerodimensional ideal with nvars(basering) = size(gls) (>1) RETURN: ideal, coordinates of one solution (if found), 0 else NOTE: if printlevel >0: displays comments (default =0) EXAMPLE: example nt_solve; shows an example " { def rn = basering; int di = size(gls); if (nvars(basering) != di){ ERROR("// wrong number of equations");} if (size(ini) != di){ ERROR("// wrong number of initial values");} int prec = system("getPrecDigits"); // precision int i1,i2,i3; int itmax, acc; intvec ipar; if ( size(#)>0 ){ i1=1; if (typeof(#[1])=="intvec") {ipar=#[1];} if (typeof(#[1])=="int") {ipar[1]=#[1];} if ( size(#)>1 ){ i1=2; if (typeof(#[2])=="int") {ipar[2]=#[2];} } } int prot = printlevel-voice+2; // prot=printlevel (default:prot=0) if (i1 < 1){itmax = 100;}else{itmax = ipar[1];} if (i1 < 2){acc = prec/2;}else{acc = ipar[2];} if ((acc <= 0)||(acc > prec-1)){acc = prec-1;} int dpl = di+1; string out; out = "ring rnewton=(real,prec),("+varstr(basering)+"),(c,dp);"; execute(out); ideal gls1=imap(rn,gls); module nt,sub; sub = transpose(jacob(gls1)); for (i1=di;i1>0;i1--){ if(sub[i1]==0){break;}} if (i1>0){ setring rn; kill rnewton; ERROR("// one var not in equation");} list direction; ideal ini1; ini1 = imap(rn,ini); number dum,y1,y2,y3,genau; genau = 0.1; dum = genau; genau = genau^acc; for (i1=di;i1>0;i1--){ sub[i1]=sub[i1]+gls1[i1]*gen(dpl);} nt = sub; for (i1=di;i1>0;i1--){ nt = subst(nt,var(i1),ini1[i1]);} // now we have in sub the general structure // and in nt the structure with subst. vars // compute initial error y1 = ml2norm(nt,genau); dbprint(prot,"// initial error = "+string(y1)); y2 = genau*y1; // begin of iteration for(i3=1;i3<=itmax;i3++){ dbprint(prot,"// iteration: "+string(i3)); // find newton direction direction=bareiss(nt,1,-1); // find dumping dum = linesearch(gls1,ini1,direction[1],y1,dum,genau); if (i3%5 == 0) { if (dum <= 0.000001) { dum = 1.0; } } dbprint(prot,"// dumping = "+string(dum)); // new value for(i1=di;i1>0;i1--){ ini1[i1]=ini1[i1]-dum*direction[1][i1];} nt = sub; for (i1=di;i1>0;i1--){ nt = subst(nt,var(i1),ini1[i1]);} y1 = ml2norm(nt,genau); dbprint(prot,"// error = "+string(y1)); if(y1y2){ "// ** WARNING: iteration bound reached with error > error bound!";} setring rn; ini = imap(rnewton,ini1); kill rnewton; return(ini); } example { "EXAMPLE:";echo=2; ring rsq = (real,40),(x,y,z,w),lp; ideal gls = x2+y2+z2-10, y2+z3+w-8, xy+yz+xz+w5 - 1,w3+y; ideal ini = 3.1,2.9,1.1,0.5; intvec ipar = 200,0; ideal sol = nt_solve(gls,ini,ipar); sol; } /////////////////////////////////////////////////////////////////////////////// static proc sqrt (number wr, number wa, number wg) { number es,we; number wb=wa; number wf=wb*wb-wr; if(wf>0){ es=wf;} else{ es=-wf;} we=wg*es; while (es>we) { wf=wf/(wb+wb); wb=wb-wf; wf=wb*wb-wr; if(wf>0){ es=wf;} else{ es=-wf;} } return(wb); } static proc il2norm (ideal H, number wg) { number wa,wb; int wi,dpl; wa = leadcoef(H[1]); wa = wa*wa; for(wi=size(H);wi>1;wi--) { wb=leadcoef(H[wi]); wa=wa+wb*wb; } return(sqrt(wa,wa,wg)); } static proc ml2norm (module H, number wg) { number wa,wb; int wi,dpl; dpl = size(H)+1; wa = leadcoef(H[1][dpl]); wa = wa*wa; for(wi=size(H);wi>1;wi--) { wb=leadcoef(H[wi][dpl]); wa=wa+wb*wb; } return(sqrt(wa,wa,wg)); } static proc linesearch(ideal nl, ideal aa, ideal bb, number z1, number tt, number gg) { int ii,d; ideal cc,jn; number ss,z2,z3,mm; mm=0.000001; ss=tt; d=size(nl); cc=aa; for(ii=d;ii>0;ii--){cc[ii]=cc[ii]-ss*bb[ii];} jn=nl; for(ii=d;ii>0;ii--){jn=subst(jn,var(ii),cc[ii]);} z2=il2norm(jn,gg); z3=-1; while(z2>=z1) { ss=0.5*ss; if(ss0;ii--) { cc[ii]=cc[ii]-ss*bb[ii]; } jn=nl; for(ii=d;ii>0;ii--){jn=subst(jn,var(ii),cc[ii]);} z3=z2; z2=il2norm(jn,gg); } if(z3<0) { while(z30;ii--) { cc[ii]=cc[ii]-ss*bb[ii]; } jn=nl; for(ii=d;ii>0;ii--){jn=subst(jn,var(ii),cc[ii]);} if(z3>0){z2=z3;} z3=il2norm(jn,gg); } } z2=z2-z1; z3=z3-z1; ss=0.25*ss*(z3-4*z2)/(z3-2*z2); if(ss>1.0){return (1.0);} if(ss0 ){ i1=1; if (typeof(#[1])=="intvec") {ipar=#[1];} if (typeof(#[1])=="int") {ipar[1]=#[1];} if ( size(#)>1 ){ i1=2; if (typeof(#[2])=="int") {ipar[2]=#[2];} } } int itb, err; if (i1 < 1) {itb = 100;} else {itb = ipar[1];} if (i1 < 2) {err = 10;} else {err = ipar[2];} if (itb == 0) { dbprint(prot,"// ** iteration bound reached with error > error bound!"); return(a); } int i,j,k; ideal p=G; matrix J=jacob(G); list h; poly hh; int fertig=1; int n=nvars(basering); for (i = 1; i <= n; i++) { for (j = n; j >= n-i+1; j--) { p[i] = subst(p[i],var(j),a[j]); for (k = n; k >= n-i+1; k--) { J[i,k] = subst(J[i,k],var(j),a[j]); } } if (J[i,n-i+1] == 0) { ERROR("// ideal not radical"); return(); } // solve linear equations hh = -p[i]; for (j = n; j >= n-i+2; j--) { hh = hh - J[i,j]*h[j]; } h[n-i+1] = number(hh/J[i,n-i+1]); } for (i = 1; i <= n; i++) { if ( absValue(h[i]) > (1/10)^err) { fertig = 0; break; } } if ( not fertig ) { if (prot > 0) { "// error:"; print(absValue(h[i])); "// iterations to be performed: "+string(itb); } for (i = 1; i <= n; i++) { a[i] = a[i] + h[i]; } ipar = itb-1,err; return(triMNewton(G,a,ipar)); } else { return(a); } } example { "EXAMPLE:"; echo = 2; ring r = (real,30),(z,y,x),(lp); ideal i = x^2-1,y^2+x4-3,z2-y4+x-1; ideal a = 2,3,4; intvec e = 20,10; ideal l = triMNewton(i,a,e); l; } ///////////////////////////////////////////////////////////////////////////////