42 #ifndef O2SCL_ASTEP_GSL_H 43 #define O2SCL_ASTEP_GSL_H 51 #include <gsl/gsl_math.h> 52 #include <gsl/gsl_odeiv.h> 54 #include <o2scl/astep.h> 55 #include <o2scl/vector.h> 56 #include <o2scl/string_conv.h> 58 #ifndef DOXYGEN_NO_O2NS 128 template<
class vec_y_t=boost::numeric::ublas::vector<
double>,
129 class vec_dydx_t=vec_y_t,
class vec_yerr_t=vec_y_t>
177 template<
class svec_t>
int set_scale(
size_t nscal,
const svec_t &scale) {
178 scale_abs.resize(nscal);
179 for(
size_t i=0;i<nscal;i++) {
180 scale_abs[i]=scale[i];
185 virtual int hadjust(
size_t dim,
unsigned int ord,
const vec_y_t &y,
186 vec_yerr_t &yerr, vec_dydx_t &yp,
double &h) {
189 const double h_old=h;
196 if (standard || scale_abs.size()==0) {
197 D0=eps_rel*(a_y*fabs(y[i])+a_dydt*fabs(h_old*yp[i]))+eps_abs;
199 D0=eps_rel*(a_y*fabs(y[i])+a_dydt*fabs(h_old*yp[i]))+
200 eps_abs*scale_abs[i%scale_abs.size()];
202 const double r=fabs(yerr[i]) / fabs(D0);
203 rmax=GSL_MAX_DBL(r, rmax);
210 double r= S / pow(rmax, 1.0/ord);
218 }
else if (rmax < 0.5) {
221 double r=S / pow(rmax, 1.0/(ord+1.0));
266 template<
class vec_y_t=boost::numeric::ublas::vector<
double>,
267 class vec_dydx_t=vec_y_t,
class vec_yerr_t=vec_y_t,
269 public astep_base<vec_y_t,vec_dydx_t,vec_yerr_t,func_t> {
275 #ifndef DOXYGEN_INTERNAL 304 int evolve_apply(
double t0,
double t1,
double &t,
double &h,
size_t nvar,
305 vec_y_t &y, vec_dydx_t &dydx, vec_y_t &yout,
306 vec_yerr_t &yerr, vec_dydx_t &dydx_out,
311 bool final_step=
false;
315 if ((dt < 0.0 && h0 > 0.0) || (dt > 0.0 && h0 < 0.0)) {
316 std::string str=
"Interval direction (t1-t0="+
o2scl::dtos(dt)+
317 ") does not match step direction (h="+
o2scl::dtos(h)+
318 " in astep_gsl::evolve_apply().";
326 bool try_step_again=
true;
327 while (try_step_again) {
328 try_step_again=
false;
330 if ((dt >= 0.0 && h0 > dt) || (dt < 0.0 && h0 < dt)) {
337 step_status=this->stepp->step(t0,h0,nvar,yout,dydx,yout,yerr,
362 if (fabs(h0) < fabs(h_old) && t_next != t_curr) {
377 if (try_step_again==
false) {
389 const size_t hadjust_status=con.hadjust
390 (nvar,this->stepp->get_order(),yout,yerr,dydx_out,h0);
392 if (hadjust_status ==
404 if (fabs(h0) < fabs(h_old) && t_next != t_curr) {
424 if (final_step==
false) {
461 virtual int astep(
double &x,
double xmax,
double &h,
462 size_t n, vec_y_t &y, vec_dydx_t &dydx_out,
463 vec_yerr_t &yerr, func_t &derivs) {
474 int ret=derivs(x,n,y,dydx_int);
475 if (ret!=0)
return ret;
477 ret=evolve_apply(x,xmax,x,h,n,y,dydx_int,yout_int,yerr,dydx_out,derivs);
480 for(
size_t i=0;i<n;i++) {
485 if (this->verbose>0) {
486 std::cout <<
"astep_gsl step:";
487 std::cout << x <<
" ";
488 for(
size_t j=0;j<n;j++) std::cout << y[j] <<
" ";
489 std::cout << std::endl;
507 size_t n, vec_y_t &y, vec_dydx_t &dydx,
508 vec_yerr_t &yerr, func_t &derivs) {
519 int ret=evolve_apply(x,xmax,x,h,n,y,dydx,yout_int,yerr,dydx_int,derivs);
522 for(
size_t i=0;i<n;i++) {
528 if (this->verbose>0) {
529 std::cout <<
"astep_gsl step:";
530 std::cout << x <<
" ";
531 for(
size_t j=0;j<n;j++) std::cout << y[j] <<
" ";
532 std::cout << std::endl;
557 virtual int astep_full(
double x,
double xmax,
double &x_out,
double &h,
558 size_t n, vec_y_t &y, vec_dydx_t &dydx,
559 vec_y_t &yout, vec_yerr_t &yerr,
560 vec_dydx_t &dydx_out, func_t &derivs) {
566 int ret=evolve_apply(x,xmax,x_out,h,n,y,dydx,yout,yerr,dydx_out,
569 if (this->verbose>0) {
570 std::cout <<
"astep_gsl step:";
571 std::cout << x_out <<
" ";
572 for(
size_t j=0;j<n;j++) std::cout << yout[j] <<
" ";
573 std::cout << std::endl;
581 #ifndef DOXYGEN_NO_O2NS
unsigned long int count
The number of steps.
The main O<span style='position: relative; top: 0.3em; font-size: 0.8em'>2</span>scl O$_2$scl names...
double a_dydt
Derivative scaling factor (default 0)
invalid argument supplied by user
double last_step
The size of the last step.
Control structure for astep_gsl.
virtual int astep_derivs(double &x, double xmax, double &h, size_t n, vec_y_t &y, vec_dydx_t &dydx, vec_yerr_t &yerr, func_t &derivs)
Make an adaptive integration step of the system derivs with derivatives.
static const size_t hadj_dec
Recommend step decrease.
unsigned long int failed_steps
The number of failed steps.
bool standard
Use standard or scaled algorithm (default true)
double a_y
Function scaling factor (default 1)
ode_control_gsl< vec_y_t, vec_dydx_t, vec_yerr_t > con
Control specification.
size_t msize
The size of the allocated vectors.
int evolve_apply(double t0, double t1, double &t, double &h, size_t nvar, vec_y_t &y, vec_dydx_t &dydx, vec_y_t &yout, vec_yerr_t &yerr, vec_dydx_t &dydx_out, func_t &derivs)
Apply the evolution for the next adaptive step.
virtual int astep_full(double x, double xmax, double &x_out, double &h, size_t n, vec_y_t &y, vec_dydx_t &dydx, vec_y_t &yout, vec_yerr_t &yerr, vec_dydx_t &dydx_out, func_t &derivs)
Make an adaptive integration step of the system derivs.
void vector_copy(const vec_t &src, vec2_t &dest)
Simple vector copy.
static const size_t hadj_nil
No adjustment required.
ubvector scale_abs
Scalings.
static const size_t hadj_inc
Recommend step increase.
std::string dtos(double x, int prec=6, bool auto_prec=false)
Convert a double to a string.
Adaptive stepper [abstract base].
virtual int astep(double &x, double xmax, double &h, size_t n, vec_y_t &y, vec_dydx_t &dydx_out, vec_yerr_t &yerr, func_t &derivs)
Make an adaptive integration step of the system derivs.
#define O2SCL_ERR(d, n)
Set an error with message d and code n.
vec_dydx_t dydx_int
Internal storage for dydx.
double eps_rel
Relative precision (default 0)
vec_y_t yout_int
Temporary storage for yout.
double eps_abs
Absolute precision (default )
problem with user-supplied function
int set_scale(size_t nscal, const svec_t &scale)
Set the scaling for each differential equation.
Adaptive ODE stepper (GSL)