Class ode_bsimp_gsl (o2scl)¶

O2scl : Class List

template<class func_t = ode_funct, class jac_func_t = ode_jac_funct, class vec_t = boost::numeric::ublas::vector<double>, class mat_t = boost::numeric::ublas::matrix<double>> class ode_bsimp_gsl¶

Bulirsch-Stoer implicit ODE stepper (GSL)

Implicit Bulirsch-Stoer method of Bader and Deuflhard ( [Bader83] ).

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

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

Note

The variable h_next was defined in the original GSL version has been removed here, as it was unused by the stepper routine.

Note

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.

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¶

size_t order¶: Order of last step.

int verbose¶: Verbose parameter.

inline int compute_weights(const ubarr &y, ubarr &w, size_t n)¶: Compute weights.

inline size_t deuf_kchoice(double eps2, size_t dimension)¶

Calculate a choice for the “order” of the method, using the Deuflhard criteria.

Used in the allocate() function.

inline 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.

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.

inline 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.

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.

inline int allocate(size_t n)¶: Allocate memory for a system of size n.

inline void free()¶: Free allocated memory.

inline ode_bsimp_gsl()¶

inline virtual ~ode_bsimp_gsl()¶

inline int reset()¶: Reset stepper.

inline 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.

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.

Public Types

typedef boost::numeric::ublas::unbounded_array<double> ubarr¶

typedef boost::numeric::ublas::vector<double> ubvector¶

typedef boost::numeric::ublas::matrix<double> ubmatrix¶

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.

(This state variable was named ‘x’ in GSL.)

Protected Static 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.