source:git/Singular/LIB/ntsolve.lib@6d37e8

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