How to correctly use Odeint library for solving ODEs in C++?

Question

I am trying to integrate odeint for a physics simulation. I build a model containing all parameters for my equation:

(model.h)

#include <vector>
#include <boost/array.hpp>
#include <string>
using namespace std;
typedef boost::array< double , 2> state_type;

class model 
{
    public:
        model(double, double, double, double, double, double, double, double, double, double,
                double, double, double, double, double, double, double, double, double);
        void setData(vector<vector<double>>);
        void operator()(const state_type&, state_type&, const double);
        //void ode_metal(const state_type&, state_type&, const double);
    private:
        double T;
        double g;
        double n_trans;
        double tau1;
        double P1;
        double T1;
        double tau2;
        double P2;
        double T2;
        double alignment;
        double d;
        double L; 
        double tau_nonrad;
        double tau_rad; 
        double F;
        double G;
        double f_p;
        double beta;
        double A;
        double pulses(double);
        vector<double> delays;
        vector<double> signals;
        //vector<double> run();
        double evaluate(vector<double>);
};

Here is my model.cpp:

#include "model.h"
#include <boost/array.hpp>
#include <boost/numeric/odeint.hpp>

using namespace std;
using namespace boost::numeric::odeint;

typedef boost::array< double , 2> state_type;
typedef runge_kutta_dopri5< double > stepper_type;
const double pi = 3.141592653589793;
const double hbar = 1.0545718001391127e-34;

model::model(double T, double g, double n_trans, double tau1, double P1, double T1, double tau2, double P2, double T2, double alignment,
                double d, double L, double tau_nonrad, double tau_rad, double F, double G, double f_p, double beta, double A) {
    this->T = T;
    this->g = g;
    this->n_trans = n_trans;
    this->tau1 = tau1;
    this->P1 = P1;
    this->T1 = T1;
    this->tau2 = tau2;
    this->P2 = P2;
    this->T2 = T2;
    this->alignment = alignment;
    this->d = d;
    this->L = L; 
    this->tau_nonrad = tau_nonrad;
    this->tau_rad = tau_rad; 
    this->F = F;
    this->G = G;
    this->f_p = f_p;
    this->beta = beta;
    this->A = A;
 }

 double model::pulses(double t) {
    return (1-(this->alignment)*abs(t-(this->tau2))/100E-12)*(this->P1)*exp(-pow((t-(this->tau1))/(this->T1),2))+(this->P2)*exp(-pow((t-(this->tau2))/(this->T2),2));
}

void model::operator()(const state_type& x, state_type& dxdt, const double t){
    double P = pulses(t)*1E10;
    double P_n = P*L*d/(hbar*(2*pi*3E8/350E-9)*L*pi*(d/2)*(d/2));
    double n_nonrad = -x[0]/tau_nonrad;
    double n_sp = -G*F*beta*x[0]/tau_rad;
    double n_st = -(3E8/16)*g*log(x[0]/n_trans)*G*F*x[1];
    dxdt[0] = P_n+n_nonrad+n_sp+n_st;
    dxdt[1] = -n_sp - n_st - f_p*x[1];
}

void model::setData(vector<vector<double>> data) {
    for(int i=0; i<data.size();i++) {
        delays.push_back(data[i][0]);
        signals.push_back(data[i][1]);
    }
    /*
    cout << "Delays: " << std::endl;
    for(int i=0; i<delays.size();i++) {
        cout << delays[i] << endl;
    }
    cout << "Signals: " << std::endl;
    for(int i=0; i<signals.size();i++) {
        cout << signals[i] << endl;
    }
    */
}

And this is my main.cpp:

#include <iostream>
#include <boost/array.hpp>
#include <boost/numeric/odeint.hpp>
#include <vector>
#include "model.h"
#include <fstream>
#include <sstream>
#include <string>

using namespace std;
using namespace boost::numeric::odeint;

double T = 200E-9;
double g = 7E5;
double n_trans = 2E24;
double tau1 = 1E-10;
double P1 = 140;
double T1 = 100E-15;
double tau2 = 5E-9;
double P2 = 40;
double T2 = 100E-15;
double alignment = 0;
double d = 200*1E-9;
double L = 6E-6;
double tau_nonrad = 7E-10;
double tau_rad = 7E-10;
double F = 1;
double G = 0.92;
double f_p = 2E12;
double beta = 1;
double A = 1;
typedef boost::array< double , 2> state_type;
typedef runge_kutta_dopri5< double > stepper_type;


