1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80
|
#include <iostream>
#include <fstream>
#include <utility>
#include <cmath>
#include <boost/numeric/odeint.hpp>
#include <boost/phoenix/core.hpp>
#include <boost/phoenix/core.hpp>
#include <boost/phoenix/operator.hpp>
using namespace std;
using namespace boost::numeric::odeint;
namespace phoenix = boost::phoenix;
const double b = 43.0e17;
//[ stiff_system_definition
typedef boost::numeric::ublas::vector< double > vector_type;
typedef boost::numeric::ublas::matrix< double > matrix_type;
struct stiff_system
{
void operator()( const vector_type &x , vector_type &dxdt , double /* t */ )
{
dxdt[ 0 ] = -b*(64.0/5)*(1 + (73.0/24)*pow(x[1],2) + (37.0/96)*pow(x[1],4) )/pow(x[0],3)*pow(1-pow(x[1],2),7.0/2);
dxdt[ 1 ] = -b*(304.0/15)*x[1]*(1 + (121.0/304)*pow(x[1],2))/ pow(x[0],4)*pow((1 - pow(x[1],2)),5.0/2);
}
};
struct stiff_system_jacobi
{
void operator()( const vector_type &x , matrix_type &J , const double /* t */, vector_type &dfdt )
{
J( 0 , 0 ) = -(3.0/x[0])*(-b*(64.0/5)*(1 + (73.0/24)*pow(x[1],2) + (37.0/96)*pow(x[1],4) )/pow(x[0],3)*pow(1-pow(x[1],2),7.0/2));
J( 0 , 1 ) = -(b/pow(x[0],3))*(64.0/5)*(111*pow(x[1],5) + 1608*pow(x[1],3) + 1256*x[1])/ 96*(pow(1 - pow(x[1],2),9.0/2));
J( 1 , 0 ) = -(4.0/x[0])*(-b*(304.0/15)*x[1]*(1 + (121.0/304)*pow(x[1],2))/ pow(x[0],4)*pow((1 - pow(x[1],2)),5.0/2));
J( 1 , 1 ) = -(-b/pow(x[0],4))*(304.0/96)*(363*pow(x[1],2) + 1762*x[1])/ (304*pow(1- x[1]*x[1],7.0/2));
dfdt[0] = 0.0;
dfdt[1] = 0.0;
}
};
int main( int argc , char **argv )
{
auto stepper = make_dense_output< rosenbrock4< double > >( 1.0e-12 , 1.0e-12 );
auto ode = make_pair( stiff_system() , stiff_system_jacobi() );
double t = 0.0;
double const end_time = 9e16;
double const dt = 1.0;
vector_type x( 20e+08 , 0.0 );
const double y_min = 7.0e07;
stepper.initialize( x , t , dt );
cout << t << " " << x[0] << " " << x[1] << endl; // or some other output
t += dt;
while( t < end_time )
{
if( t > stepper.current_time() )
{
// perform a real step
stepper.do_step( ode );
}
else
{
// perform a dense output step
stepper.calc_state( t , x );
cout << t << " " << x[0] << " " << x[1] << endl; // or some other output
t += dt;
}
if( x[1] < y_min ) // or some other condition
{
cout << "Bound reached." << endl;
break;
}
}
}
|