Changeset a5f564 in git for Singular/LIB/ntsolve.lib


Ignore:
Timestamp:
Dec 16, 2000, 2:14:43 AM (23 years ago)
Author:
Gert-Martin Greuel <greuel@…>
Branches:
(u'spielwiese', 'd0474371d8c5d8068ab70bfb42719c97936b18a6')
Children:
e35f0bd56e7ac915c46b0f6708641f7d5d3e6706
Parents:
4508ce5e48c7f51f11a43264563eaabbb7992df1
Message:
* Greuel: Syntax von nt_solve und triMNewton vereinheitlicht,
* Syntax von nt_solve erweitert, die von triMNewton geaendert


git-svn-id: file:///usr/local/Singular/svn/trunk@4921 2c84dea3-7e68-4137-9b89-c4e89433aadc
File:
1 edited

Legend:

Unmodified
Added
Removed

  • Singular/LIB/ntsolve.lib

    r4508ce5 ra5f564  
    1 ///////////////////////////////////////////////////////////////////////////////
    2 version="$Id: ntsolve.lib,v 1.7 2000-12-15 12:08:46 Singular Exp $";
     1//(GMG, last modified 16.12.00)
     2///////////////////////////////////////////////////////////////////////////////
     3version="$Id: ntsolve.lib,v 1.8 2000-12-16 01:14:43 greuel Exp $";
    34info="
    4 LIBRARY: ntsolve.lib ONE REAL SOLUTION OF POLYNOMIAL SYSTEMS (NEWTON ITERATION)
    5 AUTHORS:  Wilfred Pohl,   email: pohl@mathematik.uni-kl.de
     5LIBRARY:  ntsolve.lib     REAL NEWTON SOLVING OF POLYNOMIAL SYSTEMS
     6AUTHORS:  Wilfred Pohl, email: pohl@mathematik.uni-kl.de
    67          Dietmar Hillebrand
    78
    89PROCEDURES:
    9 nt_solve(G,..);        find one real root of 0-dimensional ideal G
    10 triMNewton(G,..);      find one real root for 0-dim. triangular system G
     10 nt_solve(G,..);        find one real root of 0-dimensional ideal G
     11 triMNewton(G,..);      find one real root for 0-dim triangular system G
    1112";
    12 ///////////////////////////////////////////////////////////////////////////////
    13 
    14 proc 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
    26 ASSUME:  gls is a zerodimensional ideal with
    27          nvars(basering) = size(gls) (> 1)
    28 RETURN:  ideal of one solution (if found)
    29          0 (else)
     13
     14///////////////////////////////////////////////////////////////////////////////
     15
     16proc nt_solve (ideal gls, ideal ini, list #)
     17"USAGE:   nt_solve(gls,ini[,ipar]); gls ideal, ini ideal, ipar list or intvec
     18         gls: contains the equations, for which a solution will be computed
     19         ini: ideal of initial values (approximate solutions to start with)
     20         ipar: control integers (default: ipar = 100,10)
     21         ipar[1] - max. number of iterations
     22         ipar[2] - accuracy (we have the l_2-norm ||.||): accept solution
     23         sol if ||gls(sol)|| < eps0*(0.1^ipar[2]) where eps0 = ||gls(ini)||
     24         is the initial error
     25ASSUME:  gls is a zerodimensional ideal with nvars(basering) = size(gls) (>1)
     26RETURN:  ideal, coordinates of one solution (if found), 0 else
     27NOTE:    printlevel >0: displays comments (default =0)
    3028EXAMPLE: example nt_solve; shows an example
    3129"
    … …  
    3432    int di = size(gls);
    3533    if (nvars(basering) != di){
    36       ERROR("wrong dimension");}
     34      ERROR("// wrong number of equations");}
    3735    if (size(ini) != di){
    38       ERROR("wrong number of initial values");}
     36      ERROR("// wrong number of initial values");}
    3937    int prec = system("getPrecDigits"); // precision
    4038
    4139    int i1,i2,i3;
    42     i1 = size(ipar);
    43     int itmax, acc, prot;
     40    int itmax, acc;
     41    intvec ipar;
     42    if ( size(#)>0 ){
     43       i1=1;
     44       if (typeof(#[1])=="intvec") {ipar=#[1];}
     45       if (typeof(#[1])=="int") {ipar[1]=#[1];}
     46       if ( size(#)>1 ){
     47          i1=2;
     48          if (typeof(#[2])=="int") {ipar[2]=#[2];}
     49       }
     50    }
     51
     52    int prot = printlevel-voice+2;  // prot=printlevel (default:prot=0)
    4453    if (i1 < 1){itmax = 100;}else{itmax = ipar[1];}
    4554    if (i1 < 2){acc = prec/2;}else{acc = ipar[2];}
    46     if (i1 < 3){prot = 0;}else{prot = ipar[3];}
    4755    if ((acc <= 0)||(acc > prec-1)){acc = prec-1;}
    4856
    4957    int dpl = di+1;
    50     string out; // for prot != 0 and more
     58    string out;
    5159    out = "ring rnewton=(real,prec),("+varstr(basering)+"),(c,dp);";
    5260    execute(out);
     61
    5362    ideal gls1=imap(rn,gls);
    5463    module nt,sub;
    … …  
    5867    if (i1>0){
    5968      setring rn; kill rnewton;
    60       ERROR("one var not in equation");}
     69      ERROR("// one var not in equation");}
     70
    6171    list direction;
    6272    ideal ini1;
    … …  
    6676    dum = genau;
    6777    genau = genau^acc;
     78
    6879    for (i1=di;i1>0;i1--){
    6980      sub[i1]=sub[i1]+gls1[i1]*gen(dpl);}
    … …  
    7182    for (i1=di;i1>0;i1--){
    7283      nt = subst(nt,var(i1),ini1[i1]);}
     84
    7385    // now we have in sub the general structure
    7486    // and in nt the structure with subst. vars
    … …  
    7688    // compute initial error
    7789    y1 = ml2norm(nt,genau);
    78     if(prot){out=" initial error = "+string(y1);out;}
     90    dbprint(prot,"// initial error = "+string(y1));
    7991    y2 = genau*y1;
    8092
    8193  // begin of iteration
    8294  for(i3=1;i3<=itmax;i3++){
    83     if(prot){out="  Nr. "+string(i3);out;}
     95     dbprint(prot,"// iteration: "+string(i3));
    8496
    8597    // find newton direction
    … …  
    95107      }
    96108    }
    97     if(prot){out="  dumping = "+string(dum);out;}
     109    dbprint(prot,"// dumping = "+string(dum));
    98110
    99111    // new value
    … …  
    104116      nt = subst(nt,var(i1),ini1[i1]);}
    105117    y1 = ml2norm(nt,genau);
    106     if(prot){out="  error = "+string(y1);out;}
     118    dbprint(prot,"// error = "+string(y1));
    107119    if(y1<y2){break;} // we are ready
    108120  }
    109121
    110122    if (y1>y2){
    111       "WARNING: no convergence";}
     123      "// ** WARNING: iteration bound reached with error > error bound!";}
    112124    setring rn;
    113125    ini = imap(rnewton,ini1);
    … …  
    121133    ideal gls =  x2+y2+z2-10, y2+z3+w-8, xy+yz+xz+w5 - 1,w3+y;
    122134    ideal ini = 3.1,2.9,1.1,0.5;
    123     intvec ipar = 200,0,1;
     135    intvec ipar = 200,0;
    124136    ideal sol = nt_solve(gls,ini,ipar);
    125137    sol;
    … …  
    234246}
    235247///////////////////////////////////////////////////////////////////////////////
    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 ///////////////////////////////////////////////////////////////////////////////
    247248//
    248249//   Multivariate Newton for triangular systems
    249 //
    250 //
    251 ///////////////////////////////////////////////////////////////////////////////
    252 
    253 proc 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.
    261 RETURN:  an improved approximation for a common zero of G;
     250//   algorithms for solving algebraic system of dimension zero
     251//   written by Dietmar Hillebrand
     252///////////////////////////////////////////////////////////////////////////////
     253
     254proc triMNewton (ideal G, ideal a, list #)
     255"USAGE:  triMNewton(G,a[,ipar]); G ideal, a ideal, ipar list or intvec
     256ASSUME:  G:   g1,..,gn, a triangular system of n equations in n vars,
     257              i.e gi=gi(var(n-i+1),..,var(n))
     258         a:   ideal of numbers, coordinates of an approximation of a common
     259              zero of G to start with (with a[i] to be substituted in var(i))
     260         ipar: control integer vector (default: ipar = 100,10)
     261         ipar[1] - max. number of iterations
     262         ipar[2] - accuracy (we have as norm |.| absolute value ):
     263         accept solution  sol if |G(sol)| < |G(a)|*(0.1^ipar[2])
     264         itb: int, the maximal number of iterations performed
     265         err: number, an error bound
     266RETURN:  an ideal, coordinates of a better approximation of a zero of G
    262267EXAMPLE: example triMNewton; shows an example
    263268"
    264269{
     270    int prot = printlevel;
     271    int i1,i2,i3;
     272    intvec ipar;
     273    if ( size(#)>0 ){
     274       i1=1;
     275       if (typeof(#[1])=="intvec") {ipar=#[1];}
     276       if (typeof(#[1])=="int") {ipar[1]=#[1];}
     277       if ( size(#)>1 ){
     278          i1=2;
     279          if (typeof(#[2])=="int") {ipar[2]=#[2];}
     280       }
     281    }
     282    int itb, err;
     283    if (i1 < 1) {itb = 100;} else {itb = ipar[1];}
     284    if (i1 < 2) {err = 10;} else {err = ipar[2];}
     285 
    265286    if (itb == 0)
    266287    {
    267         dbprint("iteration bound performed!");
    268         return(a);
     288       dbprint(prot,"// ** iteration bound reached with error > error bound!");
     289       return(a);
    269290    }
    270291
    … …  
    289310        if (J[i,n-i+1] == 0)
    290311        {
    291             print("Error: ideal not radical!");
     312            ERROR("// ideal not radical");
    292313            return();
    293314        }
    … …  
    304325    for (i = 1; i <= n; i++)
    305326    {
    306         if ( abs(h[i]) > err)
     327        if ( abs(h[i]) > (1/10)^err)
    307328        {
    308329            fertig = 0;
    … …  
    312333    if ( not fertig )
    313334    {
     335        if (prot > 0)
     336        {
     337           "// error:"; print(abs(h[i]));
     338           "// iterations to be performed: "+string(itb);
     339        }
    314340        for (i = 1; i <= n; i++)
    315341        {
    316342            a[i] = a[i] + h[i];
    317343        }
    318         return(triMNewton(G,a,err,itb-1));
     344        ipar = itb-1,err;
     345        return(triMNewton(G,a,ipar));
    319346    }
    320347    else
    … …  
    325352example
    326353{ "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);
     354   ring r = (real,30),(z,y,x),(lp);
     355   ideal i = x^2-1,y^2+x4-3,z2-y4+x-1;
     356   ideal a = 2,3,4;
     357   intvec e = 20,10;
     358   ideal l = triMNewton(i,a,e);
    333359   l;
    334360}
    335 
    336 
    337 ////////////////////////////////////////////////////////////////////////////////
     361///////////////////////////////////////////////////////////////////////////////
    338362
    339363static proc abs( number r)
    … …  
    348372   }
    349373}
    350 
     374///////////////////////////////////////////////////////////////////////////////
     375
Note: See TracChangeset for help on using the changeset viewer.