Search code examples
c++matlabodeodeint

C++ - odeint solution to system varies significantly

I am using odeint to solve for a system of 4 coupled equations, which simulate the vibration of a vehicle while it is driving. I was hoping my results would be similar to what I get in MATLAB, but that is unfortunately not happening. I have checked my equations multiple times and there is no error there, so the issue has to be happening during the integration.

I have coded out the solution in MATLAB to verify what I am getting from the C++ script. Using the same conditions, this is the solution I get from odeint: enter image description here

And this is the same solution in MATLAB: enter image description here

I did not expect the micro-oscillations seen in MATLAB to show up in the odeint results, but most of the values are not even close to correct. Am I using the wrong inegrator, or will odeint just not work for this application?

The c++ file can be found on Github, here. The class named "coupledODE" is the system of equations pertaining to my system, and odeint is impemented in the main function.

Solution

  • In the C++ code you never execute the calcRadialFreq() procedure and thus have radFreq=0 unchanged from its initialization, supplying a constant drift term but not the small oscillations.

    Incorporating this line of computation in the getRoadValues() procedure above it will make the result visually identical to the Matlab graph.

    Provided one honors the forcing frequency of 20Hz with a sampling rate of the output larger than 40Hz. A time step of dt=0.01 for 100Hz will do nicely.

    I'd also propose to compute the ODE function in small, easily to read bits. This should help with comparing different versions and catch errors. It could for instance take the form

    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  );
}