How to solve a second order differential equations with matrices by using "ODE45"?

Hi,

I am trying to solve the following 2nd order differential equation...[M]x''+[C]x'+[K]x=fsin(wt) where M,C and K are 3x3 matrices. I am using ODE45 to solve and ideally produce x on y-axis against time on x-axis. I am getting some errors which I believe is due to passing matrices through the code with different dimensions but I am unsure how to reshape them without losing my desired correct result etc.

Any help is much appreciated!

function springmassdamper
tstart = 0;
tend = 60;
n = 10;
tspan = linspace(tstart,tend,n);
xinit = [0;0];

[t,x] = ode45(@integratingfunction, tspan, xinit)

y = x(:,1);
ydot = x(:,2);
plot(tspan,x(:,1))
xlabel('Time')
ylabel('Displacement')

function dxdt = integratingfunction(t,x)
m=[1 0 0;0 1 0;0 0 1];
c=[0.002 -0.001 0;-0.002 0.002 -0.001;0 -0.001 0.001];
k=[2 -1 0;-1 2 -1;0 -1 1];
w=2*pi*50;
f=100;
F=f*sin(w*t);

y=x(1)
ydot=x(2)

dxdt = zeros(size(x));
dxdt(1) = ydot;
dxdt(2) = (1/m)*(F-c*ydot-k*y); 

Answer 1

[James Tursa]

1/m is a scalar/3x3 hence the error. Normally I would advise backslash here, but that brings up another issue. If M, C, and K are all 3x3 matrices, then I would assume x must be a 3x1 vector (at least). So you have a 2nd order equation with a three element vector, making this a 2x3 = 6th order system. Your state vector needs to be 6 elements to handle this, not 2 elements. I.e., the six elements of the 6 element state vector y would be defined as

y(1) = x1
y(2) = x2
y(3) = x3
y(4) = x1'
y(5) = x2'
y(6) = x3'

So your initial conditions must have six states:

xinit = zeros(6,1);

And inside your derivative function would be something like this:

y    = x(1:3);
ydot = x(4:6);
dxdt = zeros(size(x));
dxdt(1:3) = ydot;
dxdt(4:6) = m \ (F - c*ydot - k*y); 

Is it intended that the right hand side, f*sin(wt), is a scalar in this matrix DE?

If x really is supposed to be a scalar, then what are the 3x3 matrices for?

[George]

Hi James,

Thanks for the response, it is very useful.

I can confirm M, C and K are 3x3 matrices and x is a 3x1 vector which consists of acceleration, velocity and displacement like so - [M][x1'';x2'',x3'']+[C][x1';x2',x3']+[K][x1,x2,x3]=[F1,F2,F3] where M,C and K as stated above are 3x3 matrices which define mass, damping and stifness characteristics. The suggested 6 element state vector makes complete sense.

The RHS should be a force matrix consisting of F1, F2 and F3. I am trying to simulate a harmonic forcing function in the form 'Fsin(omega*t) on a 3DOF mass-spring-damper system. This would allow me to pass a range of frequencies using linspace through omega (2*pi*f) using the tspan.

Thanks,

George

[James Tursa]

OK. Right now you have a scalar F, so you will need to change that to have a 3x1 vector F.

[madhan ravi]

+1 James, always learn something new from your answer:)

[Christopher Lum]

George, it sounds like you are trying to solve a very classical dynamic system problem using Mathworks tools. If you are interested, I have a video showing how to do this using Simulink which is effectively a graphical wrapper for ode45 (you can find the video at https://youtu.be/Cvu2zWk3gYw ).

James, great to bump into you in this forum :).

[Satish Jawalageri]

Could I know how to solve [M]{x''}+[C]{x}={F} where [M] and [C] are 3x3 matrix and {x''},{x},{F} are 3x1 matrix and how to plot the same?

Thanks.