43 #ifndef O2SCL_GSL_BSIMP_H
44 #define O2SCL_GSL_BSIMP_H
50 #include <gsl/gsl_math.h>
52 #include <o2scl/err_hnd.h>
53 #include <o2scl/ode_step.h>
54 #include <o2scl/ode_jac_funct.h>
56 #include <o2scl/permutation.h>
58 #ifndef DOXYGEN_NO_O2NS
109 typedef boost::numeric::ublas::unbounded_array<double> ubarr;
113 #ifndef DOXYGEN_INTERAL
157 vec_t y_extrap_sequence;
182 w[i]=(u>0.0) ? u : 1.0;
194 const double safety_f=0.25;
195 const double small_eps=safety_f*eps2;
202 {2,6,10,14,22,34,50,70};
206 a_work[0]=bd_sequence[0]+1.0;
209 a_work[k+1]=a_work[k]+bd_sequence[k+1];
215 const double tmp1=a_work[k+1]-a_work[i+1];
216 const double tmp2=(a_work[i+1]-a_work[0]+1.0)*(2*k+1);
217 alpha[k][i]=pow(small_eps,tmp1/tmp2);
221 a_work[0]+=dimension;
224 a_work[k+1]=a_work[k]+bd_sequence[k+1];
228 if (a_work[k+2]>a_work[k+1]*alpha[k][k+1]) {
247 const unsigned int i_step,
const double x_i,
248 const vec_t &y_i, vec_t &y_0, vec_t &y_0_err,
257 for (j=0;j<
dim;j++) {
265 for (k=0;k<i_step;k++) {
266 double deltaloc=1.0/(x[i_step-k-1]-x_i);
267 const double f1=deltaloc*x_i;
268 const double f2=deltaloc*x[i_step-k-1];
270 for (j=0;j<
dim;j++) {
271 const double q_kj=dloc(k,j);
272 dloc(k,j)=y_0_err[j];
273 deltaloc=work[j]-q_kj;
274 y_0_err[j]=f1*deltaloc;
280 for (j=0;j<
dim;j++) {
281 dloc(i_step,j)=y_0_err[j];
296 const unsigned int n_step,
const ubarr &y,
297 const ubarr &yp_loc,
const vec_t &dfdt_loc,
298 const mat_t &dfdy_loc, vec_t &y_out) {
300 const double h=h_total/n_step;
310 const double max_sum=100.0*
dim;
317 for (i=0;i<
dim;i++) {
318 for (j=0;j<
dim;j++) {
319 a_mat(i,j)=-h*dfdy_loc(i,j);
335 for (i=0;i<
dim;i++) {
336 y_temp[i]=h*(yp_loc[i]+h*dfdt_loc[i]);
342 for (i=0;i<
dim;i++) {
343 const double di=delta_temp[i];
346 sum+=fabs(di)/weight[i];
354 status=(*funcp)(t,
dim,y_temp,y_out);
359 for (n_inter=1;n_inter<n_step;n_inter++) {
361 for (i=0;i<
dim;i++) {
362 rhs_temp[i]=h*y_out[i]-delta[i];
368 for (i=0;i<
dim;i++) {
369 delta[i]+=2.0*delta_temp[i];
371 sum+=fabs(delta[i])/weight[i];
379 status=(*funcp)(t,
dim,y_temp,y_out);
388 for (i=0;i<
dim;i++) {
389 rhs_temp[i]=h*y_out[i]-delta[i];
395 for (i=0;i<
dim;i++) {
396 y_out[i]=y_temp[i]+delta_temp[i];
397 sum+=fabs(delta_temp[i])/weight[i];
421 for(
size_t i=0;i<n;i++) yp[i]=0.0;
425 y_extrap_save.resize(n);
426 extrap_work.resize(n);
431 y_extrap_sequence.resize(n);
434 delta_temp.resize(n);
438 double sqrt_dbl_eps=sqrt(std::numeric_limits<double>::epsilon());
443 k_choice=k_choice_loc;
444 k_current=k_choice_loc;
445 order=2*k_choice_loc;
455 extrap_work.resize(0);
456 y_extrap_save.resize(0);
474 virtual ~ode_bsimp_gsl() {
483 for(
size_t i=0;i<
dim;i++) yp[i]=0.0;
500 virtual int step(
double x,
double h,
size_t n, vec_t &y, vec_t &dydx,
501 vec_t &yout, vec_t &yerr, vec_t &dydx_out,
502 func_t &derivs, jac_func_t &jac) {
512 static const int bd_sequence[
sequence_count]={2,6,10,14,22,34,50,70};
518 if (h+t_local==t_local) {
522 O2SCL_ERR(
"Stepsize underflow in ode_bsimp_gsl::step().",
535 ret=jac(t_local,
dim,y,dfdy,dfdt);
545 for (k=0;k<=k_current;k++) {
547 const unsigned int N=bd_sequence[k];
548 const double r=(h/N);
549 const double x_k=r*r;
554 int status=
step_local(t_local,h,N,y_extrap_save,yp,
555 dfdt,dfdy,y_extrap_sequence);
558 std::cout <<
"ode_bsimp_gsl: " << k <<
"/" << k_current <<
" "
559 <<
"status=" << status << std::endl;
566 for (i=0;i<
dim;i++) {
582 ret=derivs(t_local+h,
dim,y,dydx_out);
596 #ifndef DOXYGEN_NO_O2NS