Changeset 2054 for branches/newmole/source

Timestamp:

05/13/08 16:02:39 (8 months ago)

Author:

rjrw

Message:

Wrap comment in mole_solve

Cubic backtrack in newton_step

Location:

branches/newmole/source

Files:

• mole_solve.cpp (modified) (1 diff)
• newton_step.cpp (modified) (6 diffs)

branches/newmole/source/mole_solve.cpp

@@ -19 +19 @@
 #include "newton_step.h"
 #include "householder.h"
-/* Nick Abel between July and October of 2003 assisted Dr. Ferland in improving the heavy element 
- * molecular network in Cloudy. Before this routine would predict negative abundances if 
- * the fraction of carbon in the form of molecules came close to 100%. A reorganizing of 
- * the reaction network detected several bugs.  Treatment of "coupled reactions",
- * in which both densities in the reaction rate were being predicted by Cloudy, were also 
- * added.  Due to these improvements, Cloudy can now perform calculations
- * where 100% of the carbon is in the form of CO without predicting negative abundances
+/* Nick Abel between July and October of 2003 assisted Dr. Ferland in
+ * improving the heavy element molecular network in Cloudy. Before
+ * this routine would predict negative abundances if the fraction of
+ * carbon in the form of molecules came close to 100%. A reorganizing
+ * of the reaction network detected several bugs.  Treatment of
+ * "coupled reactions", in which both densities in the reaction rate
+ * were being predicted by Cloudy, were also added.  Due to these
+ * improvements, Cloudy can now perform calculations where 100% of the
+ * carbon is in the form of CO without predicting negative abundances
  *
  * Additional changes were made in November of 2003 so that our reaction 

branches/newmole/source/newton_step.cpp

@@ -32 +32 @@
                 etmp,
                 emax,
-                fstep,
-                flast,
+                f1,
+                f2,
+                error2,
                 error1,
                 error0,
@@ -52 +53 @@
         }
 
-        error0 = 1;
-        error1 = grad = 0.;
-        fstep = 1.;                     
+        error0 = error2 = grad = f1 = f2 = 0.;
 
         for (loop=0;loop<20;loop++) 
         {
-                flast = fstep;
-                
                 jacobn(b2vec,&ervals[0],&amat[0],loop==0);
 
@@ -70 +67 @@
                 }
 
-                error0 = emax = 0.;
+                error1 = emax = 0.;
                 iworst = -1;
                 for( i=0; i < n; i++ )
@@ -82 +79 @@
                                 iworst = i;
                         }
-                        error0 += etmp;
+                        error1 += etmp;
                 }
 
@@ -92 +89 @@
                         }
                         merror = solve_system(amat,b1vec,n);
-                        error1 = error0;
-                        grad = -2*error1;
+                        error0 = error1;
+                        grad = -2*error0;
+                        f1 = 1.0;
                 } else {
-                        if (error0 < (1-2e-4*fstep)*error1 || error0 < 1e-20)
+                        if (error1 < (1-2e-4*f1)*error0 || error1 < 1e-20)
                         {
                                 break;
+                        } 
+                        if (loop == 1) 
+                        {
+                                f2 = 1.;
+                                f1 = -0.5*grad/SDIV(error1-error0-grad);
                         }
-                        fstep = -0.5*grad*fstep*fstep/SDIV(error0-error1-grad*fstep);
-                        if (fstep > 0.9*flast || fstep < 0.)
-                                fstep = 0.5*flast;
-                        if (fstep < 0.1*flast)
-                                fstep = 0.1*flast;
+                        else
+                        {
+                                // Update using cubic fit, ref Dennis & Schnabel 1996
+                                double a1 = error1-error0-grad*f1;
+                                double b1 = error2-error0-grad*f2;
+                                double a = (a1/(f1*f1)-b1/(f2*f2));
+                                double b = (-f2*a1/(f1*f1)+f1*b1/(f2*f2));
+                                double ft = f1;
+                                f1 = (-b+sqrt(b*b-3*a*grad*(f1-f2)))/(3*a);
+                                f2 = ft;
+                        }
+                        error2 = error1;
+                        if (f1 > 0.5*f2 || f1 < 0.)
+                                f1 = 0.5*f2;
+                        else if (f1 < 0.1*f2)
+                                f1 = 0.1*f2;
                 }
 
@@ -113 +127 @@
                 for( i=0; i < n; i++ )
                 {
-                        b2vec[i] = b0vec[i]-b1vec[i]*fstep;                             
+                        b2vec[i] = b0vec[i]-b1vec[i]*f1;                                
                 }
         }
 
-        *error = (realnum) MIN2(error0,1e30);
+        *error = (realnum) MIN2(error1,1e30);
 
         return;