void write_cerr( const state_type &x , const double t )
{
    cout << t << '\t' << x[0] << '\t' << x[1] << endl;
}


vector<vector<double>> readData(string filename) {
    vector<vector<double>> data;
    vector<double> row;
    string line, word;

    fstream file (filename, ios::in);
    if(file.is_open())
    {
        while(getline(file, line))
        {
            row.clear();
 
            stringstream str(line);
 
            while(getline(str, word, ' ')) {
                row.push_back(stod(word));
            }
            data.push_back(row);
        }
    }
    else
        cout<<"Could not open the file\n";
    /*
    for(int i=0;i<data.size();i++)
    {
        for(int j=0;j<data[i].size();j++)
        {
            cout<<data[i][j]<<" ";
        }
        cout<<"\n";
    }
    */
    return data;
}

void write_cout( const state_type &x , double t )
{
    cout << t << '\t' << x[0] << '\t' << x[1] << endl;
}

int main() {
    model m(T, g, n_trans, tau1, P1, T1, tau2, P2, T2, alignment, d, L, tau_nonrad, tau_rad, F, G, f_p, beta, A);
    state_type x = { 1.0 , 0.0 };
    vector<vector<double>> data = readData("w1_dynamics.txt");
    m.setData(data);
    integrate_adaptive(make_controlled(1E-12, 1E-12, stepper_type()), m, x, 0.0, T, 1E-12, write_cout);
    return 0;
}

Somehow, this does not compile properly. The error is very long. Here are the first few lines:

error: function "model::operator()" cannot be called with the given argument list
            argument types are: (const state_type, double, double)
            object type is: model
          sys( x , m_dxdt.m_v , t );
          ^
