42 #ifndef O2SCL_GSL_RKF45_H
43 #define O2SCL_GSL_RKF45_H
49 #include <o2scl/err_hnd.h>
50 #include <o2scl/ode_funct.h>
51 #include <o2scl/ode_step.h>
53 #ifndef DOXYGEN_NO_O2NS
63 template<
class vec_y_t=boost::numeric::ublas::vector<
double>,
64 class vec_dydx_t=vec_y_t,
class vec_yerr_t=vec_y_t,
65 class func_t=ode_funct>
67 vec_dydx_t,vec_yerr_t,func_t> {
73 vec_dydx_t k2, k3, k4, k5, k6;
83 double ah[5], b3[2], b4[3], b5[4], b6[5];
84 double c1, c3, c4, c5, c6;
103 b4[1]=-7200.0/2197.0;
107 b5[1]=-32832.0/4104.0;
108 b5[2]=29440.0/4104.0;
111 b6[0]=-6080.0/20520.0;
112 b6[1]=41040.0/20520.0;
113 b6[2]=-28352.0/20520.0;
114 b6[3]=9295.0/20520.0;
115 b6[4]=-5643.0/20520.0;
117 c1=902880.0/7618050.0;
118 c3=3953664.0/7618050.0;
119 c4=3855735.0/7618050.0;
120 c5=-1371249.0/7618050.0;
121 c6=277020.0/7618050.0;
127 ec[4]=-2197.0/75240.0;
152 virtual int step(
double x,
double h,
size_t n, vec_y_t &y, vec_dydx_t &dydx,
153 vec_y_t &yout, vec_yerr_t &yerr, vec_dydx_t &dydx_out,
172 ytmp[i]=y[i]+ah[0]*h*dydx[i];
179 ytmp[i]=y[i]+h*(b3[0]*dydx[i]+b3[1]*k2[i]);
186 ytmp[i]=y[i]+h*(b4[0]*dydx[i]+b4[1]*k2[i]+b4[2]*k3[i]);
193 ytmp[i]=y[i]+h*(b5[0]*dydx[i]+b5[1]*k2[i]+b5[2]*k3[i]+
201 ytmp[i]=y[i]+h*(b6[0]*dydx[i]+b6[1]*k2[i]+b6[2]*k3[i]+
202 b6[3]*k4[i]+b6[4]*k5[i]);
209 yout[i]=y[i]+h*(c1*dydx[i]+c3*k3[i]+c4*k4[i]+c5*k5[i]+c6*k6[i]);
218 yerr[i]=h*(ec[1]*dydx[i]+ec[3]*k3[i]+ec[4]*k4[i]+ec[5]*k5[i]+
229 #ifndef DOXYGEN_NO_O2NS