C ++ - 系统的odeint解决方案差异很大
C++ - odeint solution to system varies significantly
我正在使用 odeint 求解一个包含 4 个耦合方程的系统,它模拟车辆在行驶时的振动。我希望我的结果与我在 MATLAB 中得到的结果相似,但不幸的是,这并没有发生。我已经多次检查我的方程式,没有错误,所以问题一定是在积分过程中发生的。
我已经在 MATLAB 中编写了解决方案来验证我从 C++ 脚本中得到了什么。使用相同的条件,这是我从 odeint 获得的解决方案:
这与 MATLAB 中的解决方案相同:
我没想到在 MATLAB 中看到的微振荡会出现在 odeint 结果中,但大多数值甚至都不接近正确。我是否使用了错误的集成器,或者 odeint 是否不适用于此应用程序?
c++文件可以在Github上找到,here.名为"coupledODE"的class是属于我的系统的方程组,odeint是impemented在主要功能。
在 C++ 代码中,您从不执行 calcRadialFreq()
过程,因此 radFreq=0
从其初始化开始就没有变化,提供了一个恒定的漂移项,但没有小的振荡。
将这一行计算合并到上面的 getRoadValues()
过程中,将使结果在视觉上与 Matlab 图完全相同。
假设一个人遵守 20Hz 的强制频率,输出采样率大于 40Hz。对于 100Hz,dt=0.01
的时间步长会很好。
我还建议以小的、易于读取的位来计算 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 );
}
我正在使用 odeint 求解一个包含 4 个耦合方程的系统,它模拟车辆在行驶时的振动。我希望我的结果与我在 MATLAB 中得到的结果相似,但不幸的是,这并没有发生。我已经多次检查我的方程式,没有错误,所以问题一定是在积分过程中发生的。
我已经在 MATLAB 中编写了解决方案来验证我从 C++ 脚本中得到了什么。使用相同的条件,这是我从 odeint 获得的解决方案:
这与 MATLAB 中的解决方案相同:
我没想到在 MATLAB 中看到的微振荡会出现在 odeint 结果中,但大多数值甚至都不接近正确。我是否使用了错误的集成器,或者 odeint 是否不适用于此应用程序?
c++文件可以在Github上找到,here.名为"coupledODE"的class是属于我的系统的方程组,odeint是impemented在主要功能。
在 C++ 代码中,您从不执行 calcRadialFreq()
过程,因此 radFreq=0
从其初始化开始就没有变化,提供了一个恒定的漂移项,但没有小的振荡。
将这一行计算合并到上面的 getRoadValues()
过程中,将使结果在视觉上与 Matlab 图完全相同。
假设一个人遵守 20Hz 的强制频率,输出采样率大于 40Hz。对于 100Hz,dt=0.01
的时间步长会很好。
我还建议以小的、易于读取的位来计算 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 );
}