Using ode45 gives exponentially increasing solution

I am trying to solve a 6 degree of freedom equation [Mhat]x'' + [Bhat]x' + [K]x = [F]sin(wt). [Mhat],[Bhat] and [K] are all 6x6 matrices. This is the ODE function I am using.

function dxdt=MDOF1(t,x)
Bhat = [0.000223552500000000,6.94950000000000e-08,5.13525000000000e-07,-4.55100000000000e-09,3.75150000000000e-05,-2.68755000000000e-09;-1.21770000000000e-07,0.00161437500000000,1.23922500000000e-05,-4.76625000000000e-05,3.22875000000000e-08,1.37145000000000e-07;4.73550000000000e-06,3.75150000000000e-06,0.719550000000000,-8.54850000000000e-07,-8.17950000000000e-06,6.42675000000000e-07;9.47100000000000e-09,-4.88925000000000e-05,-6.94950000000000e-07,1.42372500000000e-06,-1.62975000000000e-09,-8.02575000000000e-09;3.72075000000000e-05,1.97107500000000e-07,-5.28900000000000e-06,-9.25575000000000e-09,6.30375000000000e-06,3.19800000000000e-09;3.87450000000000e-09,2.59530000000000e-07,1.83885000000000e-06,-2.82900000000000e-08,1.86037500000000e-08,7.77975000000000e-08];
Mhat = [1024.33244443821,0.144525000000000,-0.141450000000000,-0.0370025000000000,-36.2850000000000,0.0145550000000000;-0.149650000000000,2745.30744443821,-0.718525000000000,-88.3550000000000,0.0192700000000000,0.192700000000000;-0.474575000000000,-0.333125000000000,1501.98244443821,0.0795400000000000,0.745175000000000,-0.0570925000000000;0.0267525000000000,-89.0725000000000,0.0799500000000000,71.8652857627223,-0.00697000000000000,-0.00697000000000000;-36.2850000000000,0.0892775000000000,0.632425000000000,0.0131200000000000,1839.42588027002,0.0879450000000000;-0.0720575000000000,0.248050000000000,-0.0427425000000000,-0.0329025000000000,0.262400000000000,2124.90165940969];
K = [0,0,0,0,0,0;0,0,0,0,0,0;0,0,6005.09482852500,-0.000341104906106250,0.000290203027676250,0;0,0,-0.000341104906106250,-75.8401802437500,1.98933466826250e-05,-0.000500836689307500;0,0,0.000290203027676250,1.98933466826250e-05,2134.38089767500,2.51222481753750e-05;0,0,0,0,0,0];
F = [0.00749847190933425 + 146.863981327050i;0.100850251293675 + 0.0135606262961325i;5871.33767882513 + 0.0419397503470725i;-0.0121112998456613 - 0.00142823795502375i;-0.0437293876182638 + 24.4266695838750i;0.00837727213389788 - 0.0102890246599238i];
freq = 0.3*ones(6,1);
absF = abs(F);

FEXT = abs(F).*(sin(freq.*t-angle(F)));
dxdt = ones(12,1);
dxdt(1:6)=x(7:12);
dxdt(7:12)= inv(Mhat)*FEXT -inv(Mhat)*K*x(1:6)-inv(Mhat)*Bhat*x(7:12);
end

This is the main program to run ode45

    t=[0 300];
    IC = zeros(12,1);
    opts = odeset('RelTol',1e-3,'AbsTol',1e-4);
    [t,x]=ode45( @MDOF1,t,IC,opts);

The issue is, the analytical solution for all 6 degrees is sinusoidal (attached for x(1))

but when I run the solver, the output is an exponentially increasing function

I'm not sure why the error is occuring. I thought it was a time step error and ran it with a maximum step size of 0.0001 s. Still the same error is occuring.

I used these same matrices to build a model on Simulink and the error carried there as well.

[Walter Roberson]

Note: in theory you can get more precision with

dxdt(7:12)= Mhat\(FEXT - K*x(1:6) - Bhat*x(7:12));

This does not, however, fix the problem that you noted.

[David Goodmanson]

Hi Aghamarshana,

the K matrix is not symmetric, which looks suspicious. However, that is probably not the real issue. The main problem could well be the damping term:

eig(-inv(Mhat)*Bhat)
ans =
  -0.00047907
  -2.2323e-07
  -5.8818e-07
   2.8822e-10
  -2.9585e-11
  -4.0592e-11

As you can see, the fourth eigenvalue is positive, which leads to exponential growth.

[Aghamarshana Meduri]

Hi David,

Thank you for identifying the issue. I have solved the problem by setting a threshold of 10^-4, less than which all values are automatically set to zero. That seems to have done the trick.

Thanks again.