C++ - решение odeint системы значительно варьируется

Question

Я использую odeint для решения системы из четырех связанных уравнений, которые имитируют вибрацию транспортного средства во время его движения. Я надеялся, что мои результаты будут похожи на то, что я получаю в MATLAB, но, к сожалению, этого не происходит. Я проверил свои уравнения несколько раз, и там нет ошибки, поэтому проблема должна происходить во время интеграции.

Я закодировал решение в MATLAB, чтобы проверить, что я получаю от сценария C++. Используя те же условия, что это решение, которое я получаю от odeint:

И это же решение в MATLAB:

Я не ожидал, что микро-колебаний видели в MATLAB, чтобы показать в результаты odeint, но большая часть значений даже не близка к правильной. Я использую неправильный inegrator, или odeint просто не работает для этого приложения?

Файл C++ можно найти на Github, here. Класс с именем "coupledODE" - это система уравнений, относящихся к моей системе, а odeint - в основной функции.

Answer 1

В коде C++ вы никогда не выполняете процедуру calcRadialFreq() и, следовательно, имеете radFreq=0 без изменений от ее инициализации, обеспечивая постоянный дрейфовый член, но не небольшие колебания.

Включая эту строку вычислений в процедуре getRoadValues() выше, она сделает результат визуально идентичным графику Matlab.

При условии, что одна из них имеет частоту вытеснения 20 Гц с частотой дискретизации выхода, превышающей 40 Гц. Шаг времени dt=0.01 для 100 Гц будет хорошо.

---

Я также хотел бы вычислить функцию ODE в небольших, легко читаемых битах. Это должно помочь в сравнении разных версий и ошибок catch. Он может, например, принять форму

void operator()(state_type &x, state_type &dxdt, double t) 
{ 
    double wave_f = car.stiffness_f*road.A*sin((road.radFreq)*t); 
    double wave_r = car.stiffness_r*road.A*sin((road.radFreq)*t-(2*pi*(car.frontLength+car.rearLength))/road.L); 

    double term1f = car.stiffness_f*x[0] + car.damping_f*x[1]; 
    double term1r = car.stiffness_r*x[0] + car.damping_r*x[1]; 
    double term2f = car.stiffness_f*x[2] + car.damping_f*x[3]; 
    double term2r = car.stiffness_r*x[2] + car.damping_r*x[3]; 
    double term3f = -term1f + term2f*car.frontLength + wave_f; 
    double term3r = -term1r - term2r*car.rearLength + wave_r; 

    dxdt[0] = x[1]; 
    dxdt[1] = (1/car.mass)*(term3f     + term3r    ); 
    dxdt[2] = x[3]; 
    dxdt[3] = (1/car.inertia)*(-term3f*car.frontLength + term3r*car.rearLength ); 
}