Running the code
m_s = 5; % Mass of the spring
m_r = 68; % Mass of the rotor
k_r = 5.45e5; % Spring constant of the rotor
c_s = 700; % Damping coefficient of the spring
c_r = 900; % Damping coefficient of the rotor
k_s = 1.5e4; % Spring constant of the spring
% Frequency range in rad/s
omega = 0:160;
% Initialize arrays to store roots
roots_real = zeros(length(omega), 5);
roots_imag = zeros(length(omega), 5);
% Initialize array to store crossing points
crossing_omega = zeros(1, 5);
crossing_from_real_to_imag = zeros(1, 5);
% Calculate and store roots for each angular frequency
for i = 1:length(omega)
a4_omega = c_s./m_s + k_r./c_r - 1i.*omega(i);
a3_omega = k_r./m_r + k_r./m_s + k_s./m_s + c_s./m_s*(k_r./c_r - 1i.*omega(i));
a2_omega = (k_r.*c_s) ./ (m_r.*m_s) + (k_r.*k_s) ./ (m_s.*c_r) - (k_r./m_r + k_r./m_s + k_r./m_s).*1i.*omega(i);
a1_omega = k_r ./m_r * (k_s./m_s - 1i.*omega(i).*c_s./m_s);
a0_omega = -1i.*omega(i).*k_s.*k_r ./(m_s.*m_r);
a = [1, a4_omega, a3_omega, a2_omega, a1_omega, a0_omega];
% Calculate roots
system_roots = roots(a);
% Store real and imaginary parts of the roots
roots_real(i, :) = real(system_roots);
roots_imag(i, :) = imag(system_roots);
end
%%
figure
plot(roots_real(:,1), roots_imag(:,1), 'ro')
hold all
plot(roots_real(:,2), roots_imag(:,2), 'bo')
plot(roots_real(:,3), roots_imag(:,3), 'mx')
plot(roots_real(:,4), roots_imag(:,4), 'k*')
plot(roots_real(:,5), roots_imag(:,5), 'kd')
ylim([-20 50])
xlim([-16 6])
ylim([-1500 1e4])
xlim([-1050 200])
xlabel('Im(s)')
ylabel('Re(s)')
figure
plot(roots_real(:,1), roots_imag(:,1), 'ro')
hold all
plot(roots_real(:,2), roots_imag(:,2), 'bo')
plot(roots_real(:,3), roots_imag(:,3), 'mx')
plot(roots_real(:,4), roots_imag(:,4), 'k*')
plot(roots_real(:,5), roots_imag(:,5), 'kd')
ylim([-20 150]) % changed the limits
xlim([-16 6])
xlabel('Im(s)')
ylabel('Re(s)')
grid
One of the challenge in a problem like this is that the order how the roots are returned from roots is not specified. Hence, when we keep appending them to the roots_ arrays for each value of omega, there's no guarantee that the columns the roots_ arrays will form continuous loci branches. We see this effect, for example, with the lone black diamond at Re(s) = -7.5 or so (btw, the axes are mislabeled). Fortunately, for this particular problem, we see the third and fifth columns of roots_ are continuous at their crossings of the jw axis.
Plot the real part of the fifth root (can do the same for the third).
figure
plot(omega,roots_real(:,5)),grid,xlabel('omega')
You could just get an answer by zooming in on the plot. To be more precise, try using fzero asusming that the roots_ can be reasonably approximated with linear interpolation between values of omega (feel free to try other options and/or use finer spacing of omega).
omegacross = fzero(@(w) interp1(omega,roots_real(:,5),w),25)
The value of the imaginary part at that frequency is
interp1(omega,roots_imag(:,5),omegacross)