model.h(15): note: this candidate was rejected because arguments do not match
          void operator()(const state_type&, state_type&, const double);
               ^
          detected during:
            instantiation of "void boost::numeric::odeint::controlled_runge_kutta<ErrorStepper, ErrorChecker, StepAdjuster, Resizer, boost::numeric::odeint::explicit_error_stepper_fsal_tag>::initialize(System, const StateIn &, boost::numeric::odeint::controlled_runge_kutta<ErrorStepper, ErrorChecker, StepAdjuster, Resizer, boost::numeric::odeint::explicit_error_stepper_fsal_tag>::time_type) [with ErrorStepper=boost::numeric::odeint::runge_kutta_dopri5<double, double, double, double,
                      boost::numeric::odeint::vector_space_algebra, boost::numeric::odeint::default_operations, boost::numeric::odeint::initially_resizer>, ErrorChecker=boost::numeric::odeint::default_error_checker<double, boost::numeric::odeint::vector_space_algebra, boost::numeric::odeint::default_operations>, StepAdjuster=boost::numeric::odeint::default_step_adjuster<double, double>, Resizer=boost::numeric::odeint::initially_resizer, System=model, StateIn=state_type]" at line 900
            instantiation of "boost::numeric::odeint::controlled_step_result={boost::numeric::odeint::controlled_step_result} boost::numeric::odeint::controlled_runge_kutta<ErrorStepper, ErrorChecker, StepAdjuster, Resizer, boost::numeric::odeint::explicit_error_stepper_fsal_tag>::try_step_v1(System, StateInOut &, boost::numeric::odeint::controlled_runge_kutta<ErrorStepper, ErrorChecker, StepAdjuster, Resizer, boost::numeric::odeint::explicit_error_stepper_fsal_tag>::time_type &,
                      boost::numeric::odeint::controlled_runge_kutta<ErrorStepper, ErrorChecker, StepAdjuster, Resizer, boost::numeric::odeint::explicit_error_stepper_fsal_tag>::time_type &) [with ErrorStepper=boost::numeric::odeint::runge_kutta_dopri5<double, double, double, double, boost::numeric::odeint::vector_space_algebra, boost::numeric::odeint::default_operations, boost::numeric::odeint::initially_resizer>, ErrorChecker=boost::numeric::odeint::default_error_checker<double,
                      boost::numeric::odeint::vector_space_algebra, boost::numeric::odeint::default_operations>, StepAdjuster=boost::numeric::odeint::default_step_adjuster<double, double>, Resizer=boost::numeric::odeint::initially_resizer, System=model, StateInOut=state_type]" at line 619
            instantiation of "boost::numeric::odeint::controlled_step_result={boost::numeric::odeint::controlled_step_result} boost::numeric::odeint::controlled_runge_kutta<ErrorStepper, ErrorChecker, StepAdjuster, Resizer, boost::numeric::odeint::explicit_error_stepper_fsal_tag>::try_step(System, StateInOut &, boost::numeric::odeint::controlled_runge_kutta<ErrorStepper, ErrorChecker, StepAdjuster, Resizer, boost::numeric::odeint::explicit_error_stepper_fsal_tag>::time_type &,
                      boost::numeric::odeint::controlled_runge_kutta<ErrorStepper, ErrorChecker, StepAdjuster, Resizer, boost::numeric::odeint::explicit_error_stepper_fsal_tag>::time_type &) [with ErrorStepper=boost::numeric::odeint::runge_kutta_dopri5<double, double, double, double, boost::numeric::odeint::vector_space_algebra, boost::numeric::odeint::default_operations, boost::numeric::odeint::initially_resizer>, ErrorChecker=boost::numeric::odeint::default_error_checker<double,
                      boost::numeric::odeint::vector_space_algebra, boost::numeric::odeint::default_operations>, StepAdjuster=boost::numeric::odeint::default_step_adjuster<double, double>, Resizer=boost::numeric::odeint::initially_resizer, System=model, StateInOut=state_type]" at line 103 of "/beegfs/ni86did/libraries/boost_gcc/include/boost/numeric/odeint/integrate/detail/integrate_adaptive.hpp"
            instantiation of "size_t={unsigned long} boost::numeric::odeint::detail::integrate_adaptive(Stepper, System, State &, Time &, Time, Time &, Observer, boost::numeric::odeint::controlled_stepper_tag) [with Stepper=boost::numeric::odeint::controlled_runge_kutta<boost::numeric::odeint::runge_kutta_dopri5<double, double, double, double, boost::numeric::odeint::vector_space_algebra, boost::numeric::odeint::default_operations, boost::numeric::odeint::initially_resizer>,
                      boost::numeric::odeint::default_error_checker<double, boost::numeric::odeint::vector_space_algebra, boost::numeric::odeint::default_operations>, boost::numeric::odeint::default_step_adjuster<double, double>, boost::numeric::odeint::initially_resizer, boost::numeric::odeint::explicit_error_stepper_fsal_tag>, System=model, State=state_type, Time=double, Observer=void (*)(const state_type &, double)]" at line 45 of
                      "/beegfs/ni86did/libraries/boost_gcc/include/boost/numeric/odeint/integrate/integrate_adaptive.hpp"
            instantiation of "size_t={unsigned long} boost::numeric::odeint::integrate_adaptive(Stepper, System, State &, Time, Time, Time, Observer) [with Stepper=boost::numeric::odeint::controlled_runge_kutta<boost::numeric::odeint::runge_kutta_dopri5<double, double, double, double, boost::numeric::odeint::vector_space_algebra, boost::numeric::odeint::default_operations, boost::numeric::odeint::initially_resizer>, boost::numeric::odeint::default_error_checker<double,
                      boost::numeric::odeint::vector_space_algebra, boost::numeric::odeint::default_operations>, boost::numeric::odeint::default_step_adjuster<double, double>, boost::numeric::odeint::initially_resizer, boost::numeric::odeint::explicit_error_stepper_fsal_tag>, System=model, State=state_type, Time=double, Observer=void (*)(const state_type &, double)]" at line 87 of "main.cpp"

Not sure what is going wrong and why my code says it is calling my code with (statetype, double, double) instead of (statetype, statetype, double). What am I doing wrong Syntax-wise? I am compiling with icpc, C++14.

I expected the code to compile properly, based on other forum entries that I saw.

Please trim your code to make it easier to find your problem. Follow these guidelines to create a [minimal reproducible example](https://stackoverflow.com/help/minimal-reproducible-example). — Community, Jun 01 '23 at 17:44

score 1 · Answer 1 · answered Jun 01 '23 at 14:21

You specified the stepper_type as double it will compile if you use

typedef runge_kutta_dopri5< state_type > stepper_type;

I think you should pass m as a reference

integrate_adaptive(make_controlled(1E-12, 1E-12, stepper_type()), 
                   std::ref(m), x, 0.0, T, 1E-12, write_cout);

See this question for more information: Destructor calls during integration with boost odeint

How to correctly use Odeint library for solving ODEs in C++?

1 Answers1