MATLAB: ODE45 Multiple Degrees of Freedom

degrees of freedomdifferential equationequation of motionodeode45second ordersystem of equations

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

Best Answer

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.
Related Question