How to set the theorem of floquet for Linear variational equation with periodic equation?

I'm trying to analyze the stability regions of the equation of motion of a system with periodic coefficients. I have only written its polynomial equation and examined its roots, but my plot did not match the reference implementation correctly.

If I wanna write the Hill function and use the ode45 code, is it possible for someone to correct my code and my approach to finding the monodromy matrix and calculating the Floquet multipliers?

clc;
clear all;

% Define parameters
N = 2; % Number of blades

I_thetaoverI_b = 2; % Moment of inertia pitch axis over I_b 
I_psioverI_b = 2; % Moment of inertia yaw axis over I_b 

C_thetaoverI_b = 0.00; % Damping coefficient over I_b 
C_psioverI_b = 0.00; % Damping coefficient over I_b 

h = 0.3; % rotor mast height, wing tip spar to rotor hub [m]
hoverR = 0.34;
R = h / hoverR; % radius [m]

gamma = 4; % lock number

V_knots = 325; % the rotor forward velocity [knots]

% Convert velocity from [knots] to [meters per second]
% 1 knot = 0.51444 meters per second
V = V_knots * 0.51444;

% Calculate angular velocity in radians per second
omega_rad_s = V / R;

% Convert angular velocity from radians per second to RPM
% 1 radian per second = (60 / (2 * pi)) RPM
Omega = omega_rad_s * (60 / (2 * pi));

freq_1_over_Omega = 1 / Omega;

%the flap moment aerodynamic coefficients for large V
M_b = -(1/10)*V;
M_u = 1/6;

%the propeller aerodynamic coefficients
H_u = V/2;

% Frequency ranges
f_pitch= 5:3:140;
f_yaw= 5:3:140;


divergence_map = [];
Rdivergence_map = [];
unstable = [];

% Modify the loop to iterate over time points
for i = 1:length(f_pitch)
    for j = 1:length(f_yaw)

        phi_steps = linspace(0, pi, 100); % Time steps within one period
        for phi = phi_steps

            % Angular frequencies for the current time point
            w_omega_pitch = 2 * pi* f_pitch(i);
            w_omega_yaw = 2 * pi * f_yaw(j);

            K_psi = (w_omega_pitch^2) * I_psioverI_b;
            K_theta = (w_omega_yaw^2) * I_thetaoverI_b;


% Define inertia matrix [M]
M_matrix = [I_thetaoverI_b + 1 + cos(2*phi), -sin(2*phi);
             -sin(2*phi), I_psioverI_b + 1 - cos(2*phi)];

% Define damping matrix [D]
D11 =  C_thetaoverI_b + h^2*gamma*H_u*(1 - cos(2*phi)) - gamma*M_b*(1 + cos(2*phi)) - (2 + 2*h*gamma*M_u)*sin(2*phi);
D12 = h^2*gamma*H_u*sin(2*phi) + gamma*M_b*sin(2*phi) - 2*(1 + cos(2*phi)) - 2*h*gamma*M_u*cos(2*phi);
D21 = h^2*gamma*H_u*sin(2*phi) + gamma*M_b*sin(2*phi) + 2*(1 - cos(2*phi)) - 2*h*gamma*M_u*cos(2*phi);
D22 = C_psioverI_b + h^2*gamma*H_u*(1 + cos(2*phi)) - gamma*M_b*(1 - cos(2*phi)) + (2 + 2*h*gamma*M_u)*sin(2*phi);

D_matrix = [D11, D12;
            D21, D22];

% Define stiffness matrix [K]
K11 = K_theta - h*gamma*V*H_u*(1 - cos(2*phi)) + gamma*V*M_u*sin(2*phi);
K12 = -h*V*gamma*H_u*sin(2*phi) + gamma*V*M_u*(1 + cos(2*phi));
K21 = -h*gamma*V*H_u*sin(2*phi) - gamma*V*M_u*(1 - cos(2*phi));
K22 = K_psi - h*gamma*V*H_u*(1 + cos(2*phi)) - gamma*V*M_u*sin(2*phi);

K_matrix = [K11, K12;
            K21, K22];

% Compute the system matrices
           
M11 = M_matrix(1, 1); M12 = M_matrix(1, 2); M21 = M_matrix(2, 1); M22 = M_matrix(2, 2);
D11 = D_matrix(1, 1); D12 = D_matrix(1, 2); D21 = D_matrix(2, 1); D22 = D_matrix(2, 2);
K11 = K_matrix(1, 1); K12 = K_matrix(1, 2); K21 = K_matrix(2, 1); K22 = K_matrix(2, 2);

