Go to the documentation of this file.
23 #ifndef O2SCL_ODE_IV_SOLVE_H
24 #define O2SCL_ODE_IV_SOLVE_H
32 #include <boost/numeric/ublas/vector.hpp>
33 #include <boost/numeric/ublas/matrix.hpp>
34 #include <boost/numeric/ublas/matrix_proxy.hpp>
36 #include <o2scl/astep.h>
37 #include <o2scl/astep_gsl.h>
39 #ifndef DOXYGEN_NO_O2NS
109 #ifndef DOXYGEN_INTERNAL
115 vec_t vtemp, vtemp2, vtemp3, vtemp4;
126 std::cout <<
type() <<
" x: " << x <<
" y: ";
127 for(
size_t i=0;i<nv;i++) std::cout << y[i] <<
" ";
128 std::cout << std::endl;
173 virtual ~ode_iv_solve() {
198 vec_t &ystart, vec_t ¥d, func_t &derivs) {
221 vec_t &ystart, vec_t ¥d, vec_t &yerr,
256 vec_t &ystart, vec_t ¥d, vec_t &yerr,
257 vec_t &dydx_end, func_t &derivs) {
259 if ((
x1>x0 && h<=0.0) || (x0>
x1 && h>=0.0)) {
260 std::string str=
"Interval direction (x1-x0="+
o2scl::dtos(
x1-x0)+
261 ") does not match step direction (h="+
o2scl::dtos(h)+
262 " in ode_iv_solve::solve_final_value().";
266 O2SCL_ERR2(
"Starting and final endpoints identical in ",
267 "ode_iv_solve::solve_final_value().",
exc_einval);
277 double x=x0, xverb=0.0, dxverb=0.0, xnext;
278 int ret=0, first_ret=0;
283 vec_t &dydx_start=vtemp;
303 for(
size_t i=0;i<n;i++) yend[i]=ystart[i];
306 derivs(x,n,yend,dydx_start);
309 while (done==
false) {
312 ret=
astp->astep_full(x,
x1,xnext,h,n,ystart,dydx_start,
313 yend,yerr,dydx_end,derivs);
318 "ode_iv_solve::solve_final_value()",ret);
319 }
else if (first_ret!=0) {
325 if (
verbose>0 && xnext>xverb) {
337 std::string str=
"Too many steps required (nsteps="+
339 ", x="+
o2scl::dtos(x)+
") in ode_iv_solve::solve_final_value().";
347 if (xnext>=
x1) done=
true;
349 if (xnext<=
x1) done=
true;
356 ret=
astp->astep_full(xnext,
x1,x,h,n,yend,dydx_end,ystart,yerr,
362 "ode_iv_solve::solve_final_value()",ret);
363 }
else if (first_ret!=0) {
381 std::string str=
"Too many steps required (ntrial="+
itos(
ntrial)+
382 ", x="+
o2scl::dtos(x)+
") in ode_iv_solve::solve_final_value().";
392 for(
size_t j=0;j<n;j++) {
394 dydx_end[j]=dydx_start[j];
400 for(
size_t j=0;j<n;j++) {
402 dydx_end[j]=dydx_start[j];
453 template<
class mat_t>
455 vec_t &x_sol, mat_t &y_sol, mat_t &yerr_sol,
456 mat_t &dydx_sol, func_t &derivs,
size_t istart=0) {
467 double dx_tab=(
x1-x0)/((
double)(n_sol-istart-1));
469 double x_verb=x0+dx_verb;
470 double x_tab=x0+dx_tab;
477 vec_t &ystart=vtemp4;
478 vec_t &dydx_start=vtemp2;
481 for(
size_t j=0;j<n;j++) ystart[j]=y_sol(istart,j);
493 derivs(x0,n,ystart,dydx_start);
497 for(
size_t j=0;j<n;j++) {
498 dydx_sol(istart,j)=dydx_start[j];
499 yerr_sol(istart,j)=0.0;
503 size_t icurr=istart+1;
505 for(
size_t j=0;j<n;j++) y_sol(icurr,j)=ystart[j];
510 while (done==
false) {
515 vec_t &dydx_out=vtemp3;
520 for(
size_t i=0;i<n;i++) yrow[i]=y_sol(icurr,i);
521 ret=
astp->astep_full(x0,
x1,xnext,h,n,yrow,dydx_start,ystart,yerr,
528 O2SCL_ERR2(
"Adaptive stepper returned non-zero in ",
530 }
else if (first_ret==0) {
538 std::string str=
"Too many steps required (ntrial="+
itos(
ntrial)+
539 ", x="+
o2scl::dtos(x0)+
") in ode_iv_solve::solve_store().";
544 if (xnext>=x_verb &&
verbose>0) {
561 for(
size_t j=0;j<n;j++) {
562 y_sol(icurr,j)=ystart[j];
563 dydx_sol(icurr,j)=dydx_out[j];
564 yerr_sol(icurr,j)=yerr[j];
572 if (xnext>=x_tab && icurr<nmax) {
576 for(
size_t j=0;j<n;j++) {
577 y_sol(icurr,j)=ystart[j];
578 dydx_sol(icurr,j)=dydx_out[j];
579 yerr_sol(icurr,j)=yerr[j];
582 y_sol(icurr+1,j)=ystart[j];
583 dydx_start[j]=dydx_out[j];
596 for(
size_t j=0;j<n;j++) {
597 dydx_start[j]=dydx_out[j];
598 y_sol(icurr,j)=ystart[j];
647 virtual const char *
type() {
return "ode_iv_solve"; }
651 typedef matrix_row_gen_ctor<boost::numeric::ublas::matrix<double> >
654 typedef std::function<int(
double,
size_t,
const solve_grid_mat_row &,
655 solve_grid_mat_row &)>
656 ode_funct_solve_grid;
671 template<
class func_t=ode_funct_solve_grid,
672 class mat_row_t=solve_grid_mat_row>
675 #ifndef DOXYGEN_INTERNAL
684 std::cout <<
type() <<
" x: " << x <<
" y: ";
685 for(
size_t i=0;i<nv;i++) std::cout << y[i] <<
" ";
686 std::cout << std::endl;
706 virtual ~ode_iv_solve_grid() {
732 template<
class vec_t,
class mat_t>
734 mat_t &ysol, mat_t &err_sol, mat_t &dydx_sol,
738 double x1=xsol[nsol-1];
741 int ret=0, first_ret=0;
747 mat_row_t y_start(n);
748 mat_row_t dydx_start(n);
751 for(
size_t j=0;j<n;j++) {
752 y_start[j]=ysol(0,j);
759 derivs(x0,n,y_start,dydx_start);
762 for(
size_t j=0;j<n;j++) {
763 dydx_sol(0,j)=dydx_start[j];
767 for(
size_t i=1;i<nsol && ret==0;i++) {
769 mat_row_t y_row(ysol,i);
770 mat_row_t dydx_row(dydx_sol,i);
771 mat_row_t yerr_row(err_sol,i);
775 while(done==
false && ret==0) {
777 ret=
astp->astep_full(x,xsol[i],xnext,h,n,y_start,dydx_start,
778 y_row,yerr_row,dydx_row,derivs);
784 "ode_iv_solve_grid::solve_grid()",ret);
785 }
else if (first_ret!=0) {
791 std::string str=
"Too many steps required (ntrial="+
itos(
ntrial)+
792 ", x="+
o2scl::dtos(x)+
") in ode_iv_solve_grid::solve_grid().";
798 for(
size_t j=0;j<n;j++) {
800 dydx_start[j]=dydx_row[j];
805 if (x>=xsol[i]) done=
true;
807 if (x<=xsol[i]) done=
true;
847 virtual const char *
type() {
return "ode_iv_solve_grid"; }
851 #ifndef DOXYGEN_NO_O2NS
int set_astep(astep_base< mat_row_t, mat_row_t, mat_row_t, func_t > &as)
Set the adaptive stepper to use.
Solve an initial-value ODE problems given an adaptive ODE stepper.
int solve_final_value(double x0, double x1, double h, size_t n, vec_t &ystart, vec_t ¥d, func_t &derivs)
Solve the initial-value problem to get the final value.
virtual const char * type()
Return the type, "ode_iv_solve".
int solve_store(double x0, double x1, double h, size_t n, size_t &n_sol, vec_t &x_sol, mat_t &y_sol, mat_t &yerr_sol, mat_t &dydx_sol, func_t &derivs, size_t istart=0)
Solve the initial-value problem and store the associated output.
std::string dtos(const fp_t &x, int prec=6, bool auto_prec=false)
Convert a double to a string.
astep_base< vec_t, vec_t, vec_t, func_t > * astp
The adaptive stepper.
@ exc_efailed
generic failure
int solve_grid(double h, size_t n, size_t nsol, vec_t &xsol, mat_t &ysol, mat_t &err_sol, mat_t &dydx_sol, func_t &derivs)
Solve the initial-value problem from x0 to x1 over a grid storing derivatives and errors.
size_t ntrial
Maximum number of applications of the adaptive stepper (default 1000)
#define O2SCL_ERR2(d, d2, n)
Set an error, two-string version.
size_t nsteps
Number of adaptive steps employed.
astep_gsl< vec_t, vec_t, vec_t, func_t > gsl_astp
The default adaptive stepper.
The main O<span style='position: relative; top: 0.3em; font-size: 0.8em'>2</span>scl O$_2$scl names...
int set_astep(astep_base< vec_t, vec_t, vec_t, func_t > &as)
Set the adaptive stepper to use.
@ exc_emaxiter
exceeded max number of iterations
int solve_final_value(double x0, double x1, double h, size_t n, vec_t &ystart, vec_t ¥d, vec_t &yerr, vec_t &dydx_end, func_t &derivs)
Solve the initial-value problem to get the final value, derivatives, and errors.
int solve_final_value(double x0, double x1, double h, size_t n, vec_t &ystart, vec_t ¥d, vec_t &yerr, func_t &derivs)
Solve the initial-value problem to get the final value with errors.
virtual const char * type()
Return the type, "ode_iv_solve".
virtual int print_iter(double x, size_t nv, mat_row_t &y)
Print out iteration information.
bool exit_on_fail
If true, stop the solution if the adaptive stepper fails (default true)
int verbose
Set output level.
static const double x1[5]
@ exc_einval
invalid argument supplied by user
#define O2SCL_CONV_RET(d, n, b)
Set a "convergence" error and return the error value.
size_t mem_size
The size of the temporary vectors.
std::function< int(double, size_t, const boost::numeric::ublas::vector< double > &, boost::numeric::ublas::vector< double > &)> ode_funct
Ordinary differential equation function in src/ode/ode_funct.h.
size_t ntrial
Maximum number of applications of the adaptive stepper (default 1000)
Adaptive ODE stepper (GSL)
bool err_nonconv
If true, call the error handler if the solution does not converge (default true)
size_t nsteps_out
Number of output points for verbose output (default 10)
bool err_nonconv
If true, call the error handler if the solution does not converge (default true)
int verbose
Set output level.
std::string itos(int x)
Convert an integer to a string.
virtual int print_iter(double x, size_t nv, vec_t &y)
Print out iteration information.
#define O2SCL_ERR(d, n)
Set an error with message d and code n.
Adaptive stepper [abstract base].
bool exit_on_fail
If true, stop the solution if the adaptive stepper fails (default true)
astep_base< mat_row_t, mat_row_t, mat_row_t, func_t > * astp
The adaptive stepper.
size_t nsteps
Number of adaptive ste!ps employed.
void free()
Free allocated memory.
void allocate(size_t n)
Allocate space for temporary vectors.
Solve an initial-value ODE problems on a grid given an adaptive ODE stepper.
std::string szttos(size_t x)
Convert a size_t to a string.
astep_gsl< mat_row_t, mat_row_t, mat_row_t, func_t > gsl_astp
The default adaptive stepper.
Documentation generated with Doxygen. Provided under the
GNU Free Documentation License (see License Information).