ODE45 Multiple Degrees of Freedom

69 views (last 30 days)
KostasK
KostasK on 29 Jul 2019
Commented: Steven Lord on 24 Nov 2020
Hi all,
I am having difficulty in modelling a 3DOF system usng ODE45 as I am not getting the correct result. Therefore I would like to ask how is it possible to model the below problem?
First of all, here is the problem. This is a typical equation of motion in matrix form, with no exitation force. So the objective is to find the displacement and velocity of the system for a time of 0 to 60 seconds. The data given is m1=m2=m3=1kg and k1=k2=k3=25N/m, and the initial conditions is that when the displacement of all carts is 0m, the velocity should be 1m/s for all.
Naturally, the above creates 3 equations of motion, and here is the code I have created below. I have been unable to find an example with a system of 3 second order ODEs, so I am suspecting that I am doing something wrong with the syntax of the code in the 'odefcn' part:
% Inputs
% Masses kg
m1 = 1 ;
m2 = 1 ;
m3 = 1 ;
% Spring coefficients N/m
k1 = 25 ;
k2 = 25 ;
k3 = 25 ;
% Matrices
% Mass
M = diag([m1 m2 m3]);
% Spring
K1 = diag([k1 + k2, k2 + k3, k3]) ;
K2 = diag([-k2, -k3], 1) ;
K3 = diag([-k2, -k3], -1) ;
K = K1 + K2 + K3 ;
% ODE Solution
% Initial Conditions
tspan = [0 60] ;
y0 = [1 1 1 0 0 0] ;
% Solution
[t, x] = ode45(@(t, x) odefcn(t, x, M, K), tspan, y0) ;
% Results
x_ = x(:, 4:end) ;
xdot_ = x(:, 1:3) ;
% Plot
figure
plot(t, x_)
grid on
xlabel('Time (s)')
% ODE Function
function dxdt = odefcn(t, x, M, K)
dxdt = zeros(6, 1) ;
dxdt(1) = x(1) ;
dxdt(2) = x(2) ;
dxdt(3) = x(3) ;
dxdt(4) = -K(1)/M(1) * x(4) - K(2)/M(1) * x(5) ;
dxdt(5) = -K(4)/M(5) * x(4) - K(5)/M(5) * x(5) - K(6)/M(5) * x(6) ;
dxdt(6) = -K(8)/M(9) * x(5) - K(8)/M(9) * x(6) ;
end

Accepted Answer

Steven Lord
Steven Lord on 29 Jul 2019
Rather than passing your mass matrix into your function, I recommend creating an options structure using odeset. Specify the Mass option in that options structure then pass it into the ODE solver. That way all your ODE function needs to do is compute -k*x. See the Solve Stiff Differential Algebraic Equation example on the ode23t documentation page for a demonstration of how to set up the options structure and call the solver.

More Answers (2)

KostasK
KostasK on 29 Jul 2019
I have found the problem, and as I suspected, it was in the function. the correct function is:
(i still don't fully understand why)
% ODE Function
function dxdt = odefcn(t, x, M, K)
dxdt = zeros(6, 1) ;
dxdt(1) = x(4) ;
dxdt(2) = x(5) ;
dxdt(3) = x(6) ;
dxdt(4) = -K(1)/M(1) * x(1) - K(2)/M(1) * x(2) ;
dxdt(5) = -K(4)/M(5) * x(1) - K(5)/M(5) * x(2) - K(6)/M(5) * x(3) ;
dxdt(6) = -K(8)/M(9) * x(2) - K(9)/M(9) * x(3) ;
end

mickael dos santos
mickael dos santos on 24 Nov 2020
hello,
i don't really understand your paragraphe about odefcn. i would like to apply your code with my matrix equation could you help me?
  1 Comment
Steven Lord
Steven Lord on 24 Nov 2020
See this documentation example of a system of ODEs that includes a mass matrix. You may be able to use it as a model for how to solve your system of ODEs.

Sign in to comment.

Community Treasure Hunt

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

Start Hunting!