ODE45 Multiple Degrees of Freedom

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

Answer 1

[Steven Lord]

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.

[KostasK]

Thanks for that, this makes things neater.

Answer 2

[KostasK]

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

Answer 3

[mickael dos santos]

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?

[Steven Lord]

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.