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 108 typedef boost::numeric::ublas::unbounded_array<double> ubarr;
112 #ifndef DOXYGEN_INTERAL 156 vec_t y_extrap_sequence;
181 w[i]=(u>0.0) ? u : 1.0;
193 const double safety_f=0.25;
194 const double small_eps=safety_f*eps2;
201 {2,6,10,14,22,34,50,70};
205 a_work[0]=bd_sequence[0]+1.0;
208 a_work[k+1]=a_work[k]+bd_sequence[k+1];
214 const double tmp1=a_work[k+1]-a_work[i+1];
215 const double tmp2=(a_work[i+1]-a_work[0]+1.0)*(2*k+1);
216 alpha[k][i]=pow(small_eps,tmp1/tmp2);
220 a_work[0]+=dimension;
223 a_work[k+1]=a_work[k]+bd_sequence[k+1];
226 for (k=0;k<sequence_max-1;k++) {
227 if (a_work[k+2]>a_work[k+1]*alpha[k][k+1]) {
246 const unsigned int i_step,
const double x_i,
247 const vec_t &y_i, vec_t &y_0, vec_t &y_0_err,
256 for (j=0;j<
dim;j++) {
264 for (k=0;k<i_step;k++) {
265 double deltaloc=1.0/(x[i_step-k-1]-x_i);
266 const double f1=deltaloc*x_i;
267 const double f2=deltaloc*x[i_step-k-1];
269 for (j=0;j<
dim;j++) {
270 const double q_kj=dloc(k,j);
271 dloc(k,j)=y_0_err[j];
272 deltaloc=work[j]-q_kj;
273 y_0_err[j]=f1*deltaloc;
279 for (j=0;j<
dim;j++) {
280 dloc(i_step,j)=y_0_err[j];
295 const unsigned int n_step,
const ubarr &y,
296 const ubarr &yp_loc,
const vec_t &dfdt_loc,
297 const mat_t &dfdy_loc, vec_t &y_out) {
299 const double h=h_total/n_step;
309 const double max_sum=100.0*
dim;
316 for (i=0;i<
dim;i++) {
317 for (j=0;j<
dim;j++) {
318 a_mat(i,j)=-h*dfdy_loc(i,j);
334 for (i=0;i<
dim;i++) {
335 y_temp[i]=h*(yp_loc[i]+h*dfdt_loc[i]);
341 for (i=0;i<
dim;i++) {
342 const double di=delta_temp[i];
345 sum+=fabs(di)/weight[i];
353 status=(*funcp)(t,
dim,y_temp,y_out);
358 for (n_inter=1;n_inter<n_step;n_inter++) {
360 for (i=0;i<
dim;i++) {
361 rhs_temp[i]=h*y_out[i]-delta[i];
367 for (i=0;i<
dim;i++) {
368 delta[i]+=2.0*delta_temp[i];
370 sum+=fabs(delta[i])/weight[i];
378 status=(*funcp)(t,
dim,y_temp,y_out);
387 for (i=0;i<
dim;i++) {
388 rhs_temp[i]=h*y_out[i]-delta[i];
394 for (i=0;i<
dim;i++) {
395 y_out[i]=y_temp[i]+delta_temp[i];
396 sum+=fabs(delta_temp[i])/weight[i];
413 d.resize(sequence_max,n);
420 for(
size_t i=0;i<n;i++) yp[i]=0.0;
424 y_extrap_save.resize(n);
425 extrap_work.resize(n);
430 y_extrap_sequence.resize(n);
433 delta_temp.resize(n);
437 double sqrt_dbl_eps=sqrt(std::numeric_limits<double>::epsilon());
442 k_choice=k_choice_loc;
443 k_current=k_choice_loc;
444 order=2*k_choice_loc;
454 extrap_work.resize(0);
455 y_extrap_save.resize(0);
482 for(
size_t i=0;i<
dim;i++) yp[i]=0.0;
499 virtual int step(
double x,
double h,
size_t n, vec_t &y, vec_t &dydx,
500 vec_t &yout, vec_t &yerr, vec_t &dydx_out,
501 func_t &derivs, jac_func_t &jac) {
511 static const int bd_sequence[
sequence_count]={2,6,10,14,22,34,50,70};
517 if (h+t_local==t_local) {
521 O2SCL_ERR(
"Stepsize underflow in ode_bsimp_gsl::step().",
534 ret=jac(t_local,dim,y,dfdy,dfdt);
544 for (k=0;k<=k_current;k++) {
546 const unsigned int N=bd_sequence[k];
547 const double r=(h/N);
548 const double x_k=r*r;
553 int status=
step_local(t_local,h,N,y_extrap_save,yp,
554 dfdt,dfdy,y_extrap_sequence);
557 std::cout <<
"ode_bsimp_gsl: " << k <<
"/" << k_current <<
" " 558 <<
"status=" << status << std::endl;
565 for (i=0;i<
dim;i++) {
576 poly_extrap(d,ex_wk,k,x_k,y_extrap_sequence,y,yerr,extrap_work);
581 ret=derivs(t_local+h,dim,y,dydx_out);
595 #ifndef DOXYGEN_NO_O2NS static const int sequence_count
Number of entries in the Bader-Deuflhard extrapolation sequence.
The main O<span style='position: relative; top: 0.3em; font-size: 0.8em'>2</span>scl O$_2$scl names...
std::function< int(double, size_t, const boost::numeric::ublas::vector< double > &, boost::numeric::ublas::matrix< double > &, boost::numeric::ublas::vector< double > &) > ode_jac_funct
Functor for ODE Jacobians.
size_t dim
Size of allocated vectors.
int verbose
Verbose parameter.
int reset()
Reset stepper.
A class for representing permutations.
int LU_decomp(const size_t N, mat_t &A, o2scl::permutation &p, int &signum)
Compute the LU decomposition of the matrix A.
virtual int step(double x, double h, size_t n, vec_t &y, vec_t &dydx, vec_t &yout, vec_t &yerr, vec_t &dydx_out, func_t &derivs, jac_func_t &jac)
Perform an integration step.
size_t order
Order of last step.
void vector_copy(const vec_t &src, vec2_t &dest)
Simple vector copy.
int poly_extrap(ubmatrix &dloc, const double x[], const unsigned int i_step, const double x_i, const vec_t &y_i, vec_t &y_0, vec_t &y_0_err, ubarr &work)
Polynomial extrapolation.
ubmatrix d
Workspace for extrapolation.
double ex_wk[sequence_max]
Workspace for extrapolation.
std::function< int(double, size_t, const boost::numeric::ublas::vector< double > &, boost::numeric::ublas::vector< double > &)> ode_funct
Ordinary differential equation function.
ubmatrix a_mat
Workspace for linear system matrix.
int allocate(size_t n)
Allocate memory for a system of size n.
int compute_weights(const ubarr &y, ubarr &w, size_t n)
Compute weights.
int step_local(const double t0, const double h_total, const unsigned int n_step, const ubarr &y, const ubarr &yp_loc, const vec_t &dfdt_loc, const mat_t &dfdy_loc, vec_t &y_out)
Basic implicit Bulirsch-Stoer step.
#define O2SCL_ERR(d, n)
Set an error with message d and code n.
static const int sequence_max
Largest index of the Bader-Deuflhard extrapolation sequence.
jac_func_t * jfuncp
Jacobian.
Bulirsch-Stoer implicit ODE stepper (GSL)
int LU_solve(const size_t N, const mat_t &LU, const o2scl::permutation &p, const vec_t &b, vec2_t &x)
Solve a linear system after LU decomposition.
void free()
Free allocated memory.
func_t * funcp
Function specifying derivatives.
size_t deuf_kchoice(double eps2, size_t dimension)
Calculate a choice for the "order" of the method, using the Deuflhard criteria.