continewz: include/continewz/ode_taylor.hpp Source File

00001 
00002 /******************************************************************************
00003 * MODULE     : ode_taylor.hpp
00004 * DESCRIPTION: Taylor solutions to dynamical systems
00005 * COPYRIGHT  : (C) 2011  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_ODE_TAYLOR_HPP
00014 #define __MMX_ODE_TAYLOR_HPP
00015 #include <analyziz/newton_polygon.hpp>
00016 #include <continewz/ode_series.hpp>
00017 
00018 namespace mmx {
00019 #define Homotopy_variant(C) typename homotopy_variant_helper<C>::HV
00020 #define TMPL template<typename C>
00021 #define R Abs_type(C)
00022 #define BR Default_radius_type(C)
00023 #define Polynomial polynomial<C>
00024 #define Series series<C>
00025 #define Taylor taylor<C>
00026 #define Function_C slp_tangent<C>
00027 #define Function_Taylor slp_tangent<Taylor >
00028 
00029 /******************************************************************************
00030 * Order and precision
00031 ******************************************************************************/
00032 
00033 extern int mmx_ode_order;
00034 extern int mmx_ode_precision;
00035 
00036 template<typename C> nat
00037 ode_order (const C& c) {
00038   if (mmx_ode_order != 0)
00039     return min ((nat) precision (c), (nat) mmx_ode_order);
00040   return precision (c);
00041 }
00042 
00043 template<typename C> xnat
00044 ode_precision (const C& c) {
00045   if (mmx_ode_precision != 0)
00046     return min ((xnat) precision (c), (xnat) mmx_ode_precision);
00047   return min ((xnat) precision (c), max ((xnat) ode_order (c), (xnat) 24));
00048 }
00049 
00050 /******************************************************************************
00051 * Extra subroutines for vectors
00052 ******************************************************************************/
00053 
00054 TMPL R
00055 max_newton_scale (const vector<Polynomial >& v, nat order) {
00056   if (N(v) == 0 || N(v[0]) <= 1)
00057     return abs (promote (1, get_format1 (CF (v))));
00058   nat i, j;
00059   double* expo_src= mmx_new<double> (order);
00060   for (nat k=0; k<order; k++) {
00061     double expo= (double) exponent (v[0][k]);
00062     for (nat i=1; i<N(v); i++)
00063       expo= max (expo, (double) exponent (v[i][k]));
00064     expo_src[k]= expo;
00065   }
00066   newton_edge<double> (i, j, expo_src, order, order >> 1);
00067   double slope= log (2.0) * (expo_src[j] - expo_src[i]) / ((double) (j-i));
00068   mmx_delete<double> (expo_src, order);
00069   return exp (as<R> (slope));
00070 }
00071 
00072 TMPL R
00073 heuristic_domain (const vector<Polynomial >& v, nat order, nat prec) {
00074   R sc1= max_newton_scale (v, order);
00075   R sc2= exp2 (promote ((int) prec, sc1) / promote ((int) order, sc1));
00076   //mmerr << "sc1 = " << sc1 << "\n";
00077   //mmerr << "sc2 = " << sc2 << "\n";
00078   return promote (1, sc1) / (sc1 * sc2);
00079 }
00080 
00081 /******************************************************************************
00082 * Integration of dynamical systems
00083 ******************************************************************************/
00084 
00085 TMPL vector<Taylor >
00086 solve_ode_certify (const Function_C& fun, const vector<ball<C> >& c,
00087                    const vector<Polynomial>& psol,
00088                    nat order, xnat prec, const R& dom2) {
00089   nat n= N(c);
00090   R dom= dom2;
00091   Function_Taylor tfun= as<Function_Taylor > (fun);
00092   nat iters= (nat) log2 ((double) (order + 16));
00093   for (nat turn=0; turn<order; turn++) {
00094     //mmerr << "dom = " << dom << "\n";
00095     vector<Taylor > tsol= fill<Taylor > (n);
00096     for (nat i=0; i<n; i++) {
00097       tsol[i]= Taylor (psol[i], order, promote (0, dom), dom);
00098       R err= as<R> (radius (c[i])) + decexp2 (sum_abs_up (tsol[i]), prec);
00099       tsol[i]= blur (tsol[i], err);
00100     }
00101     for (nat iter=0; iter<iters; iter++) {
00102       //mmerr << "---\n";
00103       //mmerr << "tsol= " << tsol << "\n";
00104       vector<Taylor> next= integrate_init (eval (tfun, tsol), center (c));
00105       //mmerr << "next= " << next << "\n";
00106       if (included (next, tsol)) return next;
00107       if (included (tsol, next)) break;
00108       for (nat i=0; i<n; i++) {
00109         R delta= sup (radius (tsol[i]), radius (next[i])) - radius (tsol[i]);
00110         R rad  = radius (tsol[i]) + incexp2 (delta, 1);
00111         tsol[i]= Taylor (psol[i], order, rad, dom);
00112       }
00113     }
00114     dom= decexp2 (dom, 1);
00115   }
00116   ERROR ("computation of Taylor bound failed");
00117 }
00118 
00119 TMPL vector<Taylor >
00120 solve_ode_taylor (const Function_C& fun, const vector<ball<C> >& c,
00121                   nat order, xnat prec, const R& dom) {
00122   Function_Taylor tfun= as<Function_Taylor > (fun);
00123   vector<Series > ssol= solve_ode_series (fun, center (c));
00124   vector<Polynomial > psol= range (ssol, 0, order);
00125   return solve_ode_certify (fun, c, psol, order, prec, dom);
00126 }
00127 
00128 TMPL vector<Taylor >
00129 solve_ode_taylor (const Function_C& fun, const vector<ball<C> >& c,
00130                   nat order, xnat prec) {
00131   Function_Taylor tfun= as<Function_Taylor > (fun);
00132   vector<Series > ssol= solve_ode_series (fun, center (c));
00133   vector<Polynomial > psol= range (ssol, 0, order);
00134   R dom= heuristic_domain (psol, order, prec);
00135   return solve_ode_certify (fun, c, psol, order, prec, dom);
00136 }
00137 
00138 TMPL vector<generic>
00139 solve_ode_taylor (const vector<multivariate<sparse_polynomial<rational> > >& f,
00140                   const vector<multivariate_coordinate<> >& vars,
00141                   const vector<C>& c) {
00142   ASSERT (N(c)>0, "at least one equation expected");
00143   nat order= ode_order (c[0]);
00144   xnat prec= ode_precision (c[0]);
00145   Function_C fun= as_slp_tangent (f, vars, CF (c));
00146   vector<ball<C> > b= as<vector<ball<C> > > (c);
00147   return as<vector<generic> > (solve_ode_taylor (fun, b, order, prec));
00148 }
00149 
00150 /******************************************************************************
00151 * First variation at solution
00152 ******************************************************************************/
00153 
00154 TMPL matrix<Taylor >
00155 first_variation_certify (const Function_C& fun, const vector<ball<C> >& c,
00156                          const vector<Taylor >& tsol2,
00157                          const matrix<Polynomial >& pvar,
00158                          nat order, xnat prec) {
00159   nat n= N(c);
00160   vector<Taylor > tsol= tsol2;
00161   R dom= domain (tsol[0]);
00162   Function_Taylor tfun= as<Function_Taylor > (fun);
00163   nat iters= (nat) log2 ((double) (order + 16));
00164   matrix<C> id= identity_matrix (n, CF (center (c)));
00165   for (nat turn=0; turn<order; turn++) {
00166     //mmerr << "dom = " << dom << "\n";
00167     matrix<Taylor > tjac= jacobian (tfun, tsol);
00168     matrix<Taylor > tvar (promote (0, tsol[0]), n, n);
00169     for (nat i=0; i<n; i++)
00170       for (nat j=0; j<n; j++) {
00171         tvar (i, j)= Taylor (pvar (i, j), order, promote (0, dom), dom);
00172         R err= decexp2 (sum_abs_up (tvar (i, j)), prec);
00173         tvar (i, j)= blur (tvar (i, j), err);
00174       }
00175     for (nat iter=0; iter<iters; iter++) {
00176       //mmerr << "---\n";
00177       //mmerr << "tvar= " << tvar << "\n";
00178       matrix<Taylor> next= integrate_init (tjac * tvar, id);
00179       //mmerr << "next= " << next << "\n";
00180       if (included (next, tvar)) return next;
00181       if (included (tvar, next)) break;
00182       for (nat i=0; i<n; i++)
00183         for (nat j=0; j<n; j++) {
00184           R delta= sup (radius (tvar (i, j)), radius (next (i, j))) -
00185                    radius (tvar (i, j));
00186           R rad  = radius (tvar (i, j)) + incexp2 (delta, 1);
00187           tvar (i, j)= Taylor (pvar (i, j), order, rad, dom);
00188         }
00189     }
00190     dom= decexp2 (dom, 1);
00191     for (nat i=0; i<n; i++)
00192       tsol[i]= restrict (tsol[i], dom);
00193   }
00194   ERROR ("computation of Taylor bound failed");
00195 }
00196 
00197 TMPL matrix<Taylor >
00198 first_variation_taylor (const Function_C& fun, const vector<ball<C> >& c,
00199                         nat order, xnat prec, const R& dom) {
00200   Function_Taylor tfun= as<Function_Taylor > (fun);
00201   vector<Series > ssol= solve_ode_series (fun, center (c));
00202   matrix<Series > svar= first_variation_series (fun, ssol);
00203   vector<Polynomial > psol= range (ssol, 0, order);
00204   matrix<Polynomial > pvar= range (svar, 0, order);
00205   vector<Taylor > tsol= solve_ode_certify (fun, c, psol, order, prec, dom);
00206   return first_variation_certify (fun, c, tsol, pvar, order, prec);
00207 }
00208 
00209 TMPL matrix<Taylor >
00210 first_variation_taylor (const Function_C& fun, const vector<ball<C> >& c,
00211                         nat order, xnat prec) {
00212   Function_Taylor tfun= as<Function_Taylor > (fun);
00213   vector<Series > ssol= solve_ode_series (fun, center (c));
00214   matrix<Series > svar= first_variation_series (fun, ssol);
00215   vector<Polynomial > psol= range (ssol, 0, order);
00216   matrix<Polynomial > pvar= range (svar, 0, order);
00217   R dom= heuristic_domain (psol, order, prec);
00218   vector<Taylor > tsol= solve_ode_certify (fun, c, psol, order, prec, dom);
00219   return first_variation_certify (fun, c, tsol, pvar, order, prec);
00220 }
00221 
00222 TMPL matrix<generic>
00223 first_variation_taylor (const vector<multivariate<sparse_polynomial<rational> > >& f,
00224                         const vector<multivariate_coordinate<> >& vars,
00225                         const vector<C>& c) {
00226   ASSERT (N(c)>0, "at least one equation expected");
00227   nat order= ode_order (c[0]);
00228   xnat prec= ode_precision (c[0]);
00229   Function_C fun= as_slp_tangent (f, vars, CF (c));
00230   vector<ball<C> > b= as<vector<ball<C> > > (c);
00231   return as<matrix<generic> > (first_variation_taylor (fun, b, order, prec));
00232 }
00233 
00234 /******************************************************************************
00235 * Integrating a dynamical system
00236 ******************************************************************************/
00237 
00238 TMPL vector<ball<C> >
00239 eval (const vector<Taylor >& v, const ball<C>& z) {
00240   vector<ball<C> > r= fill<ball<C> > (promote (0, z), N(v));
00241   for (nat i=0; i<N(v); i++)
00242     r[i]= eval (v[i], z);
00243   return r;
00244 }
00245 
00246 TMPL matrix<ball<C> >
00247 eval (const matrix<Taylor >& m, const ball<C>& z) {
00248   matrix<ball<C> > r (promote (0, z), rows (m), cols (m));
00249   for (nat i=0; i<rows(m); i++)
00250     for (nat j=0; j<cols(m); j++)
00251       r (i, j)= eval (m (i, j), z);
00252   return r;
00253 }
00254 
00255 TMPL void
00256 step_ode_taylor (const Function_C& fun,
00257                  vector<C>& c,
00258                  matrix<C>& var,
00259                  vector<ball<C> >& err,
00260                  ball<C>& target,
00261                  nat order, xnat prec)
00262 {
00263   //mmerr << "c= " << c << "\n";
00264   //mmerr << "var= " << var << "\n";
00265   //mmerr << "err= " << err << "\n";
00266   //mmerr << "target= " << target << "\n";
00267   vector<ball<C> > xc= as<vector<ball<C> > > (c);
00268   vector<ball<C> > ac= xc + as<matrix<ball<C> > > (var) * err;
00269   vector<Taylor > xsol= solve_ode_taylor (fun, xc, order, prec);
00270   matrix<Taylor > avar= first_variation_taylor (fun, ac, order, prec/2);
00271   // FIXME: cheating here, by taking prec/2
00272   //mmerr << "xsol= " << xsol << "\n";
00273   //mmerr << "avar= " << avar << "\n";
00274   R r= min (domain (xsol[0]), domain (avar (0, 0)));
00275   C delta= center (target);
00276   //mmerr << "r= " << r << "\n";
00277   //mmerr << "delta= " << delta << "\n";
00278   if (delta == promote (0, delta)) {
00279     ASSERT (abs_up (target) <= r, "integration failed");
00280     vector<ball<C> > nc  = eval (xsol, target);
00281     matrix<ball<C> > invc= invert (as<matrix<ball<C> > > (var));
00282     c= center (nc);
00283     err= err + invc * (nc - sharpen (nc));
00284     target= promote (0, target);
00285     return;
00286   }
00287   if (abs_up (sharpen (target)) > r) {
00288     R lambda= r / as<R> (abs_up (sharpen (target)));
00289     delta= promote (999, lambda) * lambda * delta / promote (1000, lambda);
00290   }
00291   //mmerr << "delta= " << delta << "\n";
00292   ball<C>          z   = as<ball<C> > (delta);
00293   vector<ball<C> > nc  = eval (xsol, z);
00294   matrix<ball<C> > nvar= eval (avar, z) * as<matrix<ball<C> > > (var);
00295   vector<ball<C> > ec  = nc - sharpen (nc);
00296   matrix<ball<C> > evar= nvar - sharpen (nvar);
00297   matrix<ball<C> > invc= invert (sharpen (nvar));
00298   c= center (nc);
00299   var= center (nvar);
00300   err= err + invc * (ec + evar * err);
00301   if (delta != center (target)) target= target - delta;
00302   else target= ball<C> (promote (0, center (target)), radius (target));
00303 }
00304 
00305 TMPL vector<ball<C> >
00306 integrate_ode_taylor (const Function_C& fun,
00307                       const vector<ball<C> >& bc,
00308                       const ball<C>& z,
00309                       nat order, xnat prec)
00310 {
00311   vector<C> c= center (bc);
00312   matrix<C> var= identity_matrix (N(c), CF(c));
00313   vector<ball<C> > err= bc - sharpen (bc);
00314   ball<C> target= z;
00315   nat counter= 0;
00316   while (!exact_eq (target, promote (0, target))) {
00317     counter++;
00318     //if ((counter % 1000) == 0)
00319     //mmerr << z - target << ", " << c << ", " << var << "\n";
00320     step_ode_taylor (fun, c, var, err, target, order, prec);
00321   }
00322   mmerr << "c= " << c << "\n";
00323   mmerr << "var= " << var << "\n";
00324   return
00325     as<vector<ball<C> > > (c) +
00326     as<matrix<ball<C> > > (var) * err;
00327 }
00328 
00329 TMPL vector<generic>
00330 integrate_ode_taylor (const vector<multivariate<sparse_polynomial<rational> > >& f,
00331                       const vector<multivariate_coordinate<> >& vars,
00332                       const vector<C>& c,
00333                       const C& z) {
00334   ASSERT (N(c)>0, "at least one equation expected");
00335   nat order= ode_order (c[0]);
00336   xnat prec= ode_precision (c[0]);
00337   Function_C fun= as_slp_tangent (f, vars, CF (c));
00338   vector<ball<C> > b= as<vector<ball<C> > > (c);
00339   return as<vector<generic> > (integrate_ode_taylor (fun, b, as<ball<C> > (z), order, prec));
00340 }
00341 
00342 } // namespace mmx
00343 #undef TMPL
00344 #undef R
00345 #undef BR
00346 #undef Polynomial
00347 #undef Series
00348 #undef Taylor
00349 #undef Function_C
00350 #undef Function_Taylor
00351 #endif // __MMX_ODE_TAYLOR_HPP