如何正确使用Odeint库在C++中求解常微分方程？

mwg9r5ms  于 2023-06-07  发布在  其他
关注(0)|答案(1)|浏览(212)

我正在尝试为物理模拟集成odeint。我建立了一个模型，包含我的方程的所有参数：
（型号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>);
};

这是我的模型。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;
    }
    */
}

这是我的主页面。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;
}

不知何故，这并没有正确编译。错误很长。下面是前几行：

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"

不知道哪里出了问题，为什么我的代码说它是用（statetype，double，double）而不是（statetype，statetype，double）调用我的代码。我在语法方面做错了什么？我用icpc编译，C++14。
根据我看到的其他论坛条目，我期望代码能够正确编译。

c++

1条答案

按热度按时间
l2osamch

l2osamch1#

您将stepper_type指定为double，如果使用

typedef runge_kutta_dopri5< state_type > stepper_type;

我认为你应该传递m作为参考

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

有关详细信息，请参阅此问题：Destructor calls during integration with boost odeint

赞(0）分享  回复(0）举报  2023-06-07
我来回答

相关问题

    查看更多

    热门标签

    更多
    JavaquerypythonNode开发语言requestUtil数据库Table后端算法LoggerMessageElementParser

    最新问答

    更多
    • 蜀ICP备13028337号-1 大数据知识库 https://www.saoniuhuo.com © All rights reserved
    • 本站内容来源互联网,如果侵犯您的权益请联系我们删除, 联系方式：448109455@qq.com