source: git/Singular/LIB/ntsolve.lib @ 1f9a84

spielwiese
Last change on this file since 1f9a84 was 1f9a84, checked in by Hans Schoenemann <hannes@…>, 13 years ago
more int division from the manual git-svn-id: file:///usr/local/Singular/svn/trunk@14203 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 8.8 KB
Line 
1///////////////////////////////////////////////////////////////////////////////
2version="$Id$";
3category="Symbolic-numerical solving";
4info="
5LIBRARY:  ntsolve.lib     Real Newton Solving of Polynomial Systems
6AUTHORS:  Wilfred Pohl, email: pohl@mathematik.uni-kl.de
7          Dietmar Hillebrand
8
9PROCEDURES:
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
12";
13
14LIB "general.lib";
15///////////////////////////////////////////////////////////////////////////////
16
17proc nt_solve (ideal gls, ideal ini, list #)
18"USAGE:   nt_solve(gls,ini[,ipar]); gls,ini= ideals, ipar=list/intvec,@*
19  gls: contains the equations, for which a solution will be computed@*
20  ini: ideal of initial values (approximate solutions to start with),@*
21  ipar: control integers (default: ipar = [100, 10])
22  @format
23 ipar[1]: max. number of iterations
24 ipar[2]: accuracy (we have the l_2-norm ||.||): accepts solution @code{sol}
25          if ||gls(sol)|| < eps0*(0.1^ipar[2])
26          where eps0 = ||gls(ini)|| is the initial error
27  @end format
28ASSUME:  gls is a zerodimensional ideal with nvars(basering) = size(gls) (>1)
29RETURN:  ideal, coordinates of one solution (if found), 0 else
30NOTE:    if printlevel >0: displays comments (default =0)
31EXAMPLE: example nt_solve; shows an example
32"
33{
34    def rn = basering;
35    int di = size(gls);
36    if (nvars(basering) != di){
37      ERROR("// wrong number of equations");}
38    if (size(ini) != di){
39      ERROR("// wrong number of initial values");}
40    int prec = system("getPrecDigits"); // precision
41
42    int i1,i2,i3;
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)
56    if (i1 < 1){itmax = 100;}else{itmax = ipar[1];}
57    if (i1 < 2){acc = prec div 2;}else{acc = ipar[2];}
58    if ((acc <= 0)||(acc > prec-1)){acc = prec-1;}
59
60    int dpl = di+1;
61    string out;
62    out = "ring rnewton=(real,prec),("+varstr(basering)+"),(c,dp);";
63    execute(out);
64
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;
72      ERROR("// one var not in equation");}
73
74    list direction;
75    ideal ini1;
76    ini1 = imap(rn,ini);
77    number dum,y1,y2,y3,genau;
78    genau = 0.1;
79    dum = genau;
80    genau = genau^acc;
81
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]);}
87
88    // now we have in sub the general structure
89    // and in nt the structure with subst. vars
90
91    // compute initial error
92    y1 = ml2norm(nt,genau);
93    dbprint(prot,"// initial error = "+string(y1));
94    y2 = genau*y1;
95
96  // begin of iteration
97  for(i3=1;i3<=itmax;i3++){
98     dbprint(prot,"// iteration: "+string(i3));
99
100    // find newton direction
101    direction=bareiss(nt,1,-1);
102
103    // find dumping
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    }
112    dbprint(prot,"// dumping = "+string(dum));
113
114    // new value
115    for(i1=di;i1>0;i1--){
116      ini1[i1]=ini1[i1]-dum*direction[1][i1];}
117    nt = sub;
118    for (i1=di;i1>0;i1--){
119      nt = subst(nt,var(i1),ini1[i1]);}
120    y1 = ml2norm(nt,genau);
121    dbprint(prot,"// error = "+string(y1));
122    if(y1<y2){break;} // we are ready
123  }
124
125    if (y1>y2){
126      "// ** WARNING: iteration bound reached with error > error bound!";}
127    setring rn;
128    ini = imap(rnewton,ini1);
129    kill rnewton;
130    return(ini);
131}
132example
133{
134    "EXAMPLE:";echo=2;
135    ring rsq = (real,40),(x,y,z,w),lp;
136    ideal gls =  x2+y2+z2-10, y2+z3+w-8, xy+yz+xz+w5 - 1,w3+y;
137    ideal ini = 3.1,2.9,1.1,0.5;
138    intvec ipar = 200,0;
139    ideal sol = nt_solve(gls,ini,ipar);
140    sol;
141}
142///////////////////////////////////////////////////////////////////////////////
143
144static proc sqrt (number wr, number wa, number wg)
145{
146  number es,we;
147  number wb=wa;
148  number wf=wb*wb-wr;
149  if(wf>0){
150    es=wf;}
151  else{
152    es=-wf;}
153  we=wg*es;
154  while (es>we)
155  {
156    wf=wf/(wb+wb);
157    wb=wb-wf;
158    wf=wb*wb-wr;
159    if(wf>0){
160      es=wf;}
161    else{
162      es=-wf;}
163  }
164  return(wb);
165}
166
167static 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
181static 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
196static
197proc linesearch(ideal nl, ideal aa, ideal bb,
198number z1, number tt, number gg)
199{
200  int ii,d;
201  ideal cc,jn;
202  number ss,z2,z3,mm;
203
204  mm=0.000001;
205  ss=tt;
206  d=size(nl);
207  cc=aa;
208  for(ii=d;ii>0;ii--){cc[ii]=cc[ii]-ss*bb[ii];}
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  {
215    ss=0.5*ss;
216    if(ss<mm){return (mm);}
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    z3=z2;
225    z2=il2norm(jn,gg);
226  }
227  if(z3<0)
228  {
229    while(z3<z2)
230    {
231      ss=ss+ss;
232      cc=aa;
233      for(ii=d;ii>0;ii--)
234      {
235        cc[ii]=cc[ii]-ss*bb[ii];
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;
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);
249}
250///////////////////////////////////////////////////////////////////////////////
251//
252//   Multivariate Newton for triangular systems
253//   algorithms for solving algebraic system of dimension zero
254//   written by Dietmar Hillebrand
255///////////////////////////////////////////////////////////////////////////////
256
257proc triMNewton (ideal G, ideal a, list #)
258"USAGE:  triMNewton(G,a[,ipar]); G,a= ideals, ipar=list/intvec
259ASSUME:  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)),@*
263  ipar: control integer vector (default: ipar = [100, 10])
264  @format
265  ipar[1]: max. number of iterations
266  ipar[2]: accuracy (we have as norm |.| absolute value ):
267           accepts solution @code{sol} if |G(sol)| < |G(a)|*(0.1^ipar[2]).
268  @end format
269RETURN:  an ideal, coordinates of a better approximation of a zero of G
270EXAMPLE: example triMNewton; shows an example
271"
272{
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];}
288
289    if (itb == 0)
290    {
291       dbprint(prot,"// ** iteration bound reached with error > error bound!");
292       return(a);
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        {
315            ERROR("// ideal not radical");
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    {
330        if ( absValue(h[i]) > (1/10)^err)
331        {
332            fertig = 0;
333            break;
334        }
335    }
336    if ( not fertig )
337    {
338        if (prot > 0)
339        {
340           "// error:"; print(absValue(h[i]));
341           "// iterations to be performed: "+string(itb);
342        }
343        for (i = 1; i <= n; i++)
344        {
345            a[i] = a[i] + h[i];
346        }
347        ipar = itb-1,err;
348        return(triMNewton(G,a,ipar));
349    }
350    else
351    {
352        return(a);
353    }
354}
355example
356{ "EXAMPLE:"; echo = 2;
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);
362   l;
363}
364///////////////////////////////////////////////////////////////////////////////
Note: See TracBrowser for help on using the repository browser.