% Find the roots of the polynomial equation

    P0 = M11*M22-M12*M21;
    P1 = (- D11*M22*1j - D22*M11*1j + M12*D21*j + D12*M21*j);
    P2 = (D11*D22*(1j)^2 - K22*M11 - K11*M22 - D12*D21*(1j)^2 + M12*K21 + M21*K12);
    P3 = (D11*K22*1j - D12*K21*1j - D21*K12*1j + D22*K11*1j);
    P4 = K11*K22 - K12*K21;
    
    P = roots([P0, P1, P2, P3, P4]);

            r = 1 * P;

%Flutter

    for k = 1:length(r)

        if (real(r(k)) > 0)
           if (imag(r(k)) <= 0)
               unstable = [unstable; phi, K_psi, K_theta];

               % Proximity check for 1/Ω divergence
               if abs(real(r(k)) - freq_1_over_Omega) < 0.1
                  Rdivergence_map = [Rdivergence_map; phi, K_psi, K_theta];

               end

           end
        end              
                       
    end

 %Divergence
 
        if (real(det(K_matrix)) < 0)
               divergence_map = [divergence_map; phi, K_psi, K_theta];
        end
 
      end 
   end
end


% Plotting section

figure;
hold on;

scatter(unstable(:,2), unstable(:,3), 'filled');
scatter(divergence_map(:,2), divergence_map(:,3), 'filled', 'r');
scatter(Rdivergence_map(:,2), Rdivergence_map(:,3), 'filled', 'g');


xlabel('K\_psi');
ylabel('K\_theta');

title('Whirl Flutter Diagram');
legend('Flutter area', 'Divergence area', '1/Ω Divergence area');

hold off;

[Umar]

Hi @Nikoo,

I did execute your code and I end up getting this plot, so this attached plot below does not match with the plot as shown in your posted comments. Am I right?

[Umar]

Hi @Nikko,

I did read your short article attachment shared in your comments and implemented the framework based on it.

% Define the system's equation of motion with periodic coefficients
A = [1+cos(20), -sin(24); -sin(20), 1-cos(24)];
ch2xH = 1; % Placeholder value for ch2xH
yMg = 1; % Placeholder value for yMg
hy = 1; % Placeholder value for hy
h2yH = 1; % Placeholder value for h2yH
YM = 1; % Placeholder value for YM
yg = 1; % Placeholder value for yg
YVM = 1; % Placeholder value for YVM
hvVH = 1; % Placeholder value for hvVH
yW1 = 1; % Placeholder value for yW1
B = [ch2xH*(1-cos(24))-yMg*(1+cos(2*pi))-(2+2*hy)*sin(24), hy*sin(2*pi)
    +yg*sin(2*pi)-2*(1+cos(24))-2*hy*cos(2); 
    h2yH*sin(20)+YM*sin(20)+2*(1-cos(24))-2*hy, cos(24)];
C = [(1+cos(20))-yMg*(1-cos(2*pi))+(2+2*hy)*sin(24); -hy*sin(20)-YVM*(1-
cos(2*pi))];
D = -hvVH*(1-cos(2*pi))+yW1*sin(2*pi)-hy*sin(20)-YVM*(1-cos(2*pi));
% Implement Floquet theory techniques to analyze stability
% Perform stability analysis using Floquet theory
eigenvalues = eig(A); % Compute eigenvalues of matrix A
T = 2*pi; % Period of the coefficients
F = expm(A*T); % Compute the monodromy matrix using Floquet theory
floquet_eigenvalues = eig(F); % Compute Floquet multipliers
% Generate plots for stability regions on the K and Ky plane
% Plot stability regions with scatter
K_range = linspace(0, 10, 100); % Define range for K
Ky_range = linspace(0, 10, 100); % Define range for Ky
[X, Y] = meshgrid(K_range, Ky_range); % Create a grid of K and Ky values
Z = zeros(size(X)); % Initialize stability matrix
figure; % Create a new figure for scatter plot
hold on; % Enable hold to overlay scatter plot on the same figure
for i = 1:length(K_range)
    for j = 1:length(Ky_range)
        K = K_range(i);
        Ky = Ky_range(j);
        A_new = A - [0, K; Ky, 0]; % Update A matrix with K and Ky values
        F_new = expm(A_new*T); % Compute new monodromy matrix
        floquet_eigenvalues_new = eig(F_new); % Compute new Floquet multipliers
        Z(j, i) = max(abs(floquet_eigenvalues_new)); % Store maximum Floquet multiplier
        
        % Use scatter plot to visualize stability regions
        scatter(K, Ky, [], Z(j, i), 'filled');
    end
