# ODE45 Multiple Degrees of Freedom

조회 수: 90(최근 30일)
Kostas 29 Jul 2019
Commented: Kostas 30 Jul 2019
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

로그인 to comment.

### 채택된 답변

Steven Lord 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.

#### 댓글 수: 1

Kostas 30 Jul 2019
Thanks for that, this makes things neater.

로그인 to comment.

Kostas 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

#### 댓글 수: 0

로그인 to comment.

이 질문에 답변하려면 로그인을(를) 수행하십시오.