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

45 次查看(过去 30 天)
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);

回答(1 个)

James Tursa
James Tursa 2020-2-24
编辑:James Tursa 2020-2-24
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?
  5 个评论
Christopher Lum
Christopher Lum 2020-2-24
编辑:Christopher Lum 2020-2-24
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
编辑:Satish Jawalageri 2020-6-2
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.

请先登录,再进行评论。

类别

Help CenterFile Exchange 中查找有关 Programming 的更多信息

产品

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by