end
colorbar; % Add colorbar to show the scale of stability
xlabel('K'); % Label x-axis
ylabel('Ky'); % Label y-axis

The above provided code implements Floquet theory techniques.

% Define the system's equation of motion with periodic coefficients
A = [1+cos(20), -sin(24); -sin(20), 1-cos(24)];
ch2xH = 1; % Placeholder value for ch2xH
yMg = 1; % Placeholder value for yMg
hy = 1; % Placeholder value for hy
h2yH = 1; % Placeholder value for h2yH
YM = 1; % Placeholder value for YM
yg = 1; % Placeholder value for yg
YVM = 1; % Placeholder value for YVM
hvVH = 1; % Placeholder value for hvVH
yW1 = 1; % Placeholder value for yW1
B = [ch2xH*(1-cos(24))-yMg*(1+cos(2*pi))-(2+2*hy)*sin(24), hy*sin(2*pi)+yg*sin(2*pi)-2*(1+cos(24))-2*hy*cos(2); 
    h2yH*sin(20)+YM*sin(20)+2*(1-cos(24))-2*hy, cos(24)];
C = [(1+cos(20))-yMg*(1-cos(2*pi))+(2+2*hy)*sin(24); -hy*sin(20)-YVM*(1-cos(2*pi))];
D = -hvVH*(1-cos(2*pi))+yW1*sin(2*pi)-hy*sin(20)-YVM*(1-cos(2*pi));

In this section, the code defines the system's equation of motion with periodic coefficients. The matrix A represents the coefficients of the system's equations, which are periodic with a period of π. The other variables (ch2xH, yMg, hy, h2yH, YM, yg, YVM, hvVH, yW1) are placeholder values for the coefficients in the equations.

% Implement Floquet theory techniques to analyze stability
% Perform stability analysis using Floquet theory
eigenvalues = eig(A); % Compute eigenvalues of matrix A
T = 2*pi; % Period of the coefficients
F = expm(A*T); % Compute the monodromy matrix using Floquet theory
floquet_eigenvalues = eig(F); % Compute Floquet multipliers

In this section, the code implements Floquet theory techniques to analyze the stability of the system. It computes the eigenvalues of matrix A using the eig function. The period of the coefficients is defined as T = 2*pi. The monodromy matrix F is computed using the expm function, which calculates the matrix exponential. Finally, the code computes the Floquet multipliers by finding the eigenvalues of matrix F.

% Generate plots for stability regions on the K and Ky plane
% Plot stability regions with scatter
K_range = linspace(0, 10, 100); % Define range for K
Ky_range = linspace(0, 10, 100); % Define range for Ky
[X, Y] = meshgrid(K_range, Ky_range); % Create a grid of K and Ky values
Z = zeros(size(X)); % Initialize stability matrix
figure; % Create a new figure for scatter plot
hold on; % Enable hold to overlay scatter plot on the same figure
for i = 1:length(K_range)
    for j = 1:length(Ky_range)
        K = K_range(i);
        Ky = Ky_range(j);
        A_new = A - [0, K; Ky, 0]; % Update A matrix with K and Ky values
        F_new = expm(A_new*T); % Compute new monodromy matrix
        floquet_eigenvalues_new = eig(F_new); % Compute new Floquet multipliers
        Z(j, i) = max(abs(floquet_eigenvalues_new)); % Store maximum Floquet multiplier
        
        % Use scatter plot to visualize stability regions
        scatter(K, Ky, [], Z(j, i), 'filled');
    end
end
colorbar; % Add colorbar to show the scale of stability
xlabel('K'); % Label x-axis
ylabel('Ky'); % Label y-axis
title('Stability Regions on K-Ky Plane'); % Set title for the plot
hold off; % Disable hold to prevent further plots on the same figure

