real part of eigenvalues vs damping plot

Hello everyone, please i wanna ask about plotting real part of eigenvalues vs damping while varying the support stiffness (k1xx); here my code :

clear all
close all
mr=306, mj=141.5, mb=7, ks=4.9*10^7, cs=17.51, ci=0, N=10000, q= 3.5*10^6;

mr = 306

mj = 141.5000

mb = 7

ks = 49000000

cs = 17.5100

ci = 0

N = 10000

kxx=2.25*10^8, kyy=2.5*10^8, cxx=2.1*10^5, cyy=2.25*10^5,kxy=0,kyx=0,cxy=0,cyx=0;

kxx = 225000000

kyy = 250000000

cxx = 210000

cyy = 225000

kxy = 0

kyx = 0

cxy = 0

k1xx=1.75*10^7, c1xx=10^6; k1yy=1.75*10^7, k1xy=0, k1yx=0,c1yy=10^6,c1xy=0,c1yx=0;

k1xx = 17500000

k1yy = 17500000

k1xy = 0

k1yx = 0

c1yy = 1000000

c1xy = 0

w=sqrt(ks/mj);
M1=[mr 0 0;0 mr 0;0 0 mj];
M2=[mj 0 0;0 mb 0;0 0 mb];
Mtot=[M1,zeros(3,3);zeros(3,3),M2];
C1=[ci+cs 0 -ci;0 ci+cs 0;-ci/2 0 (ci/2)+cxx];
C2=[0 0 0;-ci 0 0;cxy -cxx -cxy];
C3=[0 -ci/2 cyx;0 0 -cxx;0 0 -cyx];
C4=[ci/2+cyy -cyx -cyy;-cxy cxx+c1xx cxy+c1xy;-cyy cyx+c1yx cyy+c1yy];
Ctot=[C1 C2;C3 C4];
K1=[ks (w*ci+q) -ks;-(w*ci+q) ks w*ci;-ks/2 -w*ci/2 (ks/2+kxx)];
K2=[-w*ci 0 0;-ks 0 0;(kxy+w*ci/2) -kxx -kxy];
K3=[w*ci/2 -ks/2 (kyx-w*ci/2);0 0 -kxx;0 0 -kyx];
K4=[(ks/2+kyy) -kyx -kyy;-kxy (kxx+k1xx) (kxy+k1xy);-kyy (kyx+k1yx) (kyy+k1yy)];
Ktot=[K1 K2;K3 K4];
A=[zeros(6,6),Mtot;Mtot,Ctot];
B=[-Mtot,zeros(6,6);zeros(6,6),Ktot];
[V,D]=eig(A,B);
lamda=diag(D);
omega=sqrt(lamda);
k1xx=[8.75*10^6 1.75*10^7 4.37*10^7 8.75*107];
c1xx=10^5:10^7;
for i=1:length(k1xx) 
    [V,D]=eig(A,B);
    lamda=diag(D);
end 
plot(c1xx,real(lamda))

Error using plot
Specify the coordinates as vectors or matrices of the same size, or as a vector and a matrix that share the same length in at least one dimension.

figure(1)
hold on 
title('stabilité du systeme')
xlabel('Amortissement du support')
ylabel ('Partie reelle des valeurs propres')

and i get this error message: share the same length in at least one dimension.

Error in amortissement_support (line 32)

plot(c1xx,real(lamda))

[Umar]

Hi @Amira, Please let us know if you need any further assistance.

Answer 1

[Torsten]

I don't know what you try to plot, but the cell matrix lamda below should contain all information you need.

mr=306, mj=141.5, mb=7, ks=4.9*10^7, cs=17.51, ci=0, N=10000, q= 3.5*10^6;
kxx=2.25*10^8, kyy=2.5*10^8, cxx=2.1*10^5, cyy=2.25*10^5,kxy=0,kyx=0,cxy=0,cyx=0;
%k1xx=1.75*10^7, c1xx=10^6; 
k1yy=1.75*10^7, k1xy=0, k1yx=0,c1yy=10^6,c1xy=0,c1yx=0;
w=sqrt(ks/mj);
M1=[mr 0 0;0 mr 0;0 0 mj];
M2=[mj 0 0;0 mb 0;0 0 mb];
Mtot=[M1,zeros(3,3);zeros(3,3),M2];
C1=[ci+cs 0 -ci;0 ci+cs 0;-ci/2 0 (ci/2)+cxx];
C2=[0 0 0;-ci 0 0;cxy -cxx -cxy];
C3=[0 -ci/2 cyx;0 0 -cxx;0 0 -cyx];
K1=[ks (w*ci+q) -ks;-(w*ci+q) ks w*ci;-ks/2 -w*ci/2 (ks/2+kxx)];
K2=[-w*ci 0 0;-ks 0 0;(kxy+w*ci/2) -kxx -kxy];
K3=[w*ci/2 -ks/2 (kyx-w*ci/2);0 0 -kxx;0 0 -kyx];
K1xx=[8.75*10^6 1.75*10^7 4.37*10^7 8.75*10^7];
C1xx=10^5:10^7;
for iK = 1:numel(K1xx)
    k1xx = K1xx(iK);
    for iC = 1:numel(C1xx)
      c1xx = C1xx(iC);
      K4=[(ks/2+kyy) -kyx -kyy;-kxy (kxx+k1xx) (kxy+k1xy);-kyy (kyx+k1yx) (kyy+k1yy)];
      Ktot=[K1 K2;K3 K4];
      C4=[ci/2+cyy -cyx -cyy;-cxy cxx+c1xx cxy+c1xy;-cyy cyx+c1yx cyy+c1yy];
      Ctot=[C1 C2;C3 C4];
      A=[zeros(6,6),Mtot;Mtot,Ctot];
      B=[-Mtot,zeros(6,6);zeros(6,6),Ktot];
      [V,D]=eig(A,B);
      lamda{iK,iC}=real(diag(D));           
    end
