21 views (last 30 days)

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);

James Tursa
on 24 Feb 2020

Edited: James Tursa
on 24 Feb 2020

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?

James Tursa
on 24 Feb 2020

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

Christopher Lum
on 24 Feb 2020

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 :).

Sign in to comment.

Sign in to answer this question.

Opportunities for recent engineering grads.

Apply Today
## 0 Comments

Sign in to comment.