# source:git/Singular/LIB/ntsolve.lib@c948324

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