end

[Walter Roberson]

Note that your original code had

for i=1:length(k1xx) 
    [V,D]=eig(A,B);
    lamda=diag(D);
end 

which overwrites all of lambda each time through the loop.

And note that the length of the diagonal is unlikely to be the same as the length of k1xx

Answer 2

[Umar]

Hi @Amira,

I reviewed your code and did some research by going through mathworks documentations to find the root cause of your dimension mismatch error. Your original loop structure was overwriting the lambda values at each iteration without storing them, so when you tried to plot c1xx (a vector with many values) against lambda (only containing values from the final iteration), MATLAB couldn't match the dimensions. Beyond this loop issue, your eigenvalue formulation using eig(A,B) with those specific A and B matrices was incorrect for the standard mechanical system equation Mx_ddot + Cx_dot + Kx = 0. For rotor dynamics problems with gyroscopic effects like yours, MATLAB's polyeig function is the correct approach as it directly solves the quadratic eigenvalue problem (Mlambda^2 + C*lambda + K)*x = 0. I've corrected your code to properly store eigenvalues from each damping value and implemented both the polyeig method and state-space formulation for comparison.

Your results are now showing physically correct behavior. The analysis reveals that both methods produce identical eigenvalues, which validates the mathematical formulation. At low damping (c1xx = 1.00e+05), your system is completely stable with zero unstable modes. However, as damping increases to medium levels (c1xx = 4.95e+06), two modes become unstable with a maximum real part of 1.29, and at high damping (c1xx = 1.00e+07), these two modes remain unstable with a maximum real part of 2.18. This demonstrates a critical damping threshold at approximately c1xx = 2.73e+06 N·s/m where the system transitions from stable to unstable operation.

Please see attached figures, Figure 1 shows the full eigenvalue spectrum using the polyeig method across all damping values, with most modes remaining deeply stable (large negative real parts) while a green line clearly shows one mode becoming progressively unstable as damping increases. Figure 2 displays the state-space method results showing identical behavior, confirming both approaches are mathematically equivalent. Figure 3 provides a detailed comparison of the top four most significant modes from both methods in side-by-side subplots, where you can see Mode 4 (purple line) starts highly unstable at low damping but rapidly stabilizes, while Modes 1, 2, and 3 hover near the stability boundary (the red dashed line at zero), with both panels showing perfect agreement between the two solution methods.

This counterintuitive phenomenon where increasing damping destabilizes the system is not an error but rather a well-documented effect in rotor dynamics called destabilizing damping. Your system has significant gyroscopic coupling evidenced by the rotor speed parameter q = 3.50e+06 and a stiffness symmetry index of 3.52 percent, which is normal and expected for rotating machinery. The gyroscopic forces create cross-coupled stiffness terms that interact with damping in a complex way. Below the critical threshold, damping acts conventionally to stabilize the system, but above this threshold, the interaction between damping forces and gyroscopic coupling creates tangential forces that pump energy into certain whirl modes rather than dissipating it. This occurs because in rotating systems operating at supercritical speeds or with strong cross-coupling, damping forces can lead the displacement and create destabilizing effects on forward whirl modes. The system properties show your natural frequency is 588.46 rad/s (approximately 94 Hz), and the mass and stiffness matrices are well-conditioned with condition numbers around 40-80, indicating the numerical solution is reliable.

