ode45 for 4 DOF system - How to do with Matrices?

Hello! I am trying to do a ode45 code for a 4 dof vibrational system with 4x4 matrices for mass, damping, and stiffness. The vibration differential equation is Mx''+Cx'+Kx=BF(t). I have my code below for the ode45. I converted x to q to get rid of a coefficient matrix on the second derivative and hope to convert back once q is found using ode45. I am running into an error with using vertcat saying "Dimensions of arrays being concatenated are not consistent" with the error being in the dqdt line. Can anyone possibly help? Thanks!

%Part 3b
q0 = rootM*x0;
qdot0 = rootM*v0;

[t,q]=ode45( @rhs, [0 10], [q0;qdot0] );

x_ode = rootinvM*q;

plot(t,x_ode(:,1));
xlabel('t'); ylabel('x');

function dqdt=rhs(t,q)
M = [4,0,0,0;0,3,0,0;0,0,2.5,0;0,0,0,6];
rootM = sqrt(M);
rootinvM = inv(rootM);
C = [4,-1,0,0;-1,2,-1,0;0,-1,2,-1;0,0,-1,1];
K = [500,-100,0,0;-100,200,-100,0;0,-100,200,-100;0,0,-100,100];
B = [0;0;0;1];
C_tilda = rootinvM*C*rootinvM
K_tilda = rootinvM*K*rootinvM
    dqdt_1 = q(2);
    dqdt_2 = -C_tilda*q(2)-K_tilda*q(1)+rootinvM*B*cos(4*t);
    dqdt=[dqdt_1; dqdt_2];
end

[Star Strider]

I can’t run your code.

What are ‘x0’ and ‘v0’ here:

q0 = rootM*x0;
qdot0 = rootM*v0;

and is ‘rootM’ the same as you calculate in ‘rhs’? (If so, you can pass it as a parrameter, increasing the efficiency of your code).

[Jacey Allen]

I edited the question to upload my whole code, so x0 and v0 are defined above. Hope that helps.

[Indrajit Nandi]

I can't get the above code running. Is there anything I need to add to get this code running? I am a beginner user.

Answer 1

[Jacey Allen]

I have figured out how to run the ode45 with matrices.

[Christopher Smith]

I'm running into the same problem you originally had, what did you do to fix the dqdt line?

[Jacey Allen]

I watched the following youtube video and used the function method. Worked perfectly. Hopefully this helps!

https://www.youtube.com/watch?v=tMpEw3PIp1I

Answer 2

[Steven Lord]

Is there a reason you're computing the inverse of the mass matrix and performing that multiplication instead of using odeset to specify the 'Mass' matrix and solve directly?

[Jacey Allen]

I have never used the ode45 involving matrices and mass matrix so I was not sure how to go about using that.