Classes for non-adaptive integration are provided as descendents of o2scl::ode_step and classes for adaptive integration are descendants of o2scl::astep_base. To specify a set of functions to these classes, use a child of ode_funct for a generic vector type. The classes o2scl::ode_rkf45_gsl and o2scl::ode_rkck_gsl are reasonable general-purpose non-adaptive integrators and o2scl::astep_gsl is a good general-purpose adaptive method for non-stiff problems. For stiff ODE's, use o2scl::ode_bsimp_gsl (see the Stiff differential equations example ).
Solution of simple initial value problems is performed by o2scl::ode_iv_solve. This class uses an adaptive integrator (default is o2scl::astep_gsl) and does some bookkeeping to simplify the solution of initial value problems. The o2scl::ode_iv_solve class will give the final value of the functions at the endpoint, provide the functions on a user-specified grid, or it will tabulate the ODE solution for you using the grid chosen with the adaptive stepper. A example demonstrating the solution of initial value problems is given in the Ordinary differential equations example .
The solution of boundary-value problems is based on the abstract base class o2scl::ode_bv_solve. At the moment, a simple shooting method is the only implementation of this base class and is given in o2scl::ode_bv_shoot . An experimental multishooting class is given in o2scl::ode_bv_mshoot .
An application of linear solvers to solve finite-difference equations approximating a boundary value problem is given in o2scl::ode_it_solve . A example demonstrating the iterative solution of a boundary value problem is given in the Iterative solution of ODEs example .
Ordinary differential equations example
This example solves the differential equations defining the the Bessel and Airy functions with both the Cash-Karp and Prince-Dormand steppers. It demonstrates the use of o2scl::ode_rkck_gsl, o2scl::ode_rk8pd_gsl, o2scl::astep_gsl, and o2scl::ode_iv_solve and stores the results in o2scl::table objects .
The Bessel functions are defined by
The Bessel functions of the first kind,
are finite at the origin, and the example solves the
case, where
and
.
#include <boost/numeric/ublas/vector.hpp>
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/matrix_proxy.hpp>
#include <gsl/gsl_sf_bessel.h>
#include <gsl/gsl_sf_airy.h>
#include <gsl/gsl_sf_gamma.h>
#include <o2scl/test_mgr.h>
#include <o2scl/ode_funct.h>
#include <o2scl/ode_rkck_gsl.h>
#include <o2scl/ode_rk8pd_gsl.h>
#include <o2scl/astep_gsl.h>
#include <o2scl/ode_iv_solve.h>
#include <o2scl/table.h>
#ifdef O2SCL_HDF
#include <o2scl/hdf_file.h>
#include <o2scl/hdf_io.h>
#endif
using namespace std;
#ifdef O2SCL_HDF
#endif
typedef boost::numeric::ublas::matrix_row<ubmatrix> ubmatrix_row;
int derivs(
double x,
size_t nv,
const ubvector &y,
dydx[0]=y[1];
if (x==0.0) dydx[1]=0.0;
else dydx[1]=(-x*y[1]+(-x*x+alpha*alpha)*y[0])/x/x;
return 0;
}
int derivs2(
double x,
size_t nv,
const ubvector &y,
dydx[0]=y[1];
dydx[1]=y[0]*x;
return 0;
}
dydx[0]=y[1];
if (x==0.0) dydx[1]=0.0;
else dydx[1]=(-x*y[1]+(-x*x+alpha*alpha)*y[0])/x/x;
return 0;
}
int main(void) {
cout.setf(ios::scientific);
cout.setf(ios::showpos);
double x, dx=1.0e-1;
ubvector y(2), dydx(2), yout(2), yerr(2), dydx_out(2);
double alpha=1.0;
std::bind(derivs,std::placeholders::_1,std::placeholders::_2,
std::placeholders::_3,std::placeholders::_4,alpha);
ode_funct_solve_grid od3=
std::bind(derivs3,std::placeholders::_1,std::placeholders::_2,
std::placeholders::_3,std::placeholders::_4,alpha);
cout << "Bessel function, Cash-Karp: " << endl;
x=0.0;
y[0]=0.0;
y[1]=0.5;
derivs(x,2,y,dydx,alpha);
cout << " x J1(calc) J1(exact) rel. diff. "
<< "err" << endl;
while (x<1.0) {
ode.
step(x,dx,2,y,dydx,y,yerr,dydx,od);
x+=dx;
cout << x << " " << y[0] << " "
<< gsl_sf_bessel_J1(x) << " ";
cout << fabs((y[0]-gsl_sf_bessel_J1(x))/gsl_sf_bessel_J1(x)) << " ";
cout << yerr[0] << endl;
t.
test_rel(y[0],gsl_sf_bessel_J1(x),5.0e-5,
"rkck");
double line[5]={x,y[0],gsl_sf_bessel_J1(x),
fabs((y[0]-gsl_sf_bessel_J1(x))/gsl_sf_bessel_J1(x)),
yerr[0]};
}
cout << "Accuracy at end: "
<< fabs(y[0]-gsl_sf_bessel_J1(x))/gsl_sf_bessel_J1(x) << endl;
cout << endl;
cout << "Bessel function, Prince-Dormand: " << endl;
x=0.0;
y[0]=0.0;
y[1]=0.5;
derivs(x,2,y,dydx,alpha);
cout << " x J1(calc) J1(exact) rel. diff. "
<< "err" << endl;
while (x<1.0) {
ode2.
step(x,dx,2,y,dydx,y,yerr,dydx,od);
x+=dx;
cout << x << " " << y[0] << " "
<< gsl_sf_bessel_J1(x) << " ";
cout << fabs((y[0]-gsl_sf_bessel_J1(x))/gsl_sf_bessel_J1(x)) << " ";
cout << yerr[0] << endl;
t.
test_rel(y[0],gsl_sf_bessel_J1(x),5.0e-4,
"rk8pd");
double line[5]={x,y[0],gsl_sf_bessel_J1(x),
fabs((y[0]-gsl_sf_bessel_J1(x))/gsl_sf_bessel_J1(x)),
yerr[0]};
}
cout << "Accuracy at end: "
<< fabs(y[0]-gsl_sf_bessel_J1(x))/gsl_sf_bessel_J1(x) << endl;
cout << endl;
Bessel function with non-adaptive steppers
Note that with a Bessel function and a fixed step size, the Prince-Dormand stepper (even though of higher order than the Cash-Karp stepper) is actually less accurate, and seriously underestimates the error.
The Airy functions are defined by
This example solves for the Airy function of the first kind,
.
cout << "Airy function, Cash-Karp: " << endl;
x=0.0;
y[0]=1.0/pow(3.0,2.0/3.0)/gsl_sf_gamma(2.0/3.0);
y[1]=-1.0/pow(3.0,1.0/3.0)/gsl_sf_gamma(1.0/3.0);
derivs2(x,2,y,dydx);
cout << " x Ai(calc) Ai(exact) rel. diff. "
<< "err" << endl;
while (x<1.0) {
ode.
step(x,dx,2,y,dydx,y,yerr,dydx,od2);
x+=dx;
cout << x << " " << y[0] << " "
<< gsl_sf_airy_Ai(x,GSL_PREC_DOUBLE) << " ";
cout << fabs((y[0]-gsl_sf_airy_Ai(x,GSL_PREC_DOUBLE))/
gsl_sf_airy_Ai(x,GSL_PREC_DOUBLE)) << " ";
cout << yerr[0] << endl;
t.
test_rel(y[0],gsl_sf_airy_Ai(x,GSL_PREC_DOUBLE),1.0e-8,
"rkck");
double line[5]={x,y[0],gsl_sf_airy_Ai(x,GSL_PREC_DOUBLE),
fabs((y[0]-gsl_sf_airy_Ai(x,GSL_PREC_DOUBLE))/
gsl_sf_airy_Ai(x,GSL_PREC_DOUBLE)),
yerr[0]};
}
cout << "Accuracy at end: "
<< fabs(y[0]-gsl_sf_airy_Ai(x,GSL_PREC_DOUBLE))/
gsl_sf_airy_Ai(x,GSL_PREC_DOUBLE) << endl;
cout << endl;
cout << "Airy function, Prince-Dormand: " << endl;
x=0.0;
y[0]=1.0/pow(3.0,2.0/3.0)/gsl_sf_gamma(2.0/3.0);
y[1]=-1.0/pow(3.0,1.0/3.0)/gsl_sf_gamma(1.0/3.0);
derivs2(x,2,y,dydx);
cout << " x Ai(calc) Ai(exact) rel. diff. "
<< "err" << endl;
while (x<1.0) {
ode2.
step(x,dx,2,y,dydx,y,yerr,dydx,od2);
x+=dx;
cout << x << " " << y[0] << " "
<< gsl_sf_airy_Ai(x,GSL_PREC_DOUBLE) << " ";
cout << fabs((y[0]-gsl_sf_airy_Ai(x,GSL_PREC_DOUBLE))/
gsl_sf_airy_Ai(x,GSL_PREC_DOUBLE)) << " ";
cout << yerr[0] << endl;
t.
test_rel(y[0],gsl_sf_airy_Ai(x,GSL_PREC_DOUBLE),1.0e-14,
"rk8pd");
double line[5]={x,y[0],gsl_sf_airy_Ai(x,GSL_PREC_DOUBLE),
fabs((y[0]-gsl_sf_airy_Ai(x,GSL_PREC_DOUBLE))/
gsl_sf_airy_Ai(x,GSL_PREC_DOUBLE)),
yerr[0]};
}
cout << "Accuracy at end: "
<< fabs(y[0]-gsl_sf_airy_Ai(x,GSL_PREC_DOUBLE))/
gsl_sf_airy_Ai(x,GSL_PREC_DOUBLE) << endl;
cout << endl;
Airy function with non-adaptive steppers
Here the higher order stepper is more accurate.
cout.precision(5);
cout << "Adaptive stepper: " << endl;
astep_gsl<> ode3;
x=0.0;
y[0]=0.0;
y[1]=0.5;
cout << " x J1(calc) J1(exact) rel. diff.";
cout << " err_0 err_1" << endl;
int k=0;
while (x<10.0) {
int retx=ode3.astep(x,10.0,dx,2,y,dydx,yerr,od);
if (k%3==0) {
cout << retx << " " << x << " " << y[0] << " "
<< gsl_sf_bessel_J1(x) << " ";
cout << fabs((y[0]-gsl_sf_bessel_J1(x))/gsl_sf_bessel_J1(x)) << " ";
cout << yerr[0] << " " << yerr[1] << endl;
}
t.
test_rel(y[0],gsl_sf_bessel_J1(x),5.0e-3,
"astep");
t.
test_rel(y[1],0.5*(gsl_sf_bessel_J0(x)-gsl_sf_bessel_Jn(2,x)),
5.0e-3,"astep 2");
t.
test_rel(yerr[0],0.0,4.0e-6,
"astep 3");
t.
test_rel(yerr[1],0.0,4.0e-6,
"astep 4");
double line[6]={x,y[0],gsl_sf_bessel_J1(x),
fabs((y[0]-gsl_sf_bessel_J1(x))/gsl_sf_bessel_J1(x)),
yerr[0],yerr[1]};
k++;
}
cout << "Accuracy at end: "
<< fabs(y[0]-gsl_sf_bessel_J1(x))/gsl_sf_bessel_J1(x) << endl;
cout << endl;
cout << "Adaptive stepper with smaller tolerances: " << endl;
ode3.con.eps_abs=1.0e-8;
ode3.con.a_dydt=1.0;
x=0.0;
y[0]=0.0;
y[1]=0.5;
cout << " x J1(calc) J1(exact) rel. diff.";
cout << " err_0 err_1" << endl;
k=0;
while (x<10.0) {
int retx=ode3.astep(x,10.0,dx,2,y,dydx,yerr,od);
if (k%3==0) {
cout << retx << " " << x << " " << y[0] << " "
<< gsl_sf_bessel_J1(x) << " ";
cout << fabs((y[0]-gsl_sf_bessel_J1(x))/gsl_sf_bessel_J1(x)) << " ";
cout << yerr[0] << " " << yerr[1] << endl;
}
t.
test_rel(y[0],gsl_sf_bessel_J1(x),5.0e-3,
"astep");
t.
test_rel(y[1],0.5*(gsl_sf_bessel_J0(x)-gsl_sf_bessel_Jn(2,x)),
5.0e-3,"astep 2");
t.
test_rel(yerr[0],0.0,4.0e-8,
"astep 3");
t.
test_rel(yerr[1],0.0,4.0e-8,
"astep 4");
double line[6]={x,y[0],gsl_sf_bessel_J1(x),
fabs((y[0]-gsl_sf_bessel_J1(x))/gsl_sf_bessel_J1(x)),
yerr[0],yerr[1]};
k++;
}
cout << "Accuracy at end: "
<< fabs(y[0]-gsl_sf_bessel_J1(x))/gsl_sf_bessel_J1(x) << endl;
cout << endl;
cout << "Adaptive stepper, Prince-Dormand: " << endl;
ode3.set_step(ode2);
x=0.0;
y[0]=0.0;
y[1]=0.5;
cout << " x J1(calc) J1(exact) rel. diff.";
cout << " err_0 err_1" << endl;
k=0;
while (x<10.0) {
int retx=ode3.astep(x,10.0,dx,2,y,dydx,yerr,od);
if (k%3==0) {
cout << retx << " " << x << " " << y[0] << " "
<< gsl_sf_bessel_J1(x) << " ";
cout << fabs((y[0]-gsl_sf_bessel_J1(x))/gsl_sf_bessel_J1(x)) << " ";
cout << yerr[0] << " " << yerr[1] << endl;
}
t.
test_rel(y[0],gsl_sf_bessel_J1(x),5.0e-3,
"astep");
t.
test_rel(y[1],0.5*(gsl_sf_bessel_J0(x)-gsl_sf_bessel_Jn(2,x)),
5.0e-3,"astep");
t.
test_rel(yerr[0],0.0,4.0e-8,
"astep 3");
t.
test_rel(yerr[1],0.0,4.0e-8,
"astep 4");
double line[6]={x,y[0],gsl_sf_bessel_J1(x),
fabs((y[0]-gsl_sf_bessel_J1(x))/gsl_sf_bessel_J1(x)),
yerr[0],yerr[1]};
k++;
}
cout << "Accuracy at end: "
<< fabs(y[0]-gsl_sf_bessel_J1(x))/gsl_sf_bessel_J1(x) << endl;
cout << endl;
Bessel function with adaptive steppers
cout.precision(6);
cout << "Initial value solver: " << endl;
ode_iv_solve_grid<> ode4;
ode4.ntrial=10000;
const size_t ngrid=101;
ubvector xg(ngrid), yinit(2);
ubmatrix yg(ngrid,2), ypg(ngrid,2), yerrg(ngrid,2);
for(size_t i=0;i<ngrid;i++) xg[i]=((double)i)/10.0;
xg[0]=0.0;
xg[ngrid-1]=10.0;
yg(0,0)=0.0;
yg(0,1)=0.5;
ode4.solve_grid<ubvector,ubmatrix>(0.1,2,ngrid,xg,yg,yerrg,ypg,od3);
cout << " x J1(calc) J1(exact) rel. diff." << endl;
for(size_t i=1;i<ngrid;i++) {
if (i%10==0) {
cout << xg[i] << " " << yg(i,0) << " "
<< gsl_sf_bessel_J1(xg[i]) << " ";
cout << fabs(yg(i,0)-gsl_sf_bessel_J1(xg[i])) << " "
<< fabs(yerrg(i,0)) << endl;
}
t.
test_rel(yg(i,0),gsl_sf_bessel_J1(xg[i]),5.0e-7,
"astep");
t.
test_rel(yg(i,1),0.5*(gsl_sf_bessel_J0(xg[i])-
gsl_sf_bessel_Jn(2,xg[i])),5.0e-7,"astep 2");
double line[5]={xg[i],yg(i,0),gsl_sf_bessel_J1(xg[i]),
fabs(yg(i,0)-gsl_sf_bessel_J1(xg[i])),
fabs(yerrg(i,0))};
}
cout << "Accuracy at end: "
<< fabs(yg(ngrid-1,0)-gsl_sf_bessel_J1(xg[ngrid-1]))/
gsl_sf_bessel_J1(xg[ngrid-1]) << endl;
cout << endl;
#ifdef O2SCL_HDF
hdf_file hf;
hf.open_or_create("ex_ode.o2");
for(size_t i=0;i<8;i++) {
hdf_output(hf,tab[i],((
string)
"table_")+
itos(i));
}
hf.close();
#endif
cout.unsetf(ios::showpos);
return 0;
}
Bessel function with ode_iv_solve
Stiff differential equations example
This example solves the differential equations
which have the exact solution
using both the stiff stepper o2scl::ode_bsimp_gsl and the standard adaptive stepper o2scl::astep_gsl . The relative error on the adaptive stepper is orders of magnitude larger.
Comparison of non-stiff and stiff ODE integrators
#include <boost/numeric/ublas/vector.hpp>
#include <boost/numeric/ublas/matrix.hpp>
#include <o2scl/test_mgr.h>
#include <o2scl/funct.h>
#include <o2scl/ode_funct.h>
#include <o2scl/astep_gsl.h>
#include <o2scl/ode_bsimp_gsl.h>
#include <o2scl/table.h>
#ifdef O2SCL_HDF
#include <o2scl/hdf_file.h>
#include <o2scl/hdf_io.h>
#endif
using namespace std;
#ifdef O2SCL_HDF
#endif
int derivs(
double x,
size_t nv,
const ubvector &y,
dydx[0]=480.0*y[0]+980.0*y[1];
dydx[1]=-490.0*y[0]-990.0*y[1];
return 0;
}
int jac(
double x,
size_t nv,
const ubvector &y,
dfdy(0,0)=480.0;
dfdy(0,1)=980.0;
dfdy(1,0)=-490.0;
dfdy(1,1)=-990.0;
dfdx[0]=0.0;
dfdx[1]=0.0;
return 0;
}
int main(void) {
cout.setf(ios::scientific);
cout.precision(3);
ubvector y1(2), dydx1(2), yout1(2), yerr1(2), dydx_out1(2);
y1[0]=1.0;
y1[1]=0.0;
for(size_t i=1;i<=40;i++) {
gb.
step(
x1,dx,2,y1,dydx1,y1,yerr1,dydx1,od11,oj);
double exact[2]={-exp(-500.0*
x1)+2.0*exp(-10.0*
x1),
exp(-500.0*
x1)-exp(-10.0*
x1)};
cout.setf(ios::showpos);
cout <<
x1 <<
" " << y1[0] <<
" " << y1[1] <<
" "
<< yerr1[0] << " " << yerr1[1] << " " << exact[0] << " "
<< exact[1] << endl;
cout.unsetf(ios::showpos);
double line[5]={
x1,y1[0],exact[0],fabs(yerr1[0]/exact[0]),
fabs((y1[0]-exact[0])/exact[0])};
}
cout << endl;
ubvector y2(2), dydx2(2), yout2(2), yerr2(2), dydx_out2(2);
y2[0]=1.0;
y2[1]=0.0;
size_t j=0;
ga.
astep(
x2,4.0,dx,2,y2,dydx2,yerr2,od11);
double exact[2]={-exp(-500.0*
x2)+2.0*exp(-10.0*
x2),
exp(-500.0*
x2)-exp(-10.0*
x2)};
if (j%25==0) {
cout.setf(ios::showpos);
cout <<
x2 <<
" " << y2[0] <<
" " << y2[1] <<
" "
<< yerr2[0] << " " << yerr2[1] << " " << exact[0] << " "
<< exact[1] << endl;
cout.unsetf(ios::showpos);
}
j++;
double line[5]={
x2,y2[0],exact[0],fabs(yerr2[0]/exact[0]),
fabs((y2[0]-exact[0])/exact[0])};
}
cout << endl;
#ifdef O2SCL_HDF
for(size_t i=0;i<2;i++) {
hdf_output(hf,tab[i],((
string)
"table_")+
itos(i));
}
#endif
return 0;
}
Iterative solution of ODEs example
This example solves the ODEs
given the boundary conditions
by linearizing the ODEs on a mesh and using an iterative method (sometimes called relaxation). The o2scl::ode_it_solve class demonstrates how this works, but is slow for large grid sizes because the matrix is very sparse.
#include <o2scl/ode_it_solve.h>
#include <o2scl/ode_iv_solve.h>
#include <o2scl/linear_solver.h>
using namespace std;
class ode_system {
public:
double left(size_t ieq, double x, ovector_base &yleft) {
if (ieq==0) return yleft[0]-1.0;
return yleft[1]*yleft[1]+yleft[2]*yleft[2]-2.0;
}
double right(size_t ieq, double x, ovector_base &yright) {
return yright[1]-3.0;
}
double derivs(size_t ieq, double x, ovector_base &y) {
if (ieq==1) return y[0]+y[1];
else if (ieq==2) return y[0]+y[2];
return y[1];
}
int shoot_derivs(double x, size_t nv, const ovector_base &y,
ovector_base &dydx) {
dydx[0]=y[1];
dydx[1]=y[0]+y[1];
dydx[2]=y[0]+y[2];
return 0;
}
};
int main(void) {
cout.setf(ios::scientific);
omatrix_base,omatrix_row,ovector_base,omatrix_base> oit;
ode_system os;
ode_it_funct_mfptr<ode_system> f_d(&os,&ode_system::derivs);
ode_it_funct_mfptr<ode_system> f_l(&os,&ode_system::left);
ode_it_funct_mfptr<ode_system> f_r(&os,&ode_system::right);
size_t ngrid=40;
size_t neq=3;
size_t nbleft=2;
ovector x(ngrid);
omatrix y(ngrid,neq);
for(size_t i=0;i<ngrid;i++) {
x[i]=((double)i)/((double)(ngrid-1));
y[i][0]=1.0+x[i]+1.0;
y[i][1]=3.0*x[i];
y[i][2]=-0.1*x[i]-1.4;
}
int pa=0;
omatrix A(ngrid*neq,ngrid*neq);
ovector rhs(ngrid*neq), dy(ngrid*neq);
oit.
solve(ngrid,neq,nbleft,x,y,f_d,f_l,f_r,A,rhs,dy);
ode_funct_mfptr<ode_system> f_sd(&os,&ode_system::shoot_derivs);
ovector ystart(neq), yend(neq);
for(size_t i=0;i<neq;i++) ystart[i]=y[0][i];
t.
test_rel(y[ngrid-1][0],yend[0],1.0e-3,
"yb");
t.
test_rel(y[0][1],0.25951,1.0e-3,
"yc");
t.
test_rel(y[ngrid-1][1],yend[1],1.0e-3,
"yd");
t.
test_rel(y[0][2],-1.3902,1.0e-3,
"ye");
t.
test_rel(y[ngrid-1][2],yend[2],1.0e-3,
"yf");
return 0;
}