C ++ - подходящее решение для системы значительно различается

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

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

И это то же самое решение в MATLAB: введите описание изображения здесь

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

Файл c ++ можно найти на Github, здесь. Класс под названием «coupledODE» представляет собой систему уравнений относящиеся к моей системе, и odeint заблокирован в основной функции.


c++ matlab ode odeint
person cl40    schedule 02.12.2016    source источник
comment
Ссылка на github не работает. А можно ли построить решение Matlab с более высоким разрешением? Может быть, даже построить оба графика на одной диаграмме?   -  person headmyshoulder    schedule 03.12.2016
comment
@headmyshoulder исправлено, прошу прощения.   -  person cl40    schedule 03.12.2016
comment
Не могли бы вы также попытаться интегрировать решение с интеграцией_const (make_dense_output (1E-6, 1E-10, rk5 ()), ...); Просто чтобы убедиться, что это не зависит от точности.   -  person headmyshoulder    schedule 03.12.2016
comment
Ах, отсутствует коэффициент 10   -  person headmyshoulder    schedule 03.12.2016
comment
Это решение еще дальше от желаемого результата; Значит ли это, что это связано с точностью?   -  person cl40    schedule 03.12.2016
comment
Нет я так не думаю   -  person headmyshoulder    schedule 03.12.2016
comment
Вы уверены, что у вас такие же начальные условия и константы?   -  person headmyshoulder    schedule 03.12.2016
comment
Что касается пропущенного коэффициента 10, в данном случае да, он более или менее отклоняется в 10 раз. Когда я расширяю решение до 0-10 секунд, решение не достигает стабильности. Он колеблется сильнее, чем должен, или ведет себя аналогично резонансу, что неверно.   -  person cl40    schedule 03.12.2016
comment
Я могу проверить константы еще раз, но начальные условия просто нулевые, и это отображается в выводе системы из файла cpp.   -  person cl40    schedule 03.12.2016
comment
Если не считать математики, нет очевидной причины, по которой код не работает так, как задумано?   -  person cl40    schedule 03.12.2016
comment
Нет, должно работать. По крайней мере, я был бы очень удивлен, если эта ODE в Matlab будет вести себя иначе.   -  person headmyshoulder    schedule 03.12.2016
comment
Хорошо, это хоть немного обнадеживает. Возможно, порядок операций, выполняемых в уравнениях, некорректен. Я полностью прочесую это (и исходное происхождение) и посмотрю, смогу ли я что-нибудь уловить.   -  person cl40    schedule 03.12.2016
comment
Обязательно выберите модель класса 4 и установите временной шаг несколько меньше, чтобы вы достигли частоты дискретизации для быстрых колебаний, dt=0.01 или dt=0.005 должны подойти. С данными класса 4 из кода C ++ и остальными константами (сравните параметры демпфирования) из файла matlab я могу воспроизвести второй график matlab на python, как с lsoda-scipy.integrate.odeint, так и с моим собственным адаптивным вариантом RK4 - только там не происходит любое изменение размера шага, это dt=0.0033 без изменений для относительной ошибки 1e-5.   -  person Lutz Lehmann    schedule 03.12.2016
comment
Как ни странно, когда я уменьшаю временной шаг, системное решение оказывается совершенно неверным.   -  person cl40    schedule 04.12.2016

Ответы (1)


arrow_upward
1
arrow_downward

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

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

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

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

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  );
}
person Lutz Lehmann    schedule 03.12.2016
comment
Большое спасибо! Я очень ценю это. Это определенно помогает, чтобы на него смотрели другие глаза. - person cl40; 03.12.2016