Eigenvalue vibrational response Two DOF spring- mass system

A 2 degree of freedom (2 DOF) system with 3 springs only is given. I am little bit confused to write the codes. I would be glad for some help.

[ 1 0 𝒙̈(𝑡) + [ 12 −2 x(𝑡) = 0 and 𝒙̇(0) = 0 , 𝒙(0) = [1 1]^𝑇

0 4] 2 12]

So here, m1=1 and m2=4 while k1=k3=10 and k2=2. But I do not know why my code is wrong.

K1=10 ; K2=2 ; K3=10 ;
M1=1   ; M2=4 ;
M= [ M1 0 ;
     0  M2] ;
K=[K1+K2  -K2;
    -K2   K2+K3] ;
[modeShape fr]=eig (K,M); %estimation of natural frequenciesand mode shapes
A00=zeros(2); A11=eye(2);
CC=[A00 A11;-inv(M)*K A00];
global CC
max_freq=max(sqrt(diag(fr))/(2*pi)); %highest frequency in Hz
dt=1/(max_freq*20);
time=0:dt:200*dt;
y0=[1 1 0 0]; %[displ1 disp2 vel1 vel2] initial condition
plot(time,ysol(:,1:2),'linewidth',2)
xlabel ('Time')
ylabel('displacement')
ylim([-.02 .02])
grid on

[Ali]

How did you know your code is wrong? Are you getting an error message or is the plot different than what you expect?

[Joni Riihimaki]

It shows error in my y0 line and does not plot what I should get

Answer 1

[Gyan Vaibhav]

Hi Joni,

I understand that you are trying to solve a spring mass system with two degrees of freedom with 3 springs and 2 masses.

The code seems to have a good start on setting up the system for simulation, but there are a few issues and missing steps in the code that need to be addressed.

Here are a few points to be rectified:

1. To find solutions with the conditions, a solver needs to be used. “ode45” can be used for this purpose.
2. The outputs need to be stored and then the correct variable needs to be plotted. In the above code, “ysol” has been used to plot; however it was never defined, resulting in an error.

Here is the code with some corrections and improvements, to give the desired output, and better readability:

% Define system parameters 
K1 = 10; K2 = 2; K3 = 10; 
M1 = 1; M2 = 4; 

% Mass matrix 
M = [M1 0; 0 M2]; 

% Stiffness matrix 
K = [K1+K2 -K2; -K2 K2+K3]; 

% Calculate eigenvalues and eigenvectors (mode shapes) 
[modeShape, fr] = eig(K, M); 

% Calculate natural frequencies in Hz 
nat_freqs = sqrt(diag(fr)) / (2 * pi); 

% Highest frequency in Hz 
max_freq = max(nat_freqs); 

dt = 1 / (max_freq * 20); 

% Simulation time span 
time_span = 0:dt:200*dt;  

% Initial conditions [x1(0) x2(0) x1_dot(0) x2_dot(0)] 
y0 = [1 1 0 0]; 

% State-space representation 
A = [zeros(2), eye(2); -M\K, zeros(2)]; 

% ODE function 
ode_func = @(t, y) A * y; 

% Solve ODE 
[time, ysol] = ode45(ode_func, time_span, y0); 

% Plot the results 
plot(time, ysol(:, 1:2), 'linewidth', 2);  
xlabel('Time (s)'); 
ylabel('Displacement (m)'); 
ylim([-1.2 1.2]); % Adjust the y-axis limits if needed 
grid on; 
title('Displacement vs. Time for a 2 DOF System'); 
legend('Mass 1 Displacement', 'Mass 2 Displacement'); 

Please note the following corrections and additions:

1. The ODE function “ode_func” is defined as a function handle.
2. The ode45 solver is used to simulate the system response over the required time span.
3. The plot function is corrected to plot the displacements.

Refer to the following documentation to learn more about “ode45” solver.

https://www.mathworks.com/help/matlab/ref/ode45.html

Hope this helps.

Thanks

Gyan