Libraries/options for calculating second-order ODE in C++ (Visual Studio), based on MATLAB prototype?
<p>Github repo: <a href="https://github.com/armstrongc93/2dof_vibration/" rel="nofollow noreferrer">https://github.com/armstrongc93/2dof_vibration/</a></p>

<p>I have prototyped and validated the math for a 2 degree of freedom vibration analysis in MATLAB, and now need to run it using C++ for an app. I know this is not a trivial task, but it is fairly simple in terms of actual calculations. I'm just not sure of what to look for in a library. </p>

<p>I have looked at Armadillo, odeint and am currently using Eigen for part of the script (calculate natural frequencies). The documentation is a bit hard to follow, so I'm not sure if the answer is staring me in the face. The current MATLAB script I need to work into C++ is as follows (everything except the t and y variables in the anonymous function are defined numerical values, these are user inputs in the C++ script):</p>

<pre><code>% Define parameters &amp; display
m = 1000; % [kg]
J = 800; % [kg*m^2]
k_f = 20000; % [N/m]
k_r = 25000; % [N/m]
l1 = 1.0; % [m]
l2 = 1.5; % [m]
c_f = 1000; % [N*s/m]
c_r = 1000; % [N*s/m]

% Initial conditions and time
IC = [0.1 0 0 0]; % [x x' theta theta']
time = [0 5]; % [s]
fprintf('The input parameters are:\n m =%5d kg\n J =%5d kg*m^2\n',...
 m, J)
fprintf(' kf = %5d N/m\n kr = %5d N/m\n L1 = %3.1f m\n L2 = %3.1f
m\n',...
 k_f, k_r, l1, l2)
fprintf(' cf = %3.1f N*s/m\n cr = %3.1f N*s/m\n x(0) = %3.1f
m\n',c_f,...
 c_r, IC(1))
fprintf(' timespan = 0-%g s\n',time(2))

% Determine radial frequency &amp; y parameters
vel = 10; % [m/s]
L = 0.5; % [m]
A = 0; % [m], =0 for unforced conditions
h = A/2; % [m]
radFreq = (2*pi*vel)/L; % [rad/s]
fprintf('\n')
fprintf('The velocity of the vehicle is: %g m/s\n', vel)

% Solve 4 coupled ODE's
[t,sltn]=ode45(@(t,y)coupledEqtns(t,y,m,J,k_f,k_r,l1,l2,c_f,c_r,...
 radFreq,A,L),time,IC);

% Plot results
plot(t,sltn(:,1),t,sltn(:,3)) % Plot x and theta
title('Displacement and Pitch vs. Time')
xlabel('Time [s]')
ylabel('Amplitude')
legend('Displacement [m]','Pitch [rad]','location','NorthEast')
% Determine maximum values
maximumPitch = max(sltn(:,3));
maximumDisplacement = max(sltn(:,1));

% Display results
fprintf('The distance between bumps is: %3.2f m, and the height is:
%3.3f m\n'...
 ,L,h)
fprintf('The maximum pitch of the vehicle is: %5.4E rad\n',
maximumPitch)
fprintf('The maximum bounce of the vehicle is: %5.4E m\n',...
 maximumDisplacement)
</code></pre>

<p>It is calling the following "coupledEqtns" file:</p>

<pre><code>function firstODE = coupledEqtns(t,y,m,J,kf,kr,L1,L2,cf,cr,radFreq,A,L)

% Function determines first order ODE equations from second order
firstODE = [y(2);
 (1/m)*(-(cf+cr)*y(2)-(kf+kr)*y(1)-(cr*L2-cf*L1)*...
 y(4)-(kr*L2-kf*L1)*y(3)+kf*A*sin((radFreq)*t)+kr*A*sin((radFreq)...
 *t-(2*pi*(L1+L2))/L));
 y(4);
 (1/J)*(-(cf*L1^2+cr*L2^2)*y(4)-(kf*L1^2+kr*L2^2)*y(3)-(cr*L2-cf*...
 L1)*y(2)-(kr*L2-kf*L1)*y(1)+kr*L2*A*sin((radFreq)*t-
(2*pi*(L1+L2))/L)...
 -kf*L1*A*sin((radFreq)*t))];
end
</code></pre>

<p>The anonymous function is definitely what is tripping me up the most here, but this is all the math that is needed to be done. I have not used C++ for anything other than basic operations so this is definitely a bit of a challenge for me. Any and all suggestions are welcome!</p>

<p>The c++ code (from MSVS) is currently as follows for implementation of odeint; I am now getting a subscript error:</p>

<p>Class def:</p>

<pre><code> typedef std::vector&lt; double &gt; state_type;

class coupledODE : public Vehicle, public Road
{
    public:
        void operator()(state_type &amp;x, state_type &amp;dxdt, double t)
        {
            dxdt[0] = x[2];
            dxdt[1] = (1 / mass)*(-(damping_f + damping_r)*x[2] - (stiffness_f + stiffness_r)*x[1] - (damping_r*rearLength - damping_f*frontLength)*x[4] - (stiffness_r*rearLength - stiffness_f*frontLength)*x[3] + stiffness_f*A*sin((radFreq)*t) + stiffness_r*A*sin((radFreq)*t - (2 * pi*(frontLength + rearLength)) / L));
            dxdt[3] = x[4];
            dxdt[4] = (1 / inertia)*(-(damping_f*pow(frontLength, 2) + damping_r*pow(rearLength, 2))*x[4] - (stiffness_f*pow(frontLength, 2) + stiffness_r*pow(rearLength, 2))*x[3] - (stiffness_r*rearLength - stiffness_f*frontLength)*x[2] - (stiffness_r*rearLength - stiffness_f*frontLength)*x[1] + stiffness_r*rearLength*A*sin((radFreq)*t - (2 * pi*(rearLength + frontLength)) / L) - stiffness_f*frontLength*A*sin((radFreq)*t));
        }
};
</code></pre>

<p>Odeint implementation:</p>

<pre><code>    state_type x(4);
    x[0] = car.x_init1;
    x[1] = car.x_init2;
    x[2] = car.x_init3;
    x[4] = car.x_init4;

    const double timeStep = car.dt;
    double tStart = car.t0;
    double t = tStart;
    double tEnd = car.tf;

    // Initialize odeint
    runge_kutta4&lt;state_type&gt;rk4;
    for (size_t i = 0; i &lt; (tStart - tEnd) / timeStep; ++i, t += timeStep)
    {
        rk4.do_step(coupledODE(), x, t, timeStep);
    }
</code></pre>