42 #ifndef O2SCL_GSL_RKCK_H
43 #define O2SCL_GSL_RKCK_H
49 #include <o2scl/err_hnd.h>
50 #include <o2scl/ode_funct.h>
51 #include <o2scl/ode_step.h>
53 #ifndef DOXYGEN_NO_O2NS
65 template<
class vec_y_t=boost::numeric::ublas::vector<
double>,
66 class vec_dydx_t=vec_y_t,
class vec_yerr_t=vec_y_t,
68 public ode_step<vec_y_t,vec_dydx_t,vec_yerr_t,func_t> {
75 vec_dydx_t k2, k3, k4, k5, k6;
84 double ah[5], b3[2], b4[3], b5[4], b6[5], ec[7];
85 double b21, c1, c3, c4, c6;
111 b6[0]=1631.0/55296.0;
114 b6[3]=44275.0/110592.0;
118 ec[1]=37.0/378.0-2825.0/27648.0;
120 ec[3]=250.0/621.0-18575.0/48384.0;
121 ec[4]=125.0/594.0-13525.0/55296.0;
122 ec[5]=-277.0/14336.0;
123 ec[6]=512.0/1771.0-1.0/4.0;
154 virtual int step(
double x,
double h,
size_t n, vec_y_t &y,
155 vec_dydx_t &dydx, vec_y_t &yout, vec_yerr_t &yerr,
156 vec_dydx_t &dydx_out, func_t &derivs) {
173 ytmp[i]=y[i]+b21*h*dydx[i];
179 ytmp[i]=y[i]+h*(b3[0]*dydx[i]+b3[1]*k2[i]);
185 ytmp[i]=y[i]+h*(b4[0]*dydx[i]+b4[1]*k2[i]+b4[2]*k3[i]);
191 ytmp[i]=y[i]+h*(b5[0]*dydx[i]+b5[1]*k2[i]+b5[2]*k3[i]+
198 ytmp[i]=y[i]+h*(b6[0]*dydx[i]+b6[1]*k2[i]+b6[2]*k3[i]+
199 b6[3]*k4[i]+b6[4]*k5[i]);
205 yout[i]=y[i]+h*(c1*dydx[i]+c3*k3[i]+c4*k4[i]+c6*k6[i]);
214 yerr[i]=h*(ec[1]*dydx[i]+ec[3]*k3[i]+ec[4]*k4[i]+ec[5]*k5[i]+
225 #ifndef DOXYGEN_NO_O2NS