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

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.