42 #ifndef O2SCL_GSL_RK8PD_H
43 #define O2SCL_GSL_RK8PD_H
49 #include <o2scl/ode_step.h>
50 #include <o2scl/err_hnd.h>
52 #ifndef DOXYGEN_NO_O2NS
64 template<
class vec_y_t=boost::numeric::ublas::vector<
double>,
65 class vec_dydx_t=vec_y_t,
class vec_yerr_t=vec_y_t,
66 class func_t=ode_funct>
68 vec_dydx_t,vec_yerr_t,func_t> {
72 #ifdef O2SCL_NEVER_DEFINED
78 vec_dydx_t k2, k3, k4, k5, k6, k7;
79 vec_dydx_t k8, k9, k10, k11, k12, k13;
89 double Abar[13], A[12], ah[10], b21, b3[2], b4[3], b5[4], b6[5];
90 double b7[6], b8[7], b9[8], b10[9], b11[10], b12[11], b13[12];
98 Abar[0]=14005451.0/335480064.0;
103 Abar[5]=-59238493.0/1068277825.0;
104 Abar[6]=181606767.0/758867731.0;
105 Abar[7]=561292985.0/797845732.0;
106 Abar[8]=-1041891430.0/1371343529.0;
107 Abar[9]=760417239.0/1151165299.0;
108 Abar[10]=118820643.0/751138087.0;
109 Abar[11]=-528747749.0/2220607170.0;
112 A[0]=13451932.0/455176623.0;
117 A[5]=-808719846.0/976000145.0;
118 A[6]=1757004468.0/5645159321.0;
119 A[7]=656045339.0/265891186.0;
120 A[8]=-3867574721.0/1518517206.0;
121 A[9]=465885868.0/322736535.0;
122 A[10]=53011238.0/667516719.0;
132 ah[7]=5490023248.0/9719169821.0;
134 ah[9]=1201146811.0/1299019798.0;
156 b7[0]=29443841.0/614563906.0;
159 b7[3]=77736538.0/692538347.0;
160 b7[4]=-28693883.0/1125000000.0;
161 b7[5]=23124283.0/1800000000.0;
163 b8[0]=16016141.0/946692911.0;
166 b8[3]=61564180.0/158732637.0;
167 b8[4]=22789713.0/633445777.0;
168 b8[5]=545815736.0/2771057229.0;
169 b8[6]=-180193667.0/1043307555.0;
171 b9[0]=39632708.0/573591083.0;
174 b9[3]=-433636366.0/683701615.0;
175 b9[4]=-421739975.0/2616292301.0;
176 b9[5]=100302831.0/723423059.0;
177 b9[6]=790204164.0/839813087.0;
178 b9[7]=800635310.0/3783071287.0;
180 b10[0]=246121993.0/1340847787.0;
183 b10[3]=-37695042795.0/15268766246.0;
184 b10[4]=-309121744.0/1061227803.0;
185 b10[5]=-12992083.0/490766935.0;
186 b10[6]=6005943493.0/2108947869.0;
187 b10[7]=393006217.0/1396673457.0;
188 b10[8]=123872331.0/1001029789.0;
190 b11[0]=-1028468189.0/846180014.0;
193 b11[3]=8478235783.0/508512852.0;
194 b11[4]=1311729495.0/1432422823.0;
195 b11[5]=-10304129995.0/1701304382.0;
196 b11[6]=-48777925059.0/3047939560.0;
197 b11[7]=15336726248.0/1032824649.0;
198 b11[8]=-45442868181.0/3398467696.0;
199 b11[9]=3065993473.0/597172653.0;
201 b12[0]=185892177.0/718116043.0;
204 b12[3]=-3185094517.0/667107341.0;
205 b12[4]=-477755414.0/1098053517.0;
206 b12[5]=-703635378.0/230739211.0;
207 b12[6]=5731566787.0/1027545527.0;
208 b12[7]=5232866602.0/850066563.0;
209 b12[8]=-4093664535.0/808688257.0;
210 b12[9]=3962137247.0/1805957418.0;
211 b12[10]=65686358.0/487910083.0;
213 b13[0]=403863854.0/491063109.0;
216 b13[3]=-5068492393.0/434740067.0;
217 b13[4]=-411421997.0/543043805.0;
218 b13[5]=652783627.0/914296604.0;
219 b13[6]=11173962825.0/925320556.0;
220 b13[7]=-13158990841.0/6184727034.0;
221 b13[8]=3936647629.0/1978049680.0;
222 b13[9]=-160528059.0/685178525.0;
223 b13[10]=248638103.0/1413531060.0;
248 virtual int step(
double x,
double h,
size_t n, vec_y_t &y,
249 vec_dydx_t &dydx, vec_y_t &yout, vec_yerr_t &yerr,
250 vec_dydx_t &dydx_out, func_t &derivs) {
274 ytmp[i]=y[i]+b21*h*dydx[i];
280 ytmp[i]=y[i]+h*(b3[0]*dydx[i]+b3[1]*k2[i]);
286 ytmp[i]=y[i]+h*(b4[0]*dydx[i]+b4[2]*k3[i]);
292 ytmp[i]=y[i]+h*(b5[0]*dydx[i]+b5[2]*k3[i]+b5[3]*k4[i]);
298 ytmp[i]=y[i]+h*(b6[0]*dydx[i]+b6[3]*k4[i]+b6[4]*k5[i]);
304 ytmp[i]=y[i]+h*(b7[0]*dydx[i]+b7[3]*k4[i]+b7[4]*k5[i]+b7[5]*k6[i]);
310 ytmp[i]=y[i]+h*(b8[0]*dydx[i]+b8[3]*k4[i]+b8[4]*k5[i]+b8[5]*k6[i]+
317 ytmp[i]=y[i]+h*(b9[0]*dydx[i]+b9[3]*k4[i]+b9[4]*k5[i]+b9[5]*k6[i]+
318 b9[6]*k7[i]+b9[7]*k8[i]);
324 ytmp[i]=y[i]+h*(b10[0]*dydx[i]+b10[3]*k4[i]+b10[4]*k5[i]+
325 b10[5]*k6[i]+b10[6]*k7[i]+b10[7]*k8[i]+
332 ytmp[i]=y[i]+h*(b11[0]*dydx[i]+b11[3]*k4[i]+b11[4]*k5[i]+
333 b11[5]*k6[i]+b11[6]*k7[i]+b11[7]*k8[i]+
334 b11[8]*k9[i]+b11[9]*k10[i]);
340 ytmp[i]=y[i]+h*(b12[0]*dydx[i]+b12[3]*k4[i]+b12[4]*k5[i]+
341 b12[5]*k6[i]+b12[6]*k7[i]+b12[7]*k8[i]+
342 b12[8]*k9[i]+b12[9]*k10[i]+b12[10]*k11[i]);
348 ytmp[i]=y[i]+h*(b13[0]*dydx[i]+b13[3]*k4[i]+b13[4]*k5[i]+
349 b13[5]*k6[i]+b13[6]*k7[i]+b13[7]*k8[i]+
350 b13[8]*k9[i]+b13[9]*k10[i]+b13[10]*k11[i]+
359 double ksum8=Abar[0]*dydx[i]+Abar[5]*k6[i]+Abar[6]*k7[i]+
360 Abar[7]*k8[i]+Abar[8]*k9[i]+Abar[9]*k10[i]+
361 Abar[10]*k11[i]+Abar[11]*k12[i]+Abar[12]*k13[i];
363 yout[i]=y[i]+h*ksum8;
373 double ksum8=Abar[0]*dydx[i]+Abar[5]*k6[i]+Abar[6]*k7[i]+
374 Abar[7]*k8[i]+Abar[8]*k9[i]+Abar[9]*k10[i]+
375 Abar[10]*k11[i]+Abar[11]*k12[i]+Abar[12]*k13[i];
376 double ksum7=A[0]*dydx[i]+A[5]*k6[i]+A[6]*k7[i]+A[7]*k8[i]+
377 A[8]*k9[i]+A[9]*k10[i]+A[10]*k11[i]+A[11]*k12[i];
379 yerr[i]=h*(ksum7-ksum8);
389 #ifndef DOXYGEN_NO_O2NS