In this section, the code generates plots to visualize the stability regions on the K and Ky plane. It defines the range for K and Ky using the linspace function. Then, it creates a grid of K and Ky values using the meshgrid function. The stability matrix Z is initialized with zeros.

The code then enters a nested loop to iterate over each K and Ky value. For each combination, it updates the A matrix with the current K and Ky values and computes a new monodromy matrix F_new using the expm function. The Floquet multipliers for the new matrix are computed using the eig function.

The maximum Floquet multiplier for each combination of K and Ky is stored in the stability matrix Z. The scatter plot is used to visualize the stability regions, with the color of each point representing the maximum Floquet multiplier.

Finally, the code adds a colorbar to show the scale of stability, labels the x and y axes, and sets a title for the plot.

In nutshell, the provided code meets the requirements for analyzing the stability of a two-bladed rotor system using Floquet theory techniques. It computes the eigenvalues and Floquet multipliers, and generates plots to visualize the stability regions on the K and Ky plane.

title('Stability Regions on K-Ky Plane'); % Set title for the plot
hold off; % Disable hold to prevent further plots on the same figure

Please see attached plot.

<</matlabcentral/answers/uploaded_files/1746501/IMG_7318.jpeg>>

<Stability Plot>

[dpb]

@Umar I formatted the code and tried to convert the link but it 404'ed as is and tried to redirect -- that went to a shared community site but no content with my browser. Can you not just attach the image file with the paperclip to your post instead?

[Umar]

@dpb,

Thanks for your support. Next time, I will use that option. However, I would like your help on “Answer this question” button whenever I click on it, it defaults to form. May be you have to follow sequence in chronological order, like write something in form and then click “Answer this question” and then click “Submit” button.

[dpb]

I've never even tried to use a mobile device so I have no idea what the interface looks like for you...on the desktop with browser, there is an already open form into which one can begin to type; all the "Answer" button does is change focus into it for you; then you "Submit" after typing.

As noted before, with regular keyboard, I can Ctrl-E to toggle the code formatting; I gather that is not feasible on your touchscreen (although I don't mess with them, seems as though there would have to be some way to do those things, but I certainly don't know what it would be). Or, I can click on the leftmost icon in the "Code" section and the code line appears; it is not a toggle, however, you're in it until you refocus outside the area.

If you do get the form and begin typing, when you're ready, just "Submit"; don't click "Answer" at that point; that will probably bring you to a new form? Let's see what happens here if do that...

Well, in the browser if I have typed into the form and click "Answer" instead of "Submit", it brings up a warning that "Changes might be lost" and "Cancel" or "Continue" options -- either of which seems to just return to the existing form with present content still there. You can still edit in the form and then either "Submit" or "Cancel" the text entered.

Since the form is already there in the browser, I don't know that I'd ever actually tried the "Answer" button route first, but at least on the browser interface, all it seems to do is to put focus into the form for you instead of clicking into the form directly.

[Umar]

Hi @dpb, thanks for the tips, really appreciated.

[Umar]

Hi @Nikoo,

Addressing your query regarding, “Could you explain the logic behind A_new = A - [0, K; Ky, 0]; % Update A matrix with K and Ky values?”

The expression A_new = A - [0, K; Ky, 0] subtracts the values of K and Ky from the corresponding positions in the matrix A. This operation modifies the dynamics of the system represented by A, impacting stability analysis.

Now, coming back to, “I've updated my system to some extent, but I'm facing a major issue that I don't understand the mathematical logic behind it.”

Could you please elaborate further on “I'm facing a major issue that I don't understand the mathematical logic behind it”

