[380a17b] | 1 | ////////////////////////////////////////////////////////////////////////////// |
---|
[3686937] | 2 | version="version ntsolve.lib 4.0.0.0 Jun_2013 "; // $Id$ |
---|
[fd3fb7] | 3 | category="Symbolic-numerical solving"; |
---|
[5138c3e] | 4 | info=" |
---|
[8942a5] | 5 | LIBRARY: ntsolve.lib Real Newton Solving of Polynomial Systems |
---|
[a5f564] | 6 | AUTHORS: Wilfred Pohl, email: pohl@mathematik.uni-kl.de |
---|
[c948324] | 7 | Dietmar Hillebrand |
---|
[5138c3e] | 8 | |
---|
| 9 | PROCEDURES: |
---|
[9173792] | 10 | nt_solve(G,ini,[..]); find one real root of 0-dimensional ideal G |
---|
| 11 | triMNewton(G,a,[..]); find one real root for 0-dim triangular system G |
---|
[5138c3e] | 12 | "; |
---|
[a5f564] | 13 | |
---|
[6d37e8] | 14 | LIB "general.lib"; |
---|
[5138c3e] | 15 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 16 | |
---|
[a5f564] | 17 | proc nt_solve (ideal gls, ideal ini, list #) |
---|
[9173792] | 18 | "USAGE: nt_solve(gls,ini[,ipar]); gls,ini= ideals, ipar=list/intvec,@* |
---|
[731e67e] | 19 | gls: contains the equations, for which a solution will be computed@* |
---|
[9173792] | 20 | ini: ideal of initial values (approximate solutions to start with),@* |
---|
[0bc582c] | 21 | ipar: control integers (default: ipar = [100, 10]) |
---|
[9173792] | 22 | @format |
---|
| 23 | ipar[1]: max. number of iterations |
---|
[0bc582c] | 24 | ipar[2]: accuracy (we have the l_2-norm ||.||): accepts solution @code{sol} |
---|
[50cbdc] | 25 | if ||gls(sol)|| < eps0*(0.1^ipar[2]) |
---|
[9173792] | 26 | where eps0 = ||gls(ini)|| is the initial error |
---|
| 27 | @end format |
---|
[a5f564] | 28 | ASSUME: gls is a zerodimensional ideal with nvars(basering) = size(gls) (>1) |
---|
| 29 | RETURN: ideal, coordinates of one solution (if found), 0 else |
---|
[9173792] | 30 | NOTE: if printlevel >0: displays comments (default =0) |
---|
[5138c3e] | 31 | EXAMPLE: example nt_solve; shows an example |
---|
| 32 | " |
---|
| 33 | { |
---|
| 34 | def rn = basering; |
---|
| 35 | int di = size(gls); |
---|
| 36 | if (nvars(basering) != di){ |
---|
[a5f564] | 37 | ERROR("// wrong number of equations");} |
---|
[5138c3e] | 38 | if (size(ini) != di){ |
---|
[a5f564] | 39 | ERROR("// wrong number of initial values");} |
---|
[782fcd8] | 40 | int prec = system("getPrecDigits"); // precision |
---|
[5138c3e] | 41 | |
---|
| 42 | int i1,i2,i3; |
---|
[a5f564] | 43 | int itmax, acc; |
---|
| 44 | intvec ipar; |
---|
| 45 | if ( size(#)>0 ){ |
---|
| 46 | i1=1; |
---|
| 47 | if (typeof(#[1])=="intvec") {ipar=#[1];} |
---|
| 48 | if (typeof(#[1])=="int") {ipar[1]=#[1];} |
---|
| 49 | if ( size(#)>1 ){ |
---|
| 50 | i1=2; |
---|
| 51 | if (typeof(#[2])=="int") {ipar[2]=#[2];} |
---|
| 52 | } |
---|
| 53 | } |
---|
| 54 | |
---|
| 55 | int prot = printlevel-voice+2; // prot=printlevel (default:prot=0) |
---|
[5138c3e] | 56 | if (i1 < 1){itmax = 100;}else{itmax = ipar[1];} |
---|
[1f9a84] | 57 | if (i1 < 2){acc = prec div 2;}else{acc = ipar[2];} |
---|
[782fcd8] | 58 | if ((acc <= 0)||(acc > prec-1)){acc = prec-1;} |
---|
[5138c3e] | 59 | |
---|
| 60 | int dpl = di+1; |
---|
[b9b906] | 61 | string out; |
---|
[5138c3e] | 62 | out = "ring rnewton=(real,prec),("+varstr(basering)+"),(c,dp);"; |
---|
| 63 | execute(out); |
---|
[a5f564] | 64 | |
---|
[5138c3e] | 65 | ideal gls1=imap(rn,gls); |
---|
| 66 | module nt,sub; |
---|
| 67 | sub = transpose(jacob(gls1)); |
---|
| 68 | for (i1=di;i1>0;i1--){ |
---|
| 69 | if(sub[i1]==0){break;}} |
---|
| 70 | if (i1>0){ |
---|
| 71 | setring rn; kill rnewton; |
---|
[a5f564] | 72 | ERROR("// one var not in equation");} |
---|
| 73 | |
---|
[782fcd8] | 74 | list direction; |
---|
| 75 | ideal ini1; |
---|
[5138c3e] | 76 | ini1 = imap(rn,ini); |
---|
| 77 | number dum,y1,y2,y3,genau; |
---|
| 78 | genau = 0.1; |
---|
| 79 | dum = genau; |
---|
| 80 | genau = genau^acc; |
---|
[a5f564] | 81 | |
---|
[5138c3e] | 82 | for (i1=di;i1>0;i1--){ |
---|
| 83 | sub[i1]=sub[i1]+gls1[i1]*gen(dpl);} |
---|
| 84 | nt = sub; |
---|
| 85 | for (i1=di;i1>0;i1--){ |
---|
| 86 | nt = subst(nt,var(i1),ini1[i1]);} |
---|
[a5f564] | 87 | |
---|
[5138c3e] | 88 | // now we have in sub the general structure |
---|
[782fcd8] | 89 | // and in nt the structure with subst. vars |
---|
[5138c3e] | 90 | |
---|
| 91 | // compute initial error |
---|
| 92 | y1 = ml2norm(nt,genau); |
---|
[b9b906] | 93 | dbprint(prot,"// initial error = "+string(y1)); |
---|
[5138c3e] | 94 | y2 = genau*y1; |
---|
| 95 | |
---|
| 96 | // begin of iteration |
---|
| 97 | for(i3=1;i3<=itmax;i3++){ |
---|
[b9b906] | 98 | dbprint(prot,"// iteration: "+string(i3)); |
---|
[5138c3e] | 99 | |
---|
| 100 | // find newton direction |
---|
[782fcd8] | 101 | direction=bareiss(nt,1,-1); |
---|
[5138c3e] | 102 | |
---|
| 103 | // find dumping |
---|
[782fcd8] | 104 | dum = linesearch(gls1,ini1,direction[1],y1,dum,genau); |
---|
| 105 | if (i3%5 == 0) |
---|
| 106 | { |
---|
| 107 | if (dum <= 0.000001) |
---|
| 108 | { |
---|
| 109 | dum = 1.0; |
---|
| 110 | } |
---|
| 111 | } |
---|
[a5f564] | 112 | dbprint(prot,"// dumping = "+string(dum)); |
---|
[5138c3e] | 113 | |
---|
| 114 | // new value |
---|
| 115 | for(i1=di;i1>0;i1--){ |
---|
[782fcd8] | 116 | ini1[i1]=ini1[i1]-dum*direction[1][i1];} |
---|
[5138c3e] | 117 | nt = sub; |
---|
| 118 | for (i1=di;i1>0;i1--){ |
---|
| 119 | nt = subst(nt,var(i1),ini1[i1]);} |
---|
| 120 | y1 = ml2norm(nt,genau); |
---|
[a5f564] | 121 | dbprint(prot,"// error = "+string(y1)); |
---|
[5138c3e] | 122 | if(y1<y2){break;} // we are ready |
---|
| 123 | } |
---|
| 124 | |
---|
| 125 | if (y1>y2){ |
---|
[a5f564] | 126 | "// ** WARNING: iteration bound reached with error > error bound!";} |
---|
[782fcd8] | 127 | setring rn; |
---|
[5138c3e] | 128 | ini = imap(rnewton,ini1); |
---|
| 129 | kill rnewton; |
---|
| 130 | return(ini); |
---|
| 131 | } |
---|
| 132 | example |
---|
| 133 | { |
---|
| 134 | "EXAMPLE:";echo=2; |
---|
[782fcd8] | 135 | ring rsq = (real,40),(x,y,z,w),lp; |
---|
[5138c3e] | 136 | ideal gls = x2+y2+z2-10, y2+z3+w-8, xy+yz+xz+w5 - 1,w3+y; |
---|
[782fcd8] | 137 | ideal ini = 3.1,2.9,1.1,0.5; |
---|
[a5f564] | 138 | intvec ipar = 200,0; |
---|
[782fcd8] | 139 | ideal sol = nt_solve(gls,ini,ipar); |
---|
[5138c3e] | 140 | sol; |
---|
| 141 | } |
---|
| 142 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 143 | |
---|
| 144 | static proc sqrt (number wr, number wa, number wg) |
---|
| 145 | { |
---|
[782fcd8] | 146 | number es,we; |
---|
[5138c3e] | 147 | number wb=wa; |
---|
| 148 | number wf=wb*wb-wr; |
---|
[782fcd8] | 149 | if(wf>0){ |
---|
| 150 | es=wf;} |
---|
| 151 | else{ |
---|
| 152 | es=-wf;} |
---|
| 153 | we=wg*es; |
---|
| 154 | while (es>we) |
---|
[5138c3e] | 155 | { |
---|
| 156 | wf=wf/(wb+wb); |
---|
| 157 | wb=wb-wf; |
---|
| 158 | wf=wb*wb-wr; |
---|
[782fcd8] | 159 | if(wf>0){ |
---|
| 160 | es=wf;} |
---|
| 161 | else{ |
---|
| 162 | es=-wf;} |
---|
[5138c3e] | 163 | } |
---|
| 164 | return(wb); |
---|
| 165 | } |
---|
| 166 | |
---|
| 167 | static proc il2norm (ideal H, number wg) |
---|
| 168 | { |
---|
| 169 | number wa,wb; |
---|
| 170 | int wi,dpl; |
---|
| 171 | wa = leadcoef(H[1]); |
---|
| 172 | wa = wa*wa; |
---|
| 173 | for(wi=size(H);wi>1;wi--) |
---|
| 174 | { |
---|
| 175 | wb=leadcoef(H[wi]); |
---|
| 176 | wa=wa+wb*wb; |
---|
| 177 | } |
---|
| 178 | return(sqrt(wa,wa,wg)); |
---|
| 179 | } |
---|
| 180 | |
---|
| 181 | static proc ml2norm (module H, number wg) |
---|
| 182 | { |
---|
| 183 | number wa,wb; |
---|
| 184 | int wi,dpl; |
---|
| 185 | dpl = size(H)+1; |
---|
| 186 | wa = leadcoef(H[1][dpl]); |
---|
| 187 | wa = wa*wa; |
---|
| 188 | for(wi=size(H);wi>1;wi--) |
---|
| 189 | { |
---|
| 190 | wb=leadcoef(H[wi][dpl]); |
---|
| 191 | wa=wa+wb*wb; |
---|
| 192 | } |
---|
| 193 | return(sqrt(wa,wa,wg)); |
---|
| 194 | } |
---|
| 195 | |
---|
[ab907a] | 196 | static |
---|
[782fcd8] | 197 | proc linesearch(ideal nl, ideal aa, ideal bb, |
---|
[5138c3e] | 198 | number z1, number tt, number gg) |
---|
| 199 | { |
---|
| 200 | int ii,d; |
---|
[782fcd8] | 201 | ideal cc,jn; |
---|
| 202 | number ss,z2,z3,mm; |
---|
| 203 | |
---|
| 204 | mm=0.000001; |
---|
| 205 | ss=tt; |
---|
[5138c3e] | 206 | d=size(nl); |
---|
| 207 | cc=aa; |
---|
[782fcd8] | 208 | for(ii=d;ii>0;ii--){cc[ii]=cc[ii]-ss*bb[ii];} |
---|
[5138c3e] | 209 | jn=nl; |
---|
| 210 | for(ii=d;ii>0;ii--){jn=subst(jn,var(ii),cc[ii]);} |
---|
| 211 | z2=il2norm(jn,gg); |
---|
| 212 | z3=-1; |
---|
| 213 | while(z2>=z1) |
---|
| 214 | { |
---|
[782fcd8] | 215 | ss=0.5*ss; |
---|
| 216 | if(ss<mm){return (mm);} |
---|
[5138c3e] | 217 | cc=aa; |
---|
| 218 | for(ii=d;ii>0;ii--) |
---|
| 219 | { |
---|
[782fcd8] | 220 | cc[ii]=cc[ii]-ss*bb[ii]; |
---|
[5138c3e] | 221 | } |
---|
| 222 | jn=nl; |
---|
| 223 | for(ii=d;ii>0;ii--){jn=subst(jn,var(ii),cc[ii]);} |
---|
| 224 | z3=z2; |
---|
| 225 | z2=il2norm(jn,gg); |
---|
| 226 | } |
---|
| 227 | if(z3<0) |
---|
| 228 | { |
---|
| 229 | while(z3<z2) |
---|
| 230 | { |
---|
[782fcd8] | 231 | ss=ss+ss; |
---|
[5138c3e] | 232 | cc=aa; |
---|
| 233 | for(ii=d;ii>0;ii--) |
---|
| 234 | { |
---|
[782fcd8] | 235 | cc[ii]=cc[ii]-ss*bb[ii]; |
---|
[5138c3e] | 236 | } |
---|
| 237 | jn=nl; |
---|
| 238 | for(ii=d;ii>0;ii--){jn=subst(jn,var(ii),cc[ii]);} |
---|
| 239 | if(z3>0){z2=z3;} |
---|
| 240 | z3=il2norm(jn,gg); |
---|
| 241 | } |
---|
| 242 | } |
---|
| 243 | z2=z2-z1; |
---|
| 244 | z3=z3-z1; |
---|
[782fcd8] | 245 | ss=0.25*ss*(z3-4*z2)/(z3-2*z2); |
---|
| 246 | if(ss>1.0){return (1.0);} |
---|
| 247 | if(ss<mm){return (mm);} |
---|
| 248 | return(ss); |
---|
[5138c3e] | 249 | } |
---|
[c948324] | 250 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 251 | // |
---|
| 252 | // Multivariate Newton for triangular systems |
---|
[a5f564] | 253 | // algorithms for solving algebraic system of dimension zero |
---|
| 254 | // written by Dietmar Hillebrand |
---|
[c948324] | 255 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 256 | |
---|
[a5f564] | 257 | proc triMNewton (ideal G, ideal a, list #) |
---|
[9173792] | 258 | "USAGE: triMNewton(G,a[,ipar]); G,a= ideals, ipar=list/intvec |
---|
| 259 | ASSUME: G: g1,..,gn, a triangular system of n equations in n vars, i.e. |
---|
| 260 | gi=gi(var(n-i+1),..,var(n)),@* |
---|
| 261 | a: ideal of numbers, coordinates of an approximation of a common |
---|
| 262 | zero of G to start with (with a[i] to be substituted in var(i)),@* |
---|
[0bc582c] | 263 | ipar: control integer vector (default: ipar = [100, 10]) |
---|
[9173792] | 264 | @format |
---|
| 265 | ipar[1]: max. number of iterations |
---|
| 266 | ipar[2]: accuracy (we have as norm |.| absolute value ): |
---|
[0bc582c] | 267 | accepts solution @code{sol} if |G(sol)| < |G(a)|*(0.1^ipar[2]). |
---|
[9173792] | 268 | @end format |
---|
[a5f564] | 269 | RETURN: an ideal, coordinates of a better approximation of a zero of G |
---|
[c948324] | 270 | EXAMPLE: example triMNewton; shows an example |
---|
| 271 | " |
---|
| 272 | { |
---|
[a5f564] | 273 | int prot = printlevel; |
---|
| 274 | int i1,i2,i3; |
---|
| 275 | intvec ipar; |
---|
| 276 | if ( size(#)>0 ){ |
---|
| 277 | i1=1; |
---|
| 278 | if (typeof(#[1])=="intvec") {ipar=#[1];} |
---|
| 279 | if (typeof(#[1])=="int") {ipar[1]=#[1];} |
---|
| 280 | if ( size(#)>1 ){ |
---|
| 281 | i1=2; |
---|
| 282 | if (typeof(#[2])=="int") {ipar[2]=#[2];} |
---|
| 283 | } |
---|
| 284 | } |
---|
| 285 | int itb, err; |
---|
| 286 | if (i1 < 1) {itb = 100;} else {itb = ipar[1];} |
---|
| 287 | if (i1 < 2) {err = 10;} else {err = ipar[2];} |
---|
[b9b906] | 288 | |
---|
[c948324] | 289 | if (itb == 0) |
---|
| 290 | { |
---|
[a5f564] | 291 | dbprint(prot,"// ** iteration bound reached with error > error bound!"); |
---|
| 292 | return(a); |
---|
[c948324] | 293 | } |
---|
| 294 | |
---|
| 295 | int i,j,k; |
---|
| 296 | ideal p=G; |
---|
| 297 | matrix J=jacob(G); |
---|
| 298 | list h; |
---|
| 299 | poly hh; |
---|
| 300 | int fertig=1; |
---|
| 301 | int n=nvars(basering); |
---|
| 302 | |
---|
| 303 | for (i = 1; i <= n; i++) |
---|
| 304 | { |
---|
| 305 | for (j = n; j >= n-i+1; j--) |
---|
| 306 | { |
---|
| 307 | p[i] = subst(p[i],var(j),a[j]); |
---|
| 308 | for (k = n; k >= n-i+1; k--) |
---|
| 309 | { |
---|
| 310 | J[i,k] = subst(J[i,k],var(j),a[j]); |
---|
| 311 | } |
---|
| 312 | } |
---|
| 313 | if (J[i,n-i+1] == 0) |
---|
| 314 | { |
---|
[a5f564] | 315 | ERROR("// ideal not radical"); |
---|
[c948324] | 316 | return(); |
---|
| 317 | } |
---|
| 318 | |
---|
| 319 | // solve linear equations |
---|
| 320 | hh = -p[i]; |
---|
| 321 | for (j = n; j >= n-i+2; j--) |
---|
| 322 | { |
---|
| 323 | hh = hh - J[i,j]*h[j]; |
---|
| 324 | } |
---|
| 325 | h[n-i+1] = number(hh/J[i,n-i+1]); |
---|
| 326 | } |
---|
| 327 | |
---|
| 328 | for (i = 1; i <= n; i++) |
---|
| 329 | { |
---|
[6d37e8] | 330 | if ( absValue(h[i]) > (1/10)^err) |
---|
[c948324] | 331 | { |
---|
| 332 | fertig = 0; |
---|
| 333 | break; |
---|
| 334 | } |
---|
| 335 | } |
---|
| 336 | if ( not fertig ) |
---|
| 337 | { |
---|
[a5f564] | 338 | if (prot > 0) |
---|
| 339 | { |
---|
[6d37e8] | 340 | "// error:"; print(absValue(h[i])); |
---|
[a5f564] | 341 | "// iterations to be performed: "+string(itb); |
---|
| 342 | } |
---|
[c948324] | 343 | for (i = 1; i <= n; i++) |
---|
| 344 | { |
---|
| 345 | a[i] = a[i] + h[i]; |
---|
| 346 | } |
---|
[a5f564] | 347 | ipar = itb-1,err; |
---|
| 348 | return(triMNewton(G,a,ipar)); |
---|
[c948324] | 349 | } |
---|
| 350 | else |
---|
| 351 | { |
---|
| 352 | return(a); |
---|
| 353 | } |
---|
| 354 | } |
---|
| 355 | example |
---|
| 356 | { "EXAMPLE:"; echo = 2; |
---|
[a5f564] | 357 | ring r = (real,30),(z,y,x),(lp); |
---|
| 358 | ideal i = x^2-1,y^2+x4-3,z2-y4+x-1; |
---|
| 359 | ideal a = 2,3,4; |
---|
| 360 | intvec e = 20,10; |
---|
| 361 | ideal l = triMNewton(i,a,e); |
---|
[c948324] | 362 | l; |
---|
| 363 | } |
---|
[731e67e] | 364 | /////////////////////////////////////////////////////////////////////////////// |
---|