MATLAB Answers

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

21 views (last 30 days)
George on 24 Feb 2020
Edited: Christopher Lum on 24 Feb 2020
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);
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];
dxdt = zeros(size(x));
dxdt(1) = ydot;
dxdt(2) = (1/m)*(F-c*ydot-k*y);


Sign in to comment.

Answers (1)

James Tursa
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?


Show 1 older comment
James Tursa
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
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 ).
James, great to bump into you in this forum :).

Sign in to comment.

Sign in to answer this question.

Translated by