Changeset c948324 in git

Timestamp:

Dec 15, 2000, 1:08:46 PM (23 years ago)

Author:

Hans Schönemann <hannes@…>

Branches:

(u'spielwiese', '4a9821a93ffdc22a6696668bd4f6b8c9de3e6c5f')

Children:

edef3080503be9dafb0b4d33562cb634ef9d50a8

Parents:

116c0963c62cef54596fdf23c8b1827f59377a05

Message:

*hannes: moved solver.lib to ntsolve.lib


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

Location:

Singular/LIB

Files:

• ntsolve.lib (modified) (2 diffs)
• solver.lib (deleted)

Singular/LIB/ntsolve.lib

@@ -1 +1 @@
 ///////////////////////////////////////////////////////////////////////////////
-version="$Id: ntsolve.lib,v 1.6 2000-07-06 13:32:25 pohl Exp $";
+version="$Id: ntsolve.lib,v 1.7 2000-12-15 12:08:46 Singular Exp $";
 info="
 LIBRARY: ntsolve.lib ONE REAL SOLUTION OF POLYNOMIAL SYSTEMS (NEWTON ITERATION)
-AUTHOR:  Wilfred Pohl,   email: pohl@mathematik.uni-kl.de
+AUTHORS:  Wilfred Pohl,   email: pohl@mathematik.uni-kl.de
+          Dietmar Hillebrand
 
 PROCEDURES:
-nt_solve(i,..);        find one real root of 0-dimensional ideal
+nt_solve(G,..);        find one real root of 0-dimensional ideal G
+triMNewton(G,..);      find one real root for 0-dim. triangular system G
 ";
 ///////////////////////////////////////////////////////////////////////////////
@@ -233 +235 @@
 ///////////////////////////////////////////////////////////////////////////////
 
-// End: ***
+///////////////////////////////////////////////////////////////////////////////
+// solver.lib                                                                //
+// algorithms for solving algebraic system of dimension zero                 //
+// written by Dietmar Hillebrand and ...                                     //
+//                                                                           //
+///////////////////////////////////////////////////////////////////////////////
+
+
+
+///////////////////////////////////////////////////////////////////////////////
+//
+//   Multivariate Newton for triangular systems
+//
+//
+///////////////////////////////////////////////////////////////////////////////
+
+proc triMNewton (ideal G, list a, number err, int itb)
+"USAGE:  triMNewtion(G,a,err,itb);
+             ideal G=g1,..,gn, a triangular system
+                 in n vars, i.e gi=gi(var(n-i+1),..,var(n));
+             list of numbers a, an approximation of a common zero of G,
+                 (a[i] to be substituted in var(i));
+             number err, an error bound;
+             int itb, the maximal number of iterations performed.
+RETURN:  an improved approximation for a common zero of G;
+EXAMPLE: example triMNewton; shows an example
+"
+{
+    if (itb == 0)
+    {
+        dbprint("iteration bound performed!");
+        return(a);
+    }
+
+    int i,j,k;
+    ideal p=G;
+    matrix J=jacob(G);
+    list h;
+    poly hh;
+    int fertig=1;
+    int n=nvars(basering);
+
+    for (i = 1; i <= n; i++)
+    {
+        for (j = n; j >= n-i+1; j--)
+        {
+            p[i] = subst(p[i],var(j),a[j]);
+            for (k = n; k >= n-i+1; k--)
+            {
+                J[i,k] = subst(J[i,k],var(j),a[j]);
+            }
+        }
+        if (J[i,n-i+1] == 0)
+        {
+            print("Error: ideal not radical!");
+            return();
+        }
+
+        // solve linear equations
+        hh = -p[i];
+        for (j = n; j >= n-i+2; j--)
+        {
+            hh = hh - J[i,j]*h[j];
+        }
+        h[n-i+1] = number(hh/J[i,n-i+1]);
+    }
+
+    for (i = 1; i <= n; i++)
+    {
+        if ( abs(h[i]) > err)
+        {
+            fertig = 0;
+            break;
+        }
+    }
+    if ( not fertig )
+    {
+        for (i = 1; i <= n; i++)
+        {
+            a[i] = a[i] + h[i];
+        }
+        return(triMNewton(G,a,err,itb-1));
+    }
+    else
+    {
+        return(a);
+    }
+}
+example
+{ "EXAMPLE:"; echo = 2;
+   ring r=real,(z,y,x),(lp);
+   ideal i=x^2-1,y^2+x4-3,z2-y4+x-1;
+   list a=2,3,4;
+   number e=1.0e-10;
+
+   list l = triMNewton(i,a,e,20);
+   l;
+}
+
+
+////////////////////////////////////////////////////////////////////////////////
+
+static proc abs( number r)
+{
+   if (r >= 0)
+   {
+       return(r);
+   }
+   else
+   {
+       return(-r);
+   }
+}
+