'std::bad_alloc using boost library odeint

Hi I am new to c++ and boost but i was asked to solve some ode's. I tried using a stiff method, it compiled but when i run the program it says

terminate called after throwing an instance of 'std::bad_alloc'
I searched for the error but since i am really amateur i couldn't find a solution


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;
    }
}


}
On line 52 you allocate about 20 giga byte memory (filled with 0.0). See:

http://www.cplusplus.com/reference/vector/vector/vector/

(2) fill constructor
Right so what i did was create 20e8 vectors with 0 values. I thought that those were the initial conditions for x[0],x[1]. But I don't want those to have the same value. How can I initiate the calculations from x[0] = 20e8 and x[1] = 0.1 ?

Also here do you know if the numbers stand for error or somethng like that?

auto stepper = make_dense_output< rosenbrock4< double > >( 1.0e-12 , 1.0e-12 );

In the jacobian why are they setting equal to zero?
dfdt[0] = 0.0; dfdt[1] = 0.0;

Thanks again, trying to connect the dots here




How can I initiate the calculations from x[0] = 20e8 and x[1] = 0.1 ?
Use curly braces instead of parentheses. But you don't need a vector if you use just two elements. You may use an arrray or std::array.

Also here do you know if the numbers stand for error or somethng like that?
Just very small numbers. Basically next to 0 but not 0.

In the jacobian why are they setting equal to zero?
They are not used at all. So no I don't know why they exist.
Topic archived. No new replies allowed.