[John D'Errico]

@Nikoo - you got multiple answers. I'm sorry, but we can't now remove your question.

Answer 1

[Umar]

Hi @Nikko,

Please see attached plot regarding “ in which the code computes the Floquet multipliers by finding the eigenvalues of matrix F”

Answer 2

[Umar]

Brief Analysis of troubleshooting your code

I had to do a lot of debugging in order to get your code executed properly. First, I had to deal with "Arrays have incompatible sizes for this operation" error, which was to make to sure that the arrays used in the calculations have compatible sizes. In the provided code, the issue arised from the calculation of the determinant of the K_matrix array. To resolve this issue, I had to modify the code to calculate the determinant of the K_matrix array using the det function. However, before calculating the determinant, I had to ensure that the K_matrix array is square, i.e., it has an equal number of rows and columns. So bear in mind , it is essential to ensure that the arrays used in calculations have compatible sizes. Before performing operations such as addition, subtraction, or multiplication, check that the arrays have the same dimensions or are compatible for the specific operation.Additionally, it is good practice to include error handling mechanisms, such as checking the size of arrays before performing operations, as shown in the solution above. This helps to catch any potential errors and handle them gracefully, providing informative messages to aid in debugging. Then, I faced with another error “index out of bounds error”, I added a check to verify that the arrays unstable, divergence_map, and Rdivergence_map have at least 3 columns before attempting to access elements at index 2 and 3. This precautionary measure helped me to prevent index out-of-bounds errors and ensures the arrays contain the necessary data for plotting. Afterwards, I faced the “Index in position 2 exceeds array bounds" error occurs when the code tries to access an element in an array that is beyond its bounds. In the given code, the error was caused by accessing elements in the unstable, divergence_map, and Rdivergence_map arrays without initializing them or updating their values. Lastly, removed the line scatter(Rdivergence_map(:, 2), Rdivergence_map(:, 3), 'filled', 'g'); since Rdivergence_map which was not used in the code. Here is your modified code,

N = 2; % Number of blades

I_thetaoverI_b = 2; % Moment of inertia pitch axis over I_b

I_psioverI_b = 2; % Moment of inertia yaw axis over I_b

C_thetaoverI_b = 0.00; % Damping coefficient over I_b

C_psioverI_b = 0.00; % Damping coefficient over I_b

h = 0.3; % rotor mast height, wing tip spar to rotor hub

hoverR = 0.34;

R = h / hoverR;

gamma = 4; % lock number

V_knots = 325; % the rotor forward velocity [knots]

% Convert velocity from [knots] to [meters per second]

% 1 knot = 0.51444 meters per second

V = V_knots * 0.51444;

% Calculate angular velocity in radians per second

omega_rad_s = V / R;

% Convert angular velocity from radians per second to RPM

% 1 radian per second = (60 / (2 * pi)) RPM

Omega = omega_rad_s * (60 / (2 * pi));

freq_1_over_Omega = 1 / Omega;

% the flap moment aerodynamic coefficients for large V

M_b = -(1 / 10) * V;

M_u = 1 / 6;

% the propeller aerodynamic coefficients

H_u = V / 2;

%%%%%%%%%%%the flap moment aerodynamic coefficients for small V

%M_b = -1*(1 + V^2)/8 ;

%M_u = V/4;

%the propeller aerodynamic coefficients

%H_u = (V^2/2)*log(2/V);

% Frequency ranges:

f_pitch = 5:5:140; % Frequency [Cycle/s or Hz]

f_yaw = 5:5:140; % Frequency [Cycle/s or Hz]

divergence_map = [];

unstable = [];

% Modify the loop to iterate over frequency points

for i = 1:length(f_pitch)

    for j = 1:length(f_yaw)

        % Angular frequencies for the current frequency points

        w_omega_pitch = 2 * pi * f_pitch(i);

        w_omega_yaw = 2 * pi * f_yaw(j);

        K_psi = (w_omega_pitch^2) * I_psioverI_b;

        K_theta = (w_omega_yaw^2) * I_thetaoverI_b;

        phi = 0;

        % Define inertia matrix [M]

        M_matrix = [I_thetaoverI_b + 1 + cos(2 * phi), -sin(2 * phi);

            -sin(2 * phi), I_psioverI_b + 1 - cos(2 * phi)];

        % Define damping matrix [D]

D11 = h^2 * gamma * H_u * (1 - cos(2 * phi)) - gamma * M_b * (1 + cos(2 * phi)) - (2 + 2 * h * gamma * M_u) * sin(2 * phi);

        D12 = h^2 * gamma * H_u * sin(2 * phi) + gamma * M_b * sin(2 * phi) - 2 * (1 + cos(2 * phi)) - 2 * h * gamma * M_u * cos(2 * phi);

        D21 = h^2 * gamma * H_u * sin(2 * phi) + gamma * M_b * sin(2 * phi) + 2 * (1 - cos(2 * phi)) - 2 * h * gamma * M_u * cos(2 * phi);

        D22 = h^2 * gamma * H_u * (1 + cos(2 * phi)) - gamma * M_b * (1 - cos(2 * phi)) + (2 + 2 * h * gamma * M_u) * sin(2 * phi);

        D_matrix = [D11, D12;
            D21, D22];

        % Define stiffness matrix [K]

        K11 = K_theta - h * gamma * V * H_u * (1 - cos(2 * phi)) + gamma * V * M_u * sin(2 * phi);

        K12 = -h * V * gamma * H_u * sin(2 * phi) + gamma * V * M_u * (1 + cos(2 * phi));

        K21 = -h * gamma * V * H_u * sin(2 * phi) - gamma * V * M_u * (1 - cos(2 * phi));

        K22 = K_psi - h * gamma * V * H_u * (1 + cos(2 * phi)) - gamma * V * M_u * sin(2 * phi);

        K_matrix = [K11, K12;
            K21, K22];

        % Check if K_matrix is square

        if size(K_matrix, 1) == size(K_matrix, 2)

            % Calculate determinant of K_matrix

            det_K = det(K_matrix);

            % Divergence condition

            if det_K < 0

                divergence_map = [divergence_map; phi, K_psi, K_theta];

            else

                unstable = [unstable; phi, K_psi, K_theta];

            end

        else

            disp('K_matrix is not square');

        end

    end

end

% Plot the Flutter and divergence maps

figure;

hold on;

scatter(unstable(:, 2), unstable(:, 3), 'filled');

scatter(divergence_map(:, 2), divergence_map(:, 3), 'filled', 'r');

xlabel('K\_psi');

ylabel('K\_theta');

title('Whirl Flutter Diagram');

legend('Flutter area', 'Divergence area');

hold off;

Please execute the modified code to see your attached plot.

I spent enough time troubleshooting this code, so for the future reference, I will advise, try to implement small lines of code first, display results to make sure that you are satisfied with them, it’s like solving a puzzle by putting small pieces together to get the bigger picture. Also, I do admire your enthusiasm solving a complex problem. You are quite motivated and this is why I started helping you out in your previous posts to help achieve your goal and didn’t want you to give up on it that you worked so hard for.

[Umar]

Hi @Nikko,

In your code, you have already computed the eigenvalues of the matrix A using the eig(A) function. These eigenvalues represent the stability of the system at each frequency point. However, to update the A matrix, you need to compute the Floquet multipliers which can be obtained by computing the eigenvalues of the monodromy matrix, which is calculated using the expm(A*T) function. Here, T represents the period of the coefficients in your system. In your code, you have set T to 2*pi. So, to update the A matrix, you can use the following code:

% Perform stability analysis using Floquet theory

eigenvalues = eig(A); % Compute eigenvalues of matrix A

T = 2*pi; % Period of the coefficients

F = expm(A*T); % Compute the monodromy matrix using Floquet theory

floquet_eigenvalues = eig(F); % Compute Floquet multipliers

% Update A matrix

A_new = A + log(floquet_eigenvalues)/T;

Now, compute the Floquet multipliers by taking the logarithm of the eigenvalues of the monodromy matrix and dividing by the period T. Then, update the A matrix by adding the computed values. Please bear in mind that the monodromy matrix F obtained represents the linearized dynamics of the system over one period. It contains information about how the system evolves from one period to the next. The eigenvalues of F, known as Floquet multipliers, provide insights into the stability of the periodic solution. If all Floquet multipliers have a magnitude less than 1, the periodic solution is stable. I hope this helps! Let me know if you have any further questions.

[Umar]

Hi @Nikko,

What makes you think that your diagram isnt coming out correctly? Please share your analysis and opinion by elaborating on this part.

[Umar]

Hi @Nikko,

Let me digest about what you said. I have to take some time out if my schedule to review it and then get back to you. Hope, you are not in rush.

[Nikoo]

for sure. Thank you for taking the time for me.

[Umar]

Hi @Nikko,

Please see my response to your comments below.

Interpretation of NASA document and your code

Based on attached NASA document provided by you, your’s construction of the periodic matrix A aligns with the description in the NASA document regarding how to encapsulate dynamics influenced by periodic coefficients.The inclusion of `A_top` and `A_bottom` effectively captures the dynamics pertinent to pitch and yaw motions. The use of Floquet theory to compute eigenvalues is consistent with established methods for analyzing systems with periodic coefficients by calculating both `eigenvalues` and `floquet_eigenvalues`, and you are correctly assessing stability boundaries as indicated in the NASA document. Your checks for divergence conditions using `det(K_matrix) < 0` directly reflect the instability criteria mentioned in the NASA document regarding low stiffness scenarios and this matches with the proximity check for eigenvalues against a frequency is aligned with recognizing divergence-like instabilities at multiples of fundamental frequencies. I would like to share my additional insights about interaction between pitch and yaw dynamics which can lead to complex behaviors not fully captured by linear models alone; consider exploring nonlinear effects if applicable. Also, understanding real-world implications—such as how these findings could influence rotorcraft design—can add value to your analysis.

Now, coming back to, “why my chart doesn't correctly display the divergence range because it is not symmetric and starts from zero, whereas the reference does not start from zero” &, “in the flutter range, I need to obtain the boundary curve, which is not correct.”

These are the potential issues that I come up with after analysis of your recent code.

Symmetry in Divergence Range

Review your code again to make sure that the values for `K_psi` and `K_theta` used in plotting are appropriately scaled or transformed to achieve symmetry.

Starting Point of Ranges

 The issue with starting from zero may stem from how `divergence_map`, `unstable`, and `Rdivergence_map` are populated. Again, review you code to make sure that these arrays are filled based on certain conditions without normalization or scaling, their plotted values might not align with the reference chart.

Plotting Logic

Go back to your code and find out if all relevant data points are included and correctly represented. If any values fall outside of the expected range or are incorrectly indexed, it can lead to misleading visualizations.

To address these issues and improve the accuracy of your plot:

Consider normalizing your data before plotting to ensure that all values are within a consistent range. For example, you could transform `K_psi` and `K_theta` using min-max normalization or z-score normalization.

   K_psi_normalized = (K_psi - min(K_psi)) / (max(K_psi) - min(K_psi));

   K_theta_normalized = (K_theta - min(K_theta)) / (max(K_theta) - min(K_theta));

Modify the plotting section to ensure that both `K_psi` and `K_theta` are plotted relative to their respective ranges in a way that allows for better visual symmetry.

   % Ensure appropriate limits for axes

   xlim([min(K_psi) max(K_psi)]);

   ylim([min(K_theta) max(K_theta)]);

Before plotting, add debugging outputs to check the contents of your divergence maps:

   disp('Unstable Points:');

   disp(unstable);

   disp('Divergence Map:');

   disp(divergence_map);

   disp('R Divergence Map:');

   disp(Rdivergence_map);

Revisit the conditions under which you populate `divergence_map`, `unstable`, and `Rdivergence_map`. Make sure they accurately reflect the intended physical phenomena being modeled.

Addressing your query regarding, “in the flutter range, I need to obtain the boundary curve, which is not correct.”

Matrix Inversion

The lines where `M_matrix` is defined involve inverting it (`-inv(M_matrix) * K_matrix`). If `M_matrix` is singular or poorly conditioned, this could lead to numerical instability. Ensure that `M_matrix` is invertible before performing this operation.

Floquet Theory Implementation

The computation of the monodromy matrix (`F = expm(A*T)`) followed by its eigenvalues should ensure that `A` is correctly defined. Make sure that `A` contains valid entries that accurately represent the system's dynamics. Also, the update of the `A` matrix with `A_new = A + log(floquet_eigenvalues)/T;` might not be appropriate. The logarithm of eigenvalues should be taken carefully as they may not be directly applicable to modify matrix `A`. This could distort the dynamics being analyzed.

Stability Conditions

 In your flutter condition check:

     if (abs(floquet_eigenvalues_new) > 1)

Make sure that you are checking each eigenvalue individually, as `floquet_eigenvalues_new` might be a vector. It should be modified to iterate through each eigenvalue:

     for ev = floquet_eigenvalues_new'

         if abs(ev) > 1

             unstable = [unstable; K_psi, K_theta];

         end

     end

Divergence Check

 Similar to the flutter condition, ensure that you are iterating over each eigenvalue in your divergence check:

     for ev = eig(F)'

         if (abs(ev) - freq_1_over_Omega) < 0.1

             Rdivergence_map = [Rdivergence_map; K_psi, K_theta];

         end

     end

Determinant Check for Divergence

The condition for divergence based on the determinant of `K_matrix` seems appropriate; however, it might be beneficial to include checks on the physical meaning of this condition (e.g., ensuring that stiffness values are within expected ranges).

Parameter Initialization

Ensure all parameters such as `I_thetaoverI_b`, `I_psioverI_b`, etc., are initialized correctly based on the physical system being modeled.

Plotting

 After collecting data into `unstable`, `divergence_map`, and `Rdivergence_map`, ensure these arrays are not empty before plotting to avoid runtime errors.

Suggested Corrections

Here’s a refined version of critical sections of your code:

` % Ensure M_matrix is invertible before inversion

if det(M_matrix) == 0

    error('M_matrix is singular and cannot be inverted.');

end

% Calculate Floquet multipliers correctly

F = expm(A*T);

floquet_eigenvalues = eig(F);

% Stability condition check

for ev = floquet_eigenvalues'

    if abs(ev) > 1

        unstable = [unstable; K_psi, K_theta];

    end

end

% Divergence proximity check

for ev = eig(F)'

    if (abs(ev) - freq_1_over_Omega) < 0.1

        Rdivergence_map = [Rdivergence_map; K_psi, K_theta];

    end

end

% Ensure plotting only occurs if data exists

if ~isempty(unstable)

    figure;

    hold on;

    scatter(unstable(:,1), unstable(:,2), 'filled');

    scatter(divergence_map(:,1), divergence_map(:,2), 'filled', 'r');

    scatter(Rdivergence_map(:, 1), Rdivergence_map(:, 2), 'filled', 'g');

    xlabel('K\_psi');

    ylabel('K\_theta');

    title('Whirl Flutter Diagram');

    legend('Flutter area', 'Divergence area', '1/Ω Divergence area');

    hold off;

end

Validation of Results

Once you have corrected your code, validate your results against known benchmarks or theoretical predictions to ensure accuracy. Use more stable numerical methods or libraries for matrix operations, especially for large matrices or those prone to singularities.For better clarity in your plots, consider adding grid lines or different markers/colors for different conditions.

[dpb]

"... If `M_matrix` is singular or poorly conditioned, this could lead to numerical instability. ..."

MATLAB numerical inversion will provide a warning if the matrix is very badly condition (or the system is rank deficient or the like). Nevertheless, the point is valid and nothing lost by doing some specific condition testing.

"... the update of the `A` matrix with `A_new = A + log(floquet_eigenvalues)/T;` might not be appropriate. "

Unlike you, I've not even tried to read the paper, but I would think there's no chance this actually represents the transient dynamics...

if det(M_matrix) == 0

will only help if the M_matrix is identically singular; it will fail anything for which the determinant is even

>> realmin
ans =
  2.2251e-308
>> d=realmin;
>> d==0
ans =
  logical
   0
>>
end

The conditioning tests will be far more likely to be useful...the only way the above is likely to occur is if two columns of the matrix are identically multiples of one another. That can happen if it were to be constructed incorrectly, but that's about the only case.

"Use more stable numerical methods or libraries for matrix operations, ..."

Perhaps better solution techniques are available, but it is not likely you'll find a better set of methods for matrix calculations than those implemented in MATLAB. Real pros built these...

[Umar]

Hi @dpb,

Thanks for your contribution, highly appreciated.

[dpb]

"... not likely you'll find a better set of methods for matrix calculations than those implemented in MATLAB. Real pros built these..."

To get an idea of the sort of folks and on what MATLAB is based, you can look at a couple of Cleve's blog posts; these just scratch the surface, but give an idea of the sophistication of the underpinnings...

<Brief MATLAB History> and <The Gatlinburg and Householder Symposia>. Poking at <Cleve's Corner> is extremely educational although can be a deep diversion from what one's supposed to be doing! <vbg>

The ACM SIG on Programming Languages planned on a 2020 conference on the History of Programming Languages and Cleve wrote a paper on MATLAB for it which he posted piecewise as he completed sections...at some time in the distant past I collected the sequence of post links as follows...

<Part 1 Linear Algebra Handbook>

<Part 2 EISPACK>

<Part 3 LINPACK>

<Part 4 PC-MATLAB>

<Part 5 Modern MATLAB-1>

<Part 6 Modern MATLAB-2>

<Part 7 MATLAB Environment>

[Umar]

Hi @dpb,

Your insights are invaluable, and I am grateful for your guidance in navigating this complex subject. Should you have any further recommendations or discussions in mind, please feel free to share.It’s like The "rinse and repeat" instruction on shampoo bottle as you mentioned before. Still cracks me up.🤣