Large error between exact analytic answer and numeric integration

I am trying to solve the equation below in (1) and (2) using the numerical integrator (ode45, ode15s, etc.) within Matlab so I can compare it to the exact analytical steady-state solution from basic vibration analysis. This is a simplified problem to determine the cause of other issues I’m running into with other simulations.

(1) M= [3 1; 1 5],K=[7 -1; -1 10],F1= [1; 4],F2= [-1; 7]

(2 ) M(d^2x/dt^2)+Kx=F1 sin(αt)+F2cos(βt)

According to vibration theory, the exact steady-state analytical solution should be:

(3) x(t)=sin(αt)*(-Mα^2+K)^(-1)*F1+cos(βt)*(-Mβ^2+K)^(-1)*F2

Substituting in the matrices and performing the calculations, the result is:

(4) x(t)= [(-α^2-14)/(14α^4-67α^2+69); (-11α^2+29)/(14α^4-67α^2+69)] sin(αt)+[(12β^2-3)/(14α^4-67α^2+69); -22β^2+48)/(14α^4-67α^2+69))]cos(βt)

In my problem I assume α,β=1, so that the resulting solution is:

(5) x(t)=[0.8125; 1.125]sin(t)+[0.5625; 1.625]cos(t)

Approaching this problem using the numerical integrator (ode45, ode15s, etc.), the problem set up is as follows:

function dy = EoMExample(t,y)

% function to be integrated

%y(1) = x1 dy(1) = x1dot

%y(2) = x2 dy(2) = x2dot

%y(3) = x1dot dy(3) = x1dotdot

%y(4) = x2dot dy(4) = x2dotdot

%y(5) = x1dotdot dy(5) = x1dotdotdot

%y(6) = x2dotdot dy(6) = x2dotdotdot

%y(7) = x1dotdotdot

%y(8) = x2dotdotdot

dy = zeros(8,1);

dy(1) = y(3);

dy(2) = y(4);

dy(3) = (1/3)*(-y(6)-7*y(1)+y(2)+sin(t)-cos(t));

dy(4) = (1/5)* (-y(5)+y(1)-10*y(2)+4*sin(t)+7*cos(t));

dy(5) = y(7);

dy(6) = y(8);

Where the driver is:

%% Driver for integrating equations of motion

runtime = 30; %[sec] Total desired runtime for each

free stream velocity

    %Initial Conditions     
    x10             = 0;        
    x20             = 0;         
    x1dot0          = 0;    
    x2dot0          = 0;      
    x1dotdot0       = 0;        
    x2dotdot0       = 0;         
    x1dotdotdot0    = 0;            
    x2dotdotdot0    = 0;              

    % Integrate EoM for the specified timespan
    tspan1 = [0 runtime];
    x_init = [x10; x20; x1dot0; x2dot0; x1dotdot0; x2dotdot0; x1dotdotdot0; x2dotdotdot0];
    options = odeset('abstol',1e-8,'reltol',1e-5,'InitialStep',1e-4,'MaxStep',1e-3);
    [t, x] = ode45(@(t,y) EoMExample(t,y), tspan1, x_init, options);

If I run this driver and plot the results for x1 and x2, it is very chaotic and does not reach a steady-state value. The two states are not simple sinusoids as represented by the analytic solution (eqn. (5)). My question is how there could be such a major discrepancy? I have tried running the integrator for large values to see if steady-state is achieved, but it never happens. I’ve tried stiff and non-stiff integrators with the same result from both. Is there some fundamental aspect of this problem I have overlooked? Are the numerical integrators within Matlab simply not able to handle coupled 2nd order differential equations? Any and all feedback is much appreciated.

Answer 1

[Mark Shore]

I rarely work with differential equations but you should check your basic equation again.

For some reason you have a 1D harmonic oscillator with NO damping term and two sinusoidal driving forces. Normally there should be a dx/dt term (even a small one) in your equation.

Without a damping term there is nothing to stop the amplitudes from growing, generating numerical instability and never reaching a steady state.

[Thomas]

I have purposely neglected any damping so I can try to pinpoint the cause of the error. Even without damping the system should reach a steady-state oscillation which should be identical to the analytical solution of the superposition of a sine and cosine of the same frequency.

Answer 2

[Mike Hosea]

Your formulation has many problems. I don't see that you're using a mass matrix with ODE45, and you also have not inverted M, at least not correctly. Why 8 components? You're setting the 2nd and 3rd derivatives to zero. The steady state solution does not have that property, so obviously this formulation is incorrect. When I substitute a correct formulation of the problem (with only 4 components) and begin the integration on the steady-state solution, it stays there.