continewz: include/continewz/homotopy

00001 
00002 /******************************************************************************
00003 * MODULE     : homotopy_euler.hpp
00004 * DESCRIPTION: Homotopies using first order Euler method and Newton correction
00005 * COPYRIGHT  : (C) 2009  Joris van der Hoeven
00006 *******************************************************************************
00007 * This software falls under the GNU general public license and comes WITHOUT
00008 * ANY WARRANTY WHATSOEVER. See the file $TEXMACS_PATH/LICENSE for more details.
00009 * If you don't have this file, write to the Free Software Foundation, Inc.,
00010 * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
00011 ******************************************************************************/
00012 
00013 #ifndef __MMX_HOMOTOPY_EULER_HPP
00014 #define __MMX_HOMOTOPY_EULER_HPP
00015 #include <continewz/homotopy.hpp>
00016 
00017 namespace mmx {
00018 #define Function slp_tangent<C>
00019 
00020 struct homotopy_euler {};
00021 template<typename C> struct homotopy_variant_helper {
00022   typedef homotopy_euler HV; };
00023 
00024 /******************************************************************************
00025 * Subroutines
00026 ******************************************************************************/
00027 
00028 template<typename C> matrix<C>
00029 remove_column (const matrix<C>& m, nat k) {
00030   matrix<C> r (C (0), rows (m), cols (m)-1);
00031   for (nat i=0; i<rows(r); i++)
00032     for (nat j=0; j<cols(r); j++)
00033       if (j<k) r(i,j)= m(i,j);
00034       else r(i,j)= m(i,j+1);
00035   return r;
00036 }
00037 
00038 template<typename C> vector<C>
00039 remove_entry (const vector<C>& v, nat k) {
00040   return append (range (v, 0, k), range (v, k+1, N(v)));
00041 }
00042 
00043 template<typename C> vector<Abs_type(C) >
00044 mul_bound (const matrix<C>& m, const vector<C>& v) {
00045   typedef Abs_type(C) R;
00046   ASSERT (N(v) == cols (m), "dimensions don't match");
00047   vector<R> r= fill<R> (0, rows (m));
00048   for (nat i=0; i<N(r); i++) {
00049     R sum= 0;
00050     for (nat j=0; j<cols(m); j++)
00051       sum += abs (m (i, j)) * abs (v[j]);
00052     r[i]= sum;
00053   }
00054   return r;
00055 }
00056 
00057 template<typename C> xnat
00058 coherence (const vector<C>& v1, const vector<C>& v2, xnat p) {
00059   nat n= N (v1);
00060   for (nat i=0; i<n; i++) {
00061     if (v1[i] == 0) return 0;
00062     if (v1[i] == v2[i]) continue;
00063     xint e= exponent (v1[i] / (v1[i] - v2[i]));
00064     p= min (p, (xnat) max (e, (xint) 0));
00065   }
00066   return p;
00067 }
00068 
00069 template<typename F, typename C> xnat
00070 coherence (const F& f, const vector<C>& x, const vector<C>& y,
00071            const matrix<C>& J, nat k, const vector<C>& dd) {
00072   xnat p= precision (promote (1, CF(x))) >> 1;
00073   while (true) {
00074     C eta= decexp2 (promote (1, CF(x)), (xint) p);
00075     vector<C> dx =
00076       eta * append (range (dd, 0, k),
00077                     cons (promote (0, eta), range (dd, k, N(dd))));
00078     vector<C> x2 = x + dx;
00079     //mmerr << "x  = " << x << "\n";
00080     //mmerr << "x2 = " << x2 << "\n";
00081     //mmerr << "dx = " << dx << "\n";
00082     vector<C> df1= J * dx;
00083     vector<C> y2 = eval (f, x2);
00084     //mmerr << "y  = " << y << "\n";
00085     //mmerr << "y2 = " << y2 << "\n";
00086     vector<C> df2= y2 - y;
00087     //mmerr << "dy1= " << df1 << "\n";
00088     //mmerr << "dy2= " << df2 << "\n";
00089     xnat q= coherence (df1, df2, p);
00090     //mmerr << "Coherence= " << q << "\n";
00091     if (q >= (p >> 1) || p <= 10) { p= q; break; }
00092     p= p >> 1;
00093   }
00094   //mmerr << "Coherence= " << p << "\n";
00095   return p;
00096 }
00097 
00098 /******************************************************************************
00099 * Numeric homotopies
00100 ******************************************************************************/
00101 
00102 template<typename C>
00103 struct homotopy_stepper_rep<C,homotopy_euler>:
00104   public homotopy_stepper_rep<C,homotopy_abstract>
00105 {
00106 public:
00107   typedef Abs_type(C) R;
00108 
00109 public:
00110   vector<C> init;
00111   nat       k;
00112   C         t0;
00113   C         t1;
00114   xnat      coh;
00115 
00116 public:
00117   inline homotopy_stepper_rep (const Function& f):
00118     homotopy_stepper_rep<C,homotopy_abstract> (f) {}
00119 
00120   vector<C>
00121   newton_step (const vector<C>& x, const C& dt, R& sc) {
00122     //mmerr << "-----------------------\n";
00123     //mmerr << "x = " << x << "\n";
00124     //mmerr << "dt= " << dt << "\n";
00125     vector<C> y = eval (this->f, x);
00126     //mmerr << "y = " << y << "\n";
00127     matrix<C> J = jacobian (this->f, x);
00128     //mmerr << "J = " << J << "\n";
00129     matrix<C> I = invert (remove_column (J, k));
00130     //mmerr << "I = " << I << "\n";
00131     vector<C> dy= column (J, k) * dt;
00132     //mmerr << "dy= " << dy << "\n";
00133     vector<C> dd= -(I * (y + dy));
00134     vector<C> dx= append (range (dd, 0, k), cons (dt, range (dd, k, N(dd))));
00135     vector<C> x2= x + dx;
00136     //mmerr << "x2= " << x2 << "\n";
00137     sc= l2_norm (mul_bound (J, dx));
00138     //mmerr << "s1= " << l2_norm (eval (this->f, x2)) << "\n";
00139     //mmerr << "s2= " << sc << "\n";
00140     if (dt != 0) coh= min (coh, coherence (this->f, x, y, J, k, dd));
00141     return x2;
00142   }
00143 
00144   vector<C>
00145   propose_next (const vector<C>& x, const C& dt) {
00146     R s2;
00147     vector<C> x2= newton_step (x, dt, s2);
00148     vector<C> y2= eval (this->f, x2);
00149     if (promote (4, s2) * l2_norm (y2) > s2) return x;
00150     R s3;
00151     vector<C> x3= newton_step (x2, 0, s3);
00152     vector<C> y3= eval (this->f, x3);
00153     s3 += s2 / promote (100, s2);
00154     if (promote (10, s3) * l2_norm (y3) > s3) return x;
00155     R s4;
00156     vector<C> x4= newton_step (x3, 0, s4);
00157     vector<C> y4= eval (this->f, x4);
00158     s4 += s3 / promote (10000, s3);
00159     if (promote (40, s4) * l2_norm (y4) > s4) return x;
00160     return x4;
00161   }
00162 
00163   vector<C>
00164   newton_iterate (const vector<C>& x) {
00165     vector<C> x1= x;
00166     vector<C> y1= eval (this->f, x1);
00167     while (true) {
00168       R sc;
00169       vector<C> x2= newton_step (x1, 0, sc);
00170       vector<C> y2= eval (this->f, x2);
00171       if (promote (11, sc) * l2_norm (y2) >=
00172           promote (10, sc) * l2_norm (y1)) return x1;
00173       x1= x2;
00174       y1= y2;
00175     }
00176   }
00177 
00178   vector<C>
00179   hom (const vector<C>& init2, nat k2, const C& t1b) {
00180     //mmerr << "From x = " << init2 << " to x[" << k2 << "] = " << t1b
00181     //<< " [euler stepper]" << "\n";
00182     ASSERT (N(init2) == N(this->f) + 1, "homotopy stepper not initialized");
00183     ASSERT (k2 < N(init2), "coordinate out of range");
00184     init= init2;
00185     k= k2;
00186     t0= init[k];
00187     t1= t1b;
00188     coh= precision (promote (1, t0));
00189 
00190     int credit= 1000;
00191     vector<C> x= init;
00192     //C dt= (t1 - t0) / credit; // Give method time to operate at right scale
00193     C dt= t1 - t0;
00194     while (abs (t1 - x[k]) > abs (sqrt (Accuracy (C)))) {
00195       credit--;
00196       if (credit < 0) mmerr << "Failed at " << x << "\n";
00197       ASSERT (credit >= 0, "continuation failed");
00198       vector<C> x2= propose_next (x, dt);
00199       if (x2 == x) dt= decexp2 (dt, 1);
00200       else {
00201         x = x2;
00202         dt= dt * sqrt (promote (2, t0));
00203         if (abs (t1 - x[k]) < abs (dt)) dt= t1 - x[k];
00204       }
00205       if (coh <= 10) break;
00206     }
00207     //mmerr << "Steps= " << 1000-credit << "\n";
00208     //mmerr << "Coherence= " << coh << "\n";
00209     return newton_iterate (x);
00210   }
00211 };
00212 
00213 #undef Function
00214 } // namespace mmx
00215 #endif // __MMX_HOMOTOPY_EULER_HPP
include/continewz/homotopy_euler.hpp