For your analysis and any technical report, you should use the polyeig results as this is the industry standard method specifically designed for quadratic eigenvalue problems in rotor-bearing systems. Your findings demonstrate sophisticated rotor dynamics behavior that would be excellent to discuss in academic work. You can state that the stability analysis reveals a critical damping threshold at approximately 2.73e+06 N·s/m, above which gyroscopic coupling induces destabilizing forces, and that this counterintuitive behavior where increased damping reduces stability is characteristic of rotor systems with significant cross-coupled stiffness terms. For practical design purposes, you would want to maintain support damping below this critical value, likely in the range of 1-2e+06 N·s/m, to ensure stable operation across all operating conditions.

[Umar]

@Amira,please see attached since I was getting errors when trying to type corrected code over here. Also, I would like to let you know that @Torsten and @Walter Robertson are experts and have years of technical experience. I would definitely rely on their suggestions shared with you.

Good luck

clear all
close all
% System parameters
mr = 306;
mj = 141.5;
mb = 7;
ks = 4.9*10^7;
cs = 17.51;
ci = 0;
N = 10000;
q = 3.5*10^6;
kxx = 2.25*10^8;
kyy = 2.5*10^8;
cxx = 2.1*10^5;
cyy = 2.25*10^5;
kxy = 0;
kyx = 0;
cxy = 0;
cyx = 0;
k1yy = 1.75*10^7;
k1xy = 0;
k1yx = 0;
c1xy = 0;
c1yx = 0;
% Fixed support stiffness
k1xx = 1.75*10^7;
% Varying damping coefficient
c1xx_values = linspace(10^5, 10^7, 50);
c1yy = 10^6;
w = sqrt(ks/mj);
% Mass matrices
M1 = [mr 0 0; 0 mr 0; 0 0 mj];
M2 = [mj 0 0; 0 mb 0; 0 0 mb];
M = [M1, zeros(3,3); zeros(3,3), M2];
% Stiffness matrices (constant parts)
K1 = [ks (w*ci+q) -ks; -(w*ci+q) ks w*ci; -ks/2 -w*ci/2 (ks/2+kxx)];
K2 = [-w*ci 0 0; -ks 0 0; (kxy+w*ci/2) -kxx -kxy];
K3 = [w*ci/2 -ks/2 (kyx-w*ci/2); 0 0 -kxx; 0 0 -kyx];
K4 = [(ks/2+kyy) -kyx -kyy; -kxy (kxx+k1xx) (kxy+k1xy); -kyy (kyx+k1yx) (kyy+k1yy)];
K = [K1 K2; K3 K4];
% Damping matrices (constant parts)
C1 = [ci+cs 0 -ci; 0 ci+cs 0; -ci/2 0 (ci/2)+cxx];
C2 = [0 0 0; -ci 0 0; cxy -cxx -cxy];
C3 = [0 -ci/2 cyx; 0 0 -cxx; 0 0 -cyx];
% Storage for eigenvalues from both methods
num_dof = 6;
num_eigenvalues = 2 * num_dof;
real_parts_polyeig = zeros(length(c1xx_values), num_eigenvalues);
real_parts_statespace = zeros(length(c1xx_values), num_eigenvalues);
fprintf('=== ROTOR DYNAMICS EIGENVALUE ANALYSIS ===\n');

=== ROTOR DYNAMICS EIGENVALUE ANALYSIS ===

fprintf('Using BOTH polyeig (QEP) and state-space methods\n\n');

Using BOTH polyeig (QEP) and state-space methods

% Loop over damping values
for i = 1:length(c1xx_values)
c1xx = c1xx_values(i);
% Update damping matrix
C4 = [ci/2+cyy -cyx -cyy; -cxy cxx+c1xx cxy+c1xy; -cyy cyx+c1yx cyy+c1yy];
C = [C1 C2; C3 C4];
% METHOD 1: MATLAB's polyeig (RECOMMENDED for rotor dynamics)
% Solves: M*lambda^2 + C*lambda + K = 0
[~, lambda_polyeig] = polyeig(K, C, M);
[~, sort_idx] = sort(real(lambda_polyeig), 'descend');
real_parts_polyeig(i, :) = real(lambda_polyeig(sort_idx));
% METHOD 2: Standard state-space formulation (for comparison)
A_state = [zeros(num_dof), eye(num_dof); -M\K, -M\C];
lambda_state = eig(A_state);
[~, sort_idx] = sort(real(lambda_state), 'descend');
real_parts_statespace(i, :) = real(lambda_state(sort_idx));
end
fprintf('Analysis complete!\n\n');

Analysis complete!

% Plot Method 1: polyeig results
figure(1)
plot(c1xx_values, real_parts_polyeig, 'LineWidth', 1.5)
grid on
title('Stabilité du système - Méthode polyeig (QEP) - RECOMMANDÉE')
xlabel('Amortissement du support c_{1xx} (N·s/m)')
ylabel('Partie réelle des valeurs propres (1/s)')
hold on
yline(0, 'r--', 'LineWidth', 2, 'DisplayName', 'Limite de stabilité')
hold off
legend('Location', 'best')

