Using odeint function definition

Question

Pretty noob question so please bear with me.

I am following the example given here--> http://www.codeproject.com/Articles/268589/odeint-v2-Solving-ordinary-differential-equations

In particular, I am looking at this function:

void lorenz( state_type &x , state_type &dxdt , double t )
{                                                         
dxdt[0] = sigma * ( x[1] - x[0] );
dxdt[1] = R * x[0] - x[1] - x[0] * x[2];
dxdt[2] = x[0]*x[1] - b * x[2];
 } 

In my case, R takes on a series of values (vector with 100 doubles).

odeint is called as:

integrate_const( runge_kutta4< state_type >() , lorenz , x , 0.0 , 10.0 , dt );

I would like to do this for each value of R. How can I accomplish this? My knowledge of C++/OOP is limited, but I am willing to learn.

Thank you.

Answer 1

You can use the "class" version, but modify it so that it is initialized with the R value of interest to you.

class lorenz_class {
    double R_;
public:
    lorenz_class (double r) : R_(r) {}
    void operator()( state_type &x , state_type &dxdt , double t ) {
        dxdt[0] = sigma * ( x[1] - x[0] );
        dxdt[1] = R_ * x[0] - x[1] - x[0] * x[2];
        dxdt[2] = x[0]*x[1] - b * x[2];
    }
};

Then, iterate over your vector R and pass in the value to the lorenz_class instance that you pass to the integrate_const template function.

for (unsigned i = 0; i < myR.size(); ++i) {
    lorenz_class lorenz(myR[i]);
    integrate_const( runge_kutta4< state_type >() , lorenz , x , 0.0 , 10.0 , dt );
}

Answer 2

Just a little side note: The tutorial shows a very similar example of a parameter study of the Lorenz system: http://headmyshoulder.github.com/odeint-v2/doc/boost_numeric_odeint/tutorial.html. It is in the Thrust and the VexCL section and it shows how you can parallelize this problem to work on multiple CPUs or the GPU.