Public Types | Public Member Functions | Public Attributes | Protected Member Functions | Protected Attributes | Static Protected Attributes | List of all members
o2scl::ode_bsimp_gsl< func_t, jac_func_t, vec_t, mat_t > Class Template Reference

Bulirsch-Stoer implicit ODE stepper (GSL) More...

#include <ode_bsimp_gsl.h>

Public Types

typedef boost::numeric::ublas::unbounded_array< double > ubarr
 
typedef boost::numeric::ublas::vector< double > ubvector
 
typedef boost::numeric::ublas::matrix< double > ubmatrix
 

Public Member Functions

int reset ()
 Reset stepper.
 
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. More...
 

Public Attributes

int verbose
 Verbose parameter.
 

Protected Member Functions

int compute_weights (const ubarr &y, ubarr &w, size_t n)
 Compute weights.
 
size_t deuf_kchoice (double eps2, size_t dimension)
 Calculate a choice for the "order" of the method, using the Deuflhard criteria. More...
 
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. More...
 
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. More...
 
int allocate (size_t n)
 Allocate memory for a system of size n.
 
void free ()
 Free allocated memory.
 

Protected Attributes

size_t dim
 Size of allocated vectors.
 
func_t * funcp
 Function specifying derivatives.
 
jac_func_t * jfuncp
 Jacobian.
 
ubmatrix d
 Workspace for extrapolation.
 
ubmatrix a_mat
 Workspace for linear system matrix.
 
double ex_wk [sequence_max]
 Workspace for extrapolation. More...
 
size_t order
 Order of last step.
 
State info
size_t k_current
 
size_t k_choice
 
double eps
 
Workspace for extrapolation step
ubarr yp
 
ubarr y_save
 
ubarr yerr_save
 
ubarr y_extrap_save
 
vec_t y_extrap_sequence
 
ubarr extrap_work
 
vec_t dfdt
 
vec_t y_temp
 
vec_t delta_temp
 
ubarr weight
 
mat_t dfdy
 
Workspace for the basic stepper
vec_t rhs_temp
 
ubarr delta
 

Static Protected Attributes

static const int sequence_count =8
 Number of entries in the Bader-Deuflhard extrapolation sequence.
 
static const int sequence_max =7
 Largest index of the Bader-Deuflhard extrapolation sequence.
 

Detailed Description

template<class func_t = ode_funct11, class jac_func_t = ode_jac_funct11, class vec_t = boost::numeric::ublas::vector<double>, class mat_t = boost::numeric::ublas::matrix<double>>
class o2scl::ode_bsimp_gsl< func_t, jac_func_t, vec_t, mat_t >

Bader-Deuflhard implicit extrapolative stepper (Bader83).

Note
The variable h_next was defined in the original GSL version has been removed here, as it was unused by the stepper routine.
At the moment, this class retains the GSL approach to handling non-integer return values in the user-specified derivative function. If the user-specified function returns exc_efailed, then the stepper attempts to decrease the stepsize to continue. If the user-specified function returns a non-zero value other than exc_efailed, or if the Jacobian evaluation returns any non-zero value, then the stepper aborts and returns the error value without calling the error handler. This behavior may change in the future.

There is an example for the usage of this class in examples/ex_stiff.cpp documented in the Stiff differential equations example section.

Idea for Future:
More detailed documentation about the algorithm
Idea for Future:
Figure out if this should be a child of ode_step or astep. The function step_local() is actually its own ODE stepper and could be reimplemented as an object of type ode_step.
Idea for Future:
I don't like setting yerr to GSL_POSINF, there should be a better way to force an adaptive stepper which is calling this stepper to readjust the stepsize.
Idea for Future:
The functions deuf_kchoice() and compute_weights() can be moved out of this header file.
Idea for Future:
Rework internal arrays
Idea for Future:
Rework linear solver to be amenable to using a sparse matrix solver

Definition at line 104 of file ode_bsimp_gsl.h.

Member Function Documentation

◆ deuf_kchoice()

template<class func_t = ode_funct11, class jac_func_t = ode_jac_funct11, class vec_t = boost::numeric::ublas::vector<double>, class mat_t = boost::numeric::ublas::matrix<double>>
size_t o2scl::ode_bsimp_gsl< func_t, jac_func_t, vec_t, mat_t >::deuf_kchoice ( double  eps2,
size_t  dimension 
)
inlineprotected

Used in the allocate() function.

Definition at line 191 of file ode_bsimp_gsl.h.

◆ poly_extrap()

template<class func_t = ode_funct11, class jac_func_t = ode_jac_funct11, class vec_t = boost::numeric::ublas::vector<double>, class mat_t = boost::numeric::ublas::matrix<double>>
int o2scl::ode_bsimp_gsl< func_t, jac_func_t, vec_t, mat_t >::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 
)
inlineprotected

Compute the step of index i_step using polynomial extrapolation to evaulate functions by fitting a polynomial to estimates (x_i,y_i) and output the result to y_0 and y_0_err.

The index i_step begins with zero.

Definition at line 245 of file ode_bsimp_gsl.h.

◆ step()

template<class func_t = ode_funct11, class jac_func_t = ode_jac_funct11, class vec_t = boost::numeric::ublas::vector<double>, class mat_t = boost::numeric::ublas::matrix<double>>
virtual int o2scl::ode_bsimp_gsl< func_t, jac_func_t, vec_t, mat_t >::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 
)
inlinevirtual

Given initial value of the n-dimensional function in y and the derivative in dydx (which must generally be computed beforehand) at the point x, take a step of size h giving the result in yout, the uncertainty in yerr, and the new derivative in dydx_out (at x+h) using function derivs to calculate derivatives. Implementations which do not calculate yerr and/or dydx_out do not reference these variables so that a blank vec_t can be given. All of the implementations allow yout=y and dydx_out=dydx if necessary.

Definition at line 499 of file ode_bsimp_gsl.h.

◆ step_local()

template<class func_t = ode_funct11, class jac_func_t = ode_jac_funct11, class vec_t = boost::numeric::ublas::vector<double>, class mat_t = boost::numeric::ublas::matrix<double>>
int o2scl::ode_bsimp_gsl< func_t, jac_func_t, vec_t, mat_t >::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 
)
inlineprotected

Divide the step h_total into n_step smaller steps and do the Bader-Deuflhard semi-implicit iteration. This function starts at t0 with function values y, derivatives yp_loc, and information from the Jacobian to compute the final value y_out.

Definition at line 294 of file ode_bsimp_gsl.h.

Member Data Documentation

◆ ex_wk

template<class func_t = ode_funct11, class jac_func_t = ode_jac_funct11, class vec_t = boost::numeric::ublas::vector<double>, class mat_t = boost::numeric::ublas::matrix<double>>
double o2scl::ode_bsimp_gsl< func_t, jac_func_t, vec_t, mat_t >::ex_wk[sequence_max]
protected

(This state variable was named 'x' in GSL.)

Definition at line 141 of file ode_bsimp_gsl.h.


The documentation for this class was generated from the following file:

Documentation generated with Doxygen. Provided under the GNU Free Documentation License (see License Information).