% Plot Method 2: state-space results (for comparison)
figure(2)
plot(c1xx_values, real_parts_statespace, 'LineWidth', 1.5)
grid on
title('Stabilité du système - Méthode State-Space (comparaison)')
xlabel('Amortissement du support c_{1xx} (N·s/m)')
ylabel('Partie réelle des valeurs propres (1/s)')
hold on
yline(0, 'r--', 'LineWidth', 2, 'DisplayName', 'Limite de stabilité')
hold off
legend('Location', 'best')

% Compare the two methods
figure(3)
subplot(2,1,1)
plot(c1xx_values, real_parts_polyeig(:,1:4), 'LineWidth', 2)
grid on
title('Top 4 Modes - polyeig (QEP)')
xlabel('c_{1xx} (N·s/m)')
ylabel('Re(λ) (1/s)')
yline(0, 'r--', 'LineWidth', 1.5)
legend('Mode 1','Mode 2','Mode 3','Mode 4','Stabilité')
subplot(2,1,2)
plot(c1xx_values, real_parts_statespace(:,1:4), 'LineWidth', 2)
grid on
title('Top 4 Modes - State-Space')
xlabel('c_{1xx} (N·s/m)')
ylabel('Re(λ) (1/s)')
yline(0, 'r--', 'LineWidth', 1.5)
legend('Mode 1','Mode 2','Mode 3','Mode 4','Stabilité')

% Stability analysis
fprintf('=== STABILITY COMPARISON ===\n\n');

=== STABILITY COMPARISON ===

methods = {'polyeig (QEP)', 'State-Space'};
real_parts_all = {real_parts_polyeig, real_parts_statespace};
for method_idx = 1:2
fprintf('--- %s Method ---\n', methods{method_idx});
real_parts = real_parts_all{method_idx};
% Test at low, medium, high damping
test_indices = [1, round(length(c1xx_values)/2), length(c1xx_values)];
test_labels = {'Low', 'Medium', 'High'};
for idx = 1:length(test_indices)
i = test_indices(idx);
num_unstable = sum(real_parts(i,:) > 1e-6);
fprintf('%s damping (c1xx = %.2e):\n', test_labels{idx}, c1xx_values(i));
fprintf(' Unstable modes: %d / %d\n', num_unstable, num_eigenvalues);
if num_unstable > 0
fprintf(' Max real part: %.4e\n', max(real_parts(i,:)));
end
end
fprintf('\n');
end

--- polyeig (QEP) Method ---

Low damping (c1xx = 1.00e+05):

 Unstable modes: 0 / 12

Medium damping (c1xx = 4.95e+06):

 Unstable modes: 2 / 12

 Max real part: 1.2896e+00

High damping (c1xx = 1.00e+07):

 Unstable modes: 2 / 12

 Max real part: 2.1810e+00

--- State-Space Method ---

Low damping (c1xx = 1.00e+05):

 Unstable modes: 0 / 12

Medium damping (c1xx = 4.95e+06):

 Unstable modes: 2 / 12

 Max real part: 1.2896e+00

High damping (c1xx = 1.00e+07):

 Unstable modes: 2 / 12

 Max real part: 2.1810e+00

% Additional diagnostics
fprintf('=== SYSTEM PROPERTIES ===\n');

=== SYSTEM PROPERTIES ===

fprintf('Rotor speed parameter (q): %.2e\n', q);

Rotor speed parameter (q): 3.50e+06

fprintf('Natural frequency (w): %.2f rad/s\n', w);

Natural frequency (w): 588.46 rad/s

fprintf('Mass matrix condition number: %.2e\n', cond(M));

Mass matrix condition number: 4.37e+01

fprintf('Stiffness matrix condition number: %.2e\n', cond(K));

Stiffness matrix condition number: 8.11e+01

% Check for cross-coupling
K_symmetric = (K + K')/2;
K_skew = (K - K')/2;
fprintf('Stiffness symmetry index: %.2e\n', norm(K_skew,'fro')/norm(K,'fro'));

Stiffness symmetry index: 3.52e-02

fprintf('\nNote: Non-zero symmetry index indicates gyroscopic coupling\n');

Note: Non-zero symmetry index indicates gyroscopic coupling

fprintf(' This is NORMAL and EXPECTED in rotor dynamics!\n\n');

 This is NORMAL and EXPECTED in rotor dynamics!

fprintf('=== RECOMMENDATION ===\n');

=== RECOMMENDATION ===

fprintf('Use the polyeig (QEP) method results for your analysis.\n');

Use the polyeig (QEP) method results for your analysis.

fprintf('This method is specifically designed for quadratic eigenvalue\n');

This method is specifically designed for quadratic eigenvalue

fprintf('problems and handles gyroscopic effects correctly.\n')

problems and handles gyroscopic effects correctly.