How can I evaluate the transition matrix using the fourth order Runge-Kutta ODEs

I have a two-degree-of-freedom system. for solving the system I am trying to use RK4. Can anyone help me with correctly writing the transition matrix in my matlab code according to the definition in this reference?

% Time settings
t0 = 0;
tf = pi / omega_rad_s^2;
N = 100;
ts = (tf - t0) / N;
t = linspace(t0, tf, N+1);

%Initial condition
y01 = 1;
y02 = 0;
y03 = 0;
y04 = 0;

y = [y01; y02; y03; y04]; % IC vector

% Solve the system using RK4
    for i = 1:N

        k1 = Function_system(t(i), y(:, i));
        k2 = Function_system(t(i) + ts/2, y(:, i) + ts/2 * k1);
        k3 = Function_system(t(i) + ts/2, y(:, i) + ((sqrt(2)-1)/2) * ts * k1 + (1 - (1/sqrt(2))) * ts * k2);
        k4 = Function_system(t(i) + ts, y(:, i) - (1/sqrt(2)) * ts * k2 + (1 + 1/sqrt(2)) * ts * k3);
        y(:, i+1) = y(:, i) + ts/6 * (k1 + 2*(1 - (1/sqrt(2)))*k2 + 2*(1 + (1/sqrt(2)))*k3 + k4);

    end

[Torsten]

Here your function "Function_system" has three inputs:

k1 = Function_system(t(i), y(:, i), K_theta);

Here it has two:

k2 = Function_system(t(i) + ts/2, y(:, i) + ts/2 * k1);

What is correct ?

[Nikoo]

k1 = Function_system(t(i), y(:, i));
k2 = Function_system(t(i) + ts/2, y(:, i) + ts/2 * k1);

All should have 2

Answer 1

[Torsten]

Seems to work. What's your problem ?

Don't use the transition matrix to integrate in one pass. If you have to, generate it with symbolic inputs to "Function_system".

% Time settings
t0 = 0;
tf = 3;
N = 100;
ts = (tf - t0) / N;
t = linspace(t0, tf, N+1);
%Initial condition
y01 = 1;
y02 = 1;
y03 = 1;
y04 = 1;
y = [y01; y02; y03; y04]; % IC vector
% Solve the system using RK4
    for i = 1:N
        k1 = Function_system(t(i), y(:, i));
        k2 = Function_system(t(i) + ts/2, y(:, i) + ts/2 * k1);
        k3 = Function_system(t(i) + ts/2, y(:, i) + ((sqrt(2)-1)/2) * ts * k1 + (1 - (1/sqrt(2))) * ts * k2);
        k4 = Function_system(t(i) + ts, y(:, i) - (1/sqrt(2)) * ts * k2 + (1 + 1/sqrt(2)) * ts * k3);
        y(:, i+1) = y(:, i) + ts/6 * (k1 + 2*(1 - (1/sqrt(2)))*k2 + 2*(1 + (1/sqrt(2)))*k3 + k4);
    end
    plot(t,y(3,:))

function dydt = Function_system(t,y)
  dydt = [-y(1),-2*y(2),y(3),2*y(4)].';
end

[Nikoo]

Yes, that works, but I need its transition matrix(4x4 for 2 DOF system) to furthur examine its eigenvalues. could you write it symbolically for me to see how it can be done?

[Torsten]

So your system is of the form

y' = A(t)*y

?

Here is an example for a system of two equations:

syms t dt
A(t) = [1 t; cos(t) 2];
n = size(A,1);
E = A(t+dt/2)*(eye(n)+dt/2*A(t));
F = A(t+dt/2)*(eye(n)+(-1/2+1/sym(sqrt(2)))*dt*A(t)+(1-1/sym(sqrt(2)))*dt*E);
G = A(t+dt)*(eye(n)-dt/sym(sqrt(2))*E+(1+1/sym(sqrt(2)))*dt*F);
K = eye(n)+dt/6*(A(t)+2*(1-1/sym(sqrt(2)))*E+2*(1+1/sym(sqrt(2)))*F+G);
K = simplify(K);
eig(K)

[Nikoo]

% Time settings
t0 = 0;
tf = pi / omega_rad_s^2;
N = 100;
ts = (tf - t0) / N;
t = linspace(t0, tf, N+1);
%Initial condition
y01 = 1;
y02 = 0;
y03 = 0;
y04 = 0;
y = [y01; y02; y03; y04]; % IC vector
% Solve the system using RK4
    for i = 1:N
        k1 = Function_system(t(i), y(:, i));
        k2 = Function_system(t(i) + ts/2, y(:, i) + ts/2 * k1);
        k3 = Function_system(t(i) + ts/2, y(:, i) + ((sqrt(2)-1)/2) * ts * k1 + (1 - (1/sqrt(2))) * ts * k2);
        k4 = Function_system(t(i) + ts, y(:, i) - (1/sqrt(2)) * ts * k2 + (1 + 1/sqrt(2)) * ts * k3);
        y(:, i+1) = y(:, i) + ts/6 * (k1 + 2*(1 - (1/sqrt(2)))*k2 + 2*(1 + (1/sqrt(2)))*k3 + k4);
    end

function [dy_dt] = function_system(t, y)

    % Define inertia matrix [M]
    M = [M11 M12; M21 M22];

    % Define damping matrix [D]
    D = [D11 D12; D21 D22];

    % Define stiffness matrix [K]
    K = [K11 K12; K21 K22];

    % Calculate the inverse matrix-vector product and assign correctly
    acceleration = -inv(M) * (D * [y(2); y(4)] + K * [y(1); y(3)]);

    % System equations
    dy_dt = [y(2) ; acceleration(1); y(4); acceleration(2)] ;

end

I don't have A separately. I wrote the function for my system symbolicaly. how can i define A? And should the transition matrix be written inside the loop?

[Walter Roberson]

Note that

acceleration = -inv(M) * (D * [y(2); y(4)] + K * [y(1); y(3)]);

is better written

acceleration = -M \ (D * [y(2); y(4)] + K * [y(1); y(3)]);

[Torsten]

syms M11 M12 M21 M22 D11 D12 D21 D22 K11 K12 K21 K22 
syms y [1 4]
    % Define inertia matrix [M]
    M = [M11 M12; M21 M22];
    % Inverse of M
    Minv = inv(M);
    % Define damping matrix [D]
    D = [D11 D12; D21 D22];
    % Define stiffness matrix [K]
    K = [K11 K12; K21 K22];
    % Calculate the inverse matrix-vector product and assign correctly
    acceleration = -Minv * (D * [y(2); y(4)] + K * [y(1); y(3)]);
    % System equations
    dy_dt = [y(2) ; acceleration(1); y(4); acceleration(2)] ;
b = zeros(4,1);
vars = y;
[A,~] = equationsToMatrix(dy_dt==b,vars)

After you get A here, you can use "subs" to give numerical values to all the symbolic variables M, D and K. The following code can also be performed numerically by giving a value to "dt" and removing the "sym" commands.

syms dt
n = size(A,1);
E = A*(eye(n)+dt/2*A);
F = A*(eye(n)+(-1/2+1/sym(sqrt(2)))*dt*A+(1-1/sym(sqrt(2)))*dt*E);
G = A*(eye(n)-dt/sym(sqrt(2))*E+(1+1/sym(sqrt(2)))*dt*F);
K = eye(n)+dt/6*(A+2*(1-1/sym(sqrt(2)))*E+2*(1+1/sym(sqrt(2)))*F+G);
K = simplify(K)

[Torsten]

You can directly use my code - so why do you still use your version ?

Your A is constant and does not depend on t - thus A(t) = A(t+dt/2) = A(t+dt) == A for all values of t.

[Nikoo]

I tried to write it without using functions, as you suggested but I belive A is not constant : Am I making a mistake?

% Time settings

t0 = 0;
tf = pi ;
N = 1000;
ts = (tf - t0) / N;
t = linspace(t0, tf, N+1);

%Initial condition
I = eye(4);

%%%%%% Compute the monodromy matrix %%%%%%%
% Solve the system using RK4

 for i = 1:N
    MO = I;

    ct = cos(2*t);
    st = sin(2*t);

    % Define inertia matrix [M]
    M11 = 3 + 1 + ct;
    M12 = -st;
    M21 = -st;
    M22 = 3 + 1 - ct;

    M  = @(t)[M11(t) M12(t); M21(t) M22(t)];

    % Define damping matrix [D]
    D11 =  (1 - ct) - (1 + ct) - 2 * st;
    D12 =  st + 0.4*st - 2*(1 + ct) - 2*ct;
    D21 =  st + 0.4st + 2*(1 - ct) - 2*ct;
    D22 =  (1 + ct) - (1 - ct) + 2 * st;

    D = @(t)[D11(t) D12(t); D21(t) D22(t)];

    % Define stiffness matrix [K]
    K11 = - (1 - ct) + st;
    K12 = -st + 0.5*(1 + ct);
    K21 = -st - 0.5*(1 - ct);
    K22 = - (1 + ct) - st;

    K = @(t)[K11(t) K12(t); K21(t) K22(t)];

    % Compute the terms -M \ K and -M \ C 
    MK =@(t) -M(t) \ K(t); % Solves M * y = K for y
    MC =@(t) -M(t) \ C(t); % Solves M * y = C for y

    % Compute A(phi)
    A = @(t)[ 0, 1, 0, 0;
    MK(1,1), MC(1,1), MK(1,2), MC(1,2);
    0, 0, 0, 1;
    MK(2,1), MC(2,1), MK(2,2), MC(2,2)];

    % Compute A at specific points
    A_mid = A(t + ts/2);    % A(phi + h/2)
    A_next = A(t + ts);     % A(phi + h)

    % Compute intermediate matrices E, F, G
    E = A_mid * (I + (1/2 * ts * A));
    F = A_mid * (I + (-1/2 + 1/sqrt(2)) * ts * A + (1 - 1/sqrt(2)) * ts * E);
    G = A_next * (I - ts/sqrt(2) * E + (1 + 1/sqrt(2)) * ts * F); 

    % Update MO based on the approximations above
    MO = MO * (I + ts/6 * (A + 2 * (1 - 1/sqrt(2)) * E + 2 * (1 + 1/sqrt(2)) * F ) + G);

end

% computation of floquet multipliers
eigenvalues= eig(MO);

[Torsten]

I gave you code that handles the case that you have a constant matrix A and the case that A depends on t. Simply use it.

[Torsten]

syms y [1 4]
syms t dt
ct = cos(2*t);
st = sin(2*t);
% Define inertia matrix [M]
M11 = 3 + 1 + ct;
M12 = -st;
M21 = -st;
M22 = 3 + 1 - ct;
M  = [M11 M12;M21 M22];
Minv = inv(M);
% Define damping matrix [D]
D11 =  (1 - ct) - (1 + ct) - 2 * st;
D12 =  st + 0.4*st - 2*(1 + ct) - 2*ct;
D21 =  st + 0.4*st + 2*(1 - ct) - 2*ct;
D22 =  (1 + ct) - (1 - ct) + 2 * st;
D = [D11 D12; D21 D22];
% Define stiffness matrix [K]
K11 = - (1 - ct) + st;
K12 = -st + 0.5*(1 + ct);
K21 = -st - 0.5*(1 - ct);
K22 = - (1 + ct) - st;
K = [K11 K12; K21 K22];
% Calculate the inverse matrix-vector product and assign correctly
acceleration = -Minv * (D * [y(2); y(4)] + K * [y(1); y(3)]);
% System equations
dy_dt = [y(2) ; acceleration(1); y(4); acceleration(2)] ;
b = zeros(4,1);
vars = y;
[A(t),~] = equationsToMatrix(dy_dt==b,vars);
n = size(A,1);
E = A(t+dt/2)*(eye(n)+dt/2*A(t));
F = A(t+dt/2)*(eye(n)+(-1/2+1/sym(sqrt(2)))*dt*A(t)+(1-1/sym(sqrt(2)))*dt*E);
G = A(t+dt)*(eye(n)-dt/sym(sqrt(2))*E+(1+1/sym(sqrt(2)))*dt*F);
K = eye(n)+dt/6*(A(t)+2*(1-1/sym(sqrt(2)))*E+2*(1+1/sym(sqrt(2)))*F+G);
% Time settings
t0 = 0;
tf = pi ;
N = 20;
ts = (tf - t0) / N;
T = linspace(t0, tf, N+1);
K = subs(K,dt,ts);
% I wouldn't trust in the result of this N-fold matrix multiplication
Knum = double(subs(K,t,0));
for i = 1:N-1
  Knum = double(subs(K,t,i*ts))*Knum;   
end
format long
eig(Knum)

ans = 4×1

1.0e+12 * 0.000000000000001 5.859349701937719 0.000000000000000 -0.000000000000000

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

[Nikoo]

Thank you for your guidance. I really appriciate it.

I have one question:

Can the RK4 method be replaced with ode45, or is RK4 better suited for periodic systems?

[Torsten]

What is your aim ? Solving the system from above or examining the RK4 scheme ?

If you want to solve the system, use "ode45". But I thought your task is to examine the eigenvalues of K derived from the RK4 scheme ...

[Nikoo]

Yes, my goal is to find the transition matrix and analyze its eigenvalues.

I thought that ODE45 might be a suitable function, but after doing some research, I think using RK4 might give me more accurate results.

[Torsten]

I thought that ODE45 might be a suitable function, but after doing some research, I think using RK4 might give me more accurate results.

No, but how to find the numerical transition matrix for "ODE45" ? Do you know how to compute K for the scheme implemented in "ODE45" ?

Answer 2

[Torsten]

If it's still of interest: this code seems to work correctly. I don't know why using

n = size(A,1);
E = A(t+dt/2)*(eye(n)+dt/2*A(t));
F = A(t+dt/2)*(eye(n)+(-1/2+1/sym(sqrt(2)))*dt*A(t)+(1-1/sym(sqrt(2)))*dt*E);
G = A(t+dt)*(eye(n)-dt/sym(sqrt(2))*E+(1+1/sym(sqrt(2)))*dt*F);
K = eye(n)+dt/6*(A(t)+2*(1-1/sym(sqrt(2)))*E+2*(1+1/sym(sqrt(2)))*F+G);

gives wrong results.

format long
syms y [4 1]
syms t dt
ct = cos(2*t);
st = sin(2*t);
% Define inertia matrix [M]
M11 = 3 + 1 + ct;
M12 = -st;
M21 = -st;
M22 = 3 + 1 - ct;
M  = [M11 M12;M21 M22];
Minv = inv(M);
% Define damping matrix [D]
D11 =  (1 - ct) - (1 + ct) - 2 * st;
D12 =  st + 0.4*st - 2*(1 + ct) - 2*ct;
D21 =  st + 0.4*st + 2*(1 - ct) - 2*ct;
D22 =  (1 + ct) - (1 - ct) + 2 * st;
D = [D11 D12; D21 D22];
% Define stiffness matrix [K]
K11 = - (1 - ct) + st;
K12 = -st + 0.5*(1 + ct);
K21 = -st - 0.5*(1 - ct);
K22 = - (1 + ct) - st;
K = [K11 K12; K21 K22];
% Calculate the inverse matrix-vector product and assign correctly
acceleration = -Minv * (D * [y(2); y(4)] + K * [y(1); y(3)]);
% System equations
dy_dt = [y(2) ; acceleration(1); y(4); acceleration(2)] ;
b = zeros(4,1);
vars = y;
[A(t),~] = equationsToMatrix(dy_dt==b,vars);
% Compute transition matrix K
Function_system = @(t,y)A(t)*y;
k1 = Function_system(t, y);
k2 = Function_system(t + dt/2, y + dt/2 * k1);
k3 = Function_system(t + dt/2, y + ((sqrt(2)-1)/2) * dt * k1 + (1 - (1/sqrt(2))) * dt * k2);
k4 = Function_system(t + dt, y - (1/sqrt(2)) * dt * k2 + (1 + 1/sqrt(2)) * dt * k3);
ynew = y + dt/6 * (k1 + 2*(1 - (1/sqrt(2)))*k2 + 2*(1 + (1/sqrt(2)))*k3 + k4);
b = zeros(4,1);
vars = y;
[K,~] = equationsToMatrix(ynew==b,vars);
% Compute numerical solution for tend = pi and y0 = [1;1;1;1] using the
% transition matrix K
t0 = 0;
tf = pi ;
N = 20;
ts = (tf - t0) / N;
T = linspace(t0, tf, N+1);
K = subs(K,dt,ts);
Knum = double(subs(K,t,0));
for i = 1:N-1
  Knum = double(subs(K,t,i*ts))*Knum;   
end
y_at_pi_with_transition_matrix = Knum*[1;1;1;1]

y_at_pi_with_transition_matrix = 4×1

9.754957068626140 2.363232167775414 -0.641295161941964 -2.287853242583158

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

% Compute eigenvalues
eig(Knum)

ans = 
  3.282252311329557 + 1.080784504109312i
  3.282252311329557 - 1.080784504109312i
 -1.030044815844306 + 0.000000000000000i
 -0.039890272044380 + 0.000000000000000i

% Compute numerical solution for tend = pi and y0 = [1;1;1;1] with ode45 to check the 
% result with transition matrix
fun = @(t,y)double(A(t))*y;
y0 = [1;1;1;1];
[T,Y] = ode45(fun,[0 pi],y0);
y_at_pi_with_ode45 = Y(end,:).'

y_at_pi_with_ode45 = 4×1

9.755085521136586 2.363155160915888 -0.641612321978496 -2.288142249726387

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

[Torsten]

Use

n = size(A(0),1);

instead of

n = size(A,1);

and the older code will work, too.

[Nikoo]

Hi, thanks a lot, Torsten, for your help.

I have a few questions: I can't run these two functions:

[A(t),~] = equationsToMatrix(dy_dt==b,vars);
[K,~] = equationsToMatrix(ynew==b,vars);

and they give me an error saying they're not (equationsToMatrix) in my toolbox. Is there another way to write them?

Also, I didn't understand why the double function is being used. Could you explain?

[Walter Roberson]

equationsToMatrix is part of the Symbolic Toolbox, which you must have or else

syms y [4 1]

would not work.

In particular, we could consider the possibility that you have access to syms because you (hypothetically) have Maple installed -- but the Maple version of syms does not support the [4 1] notation.

[Torsten]

@Walter Roberson

Code from @Nikoo never used symbolic computations ...

@Nikoo

format long
t0 = 0;
tf = pi ;
N = 20;
dt = (tf - t0) / N;
T = linspace(t0, tf, N+1);
Kmat = K(0,dt);
for i = 1:N-1
  Kmat = K(i*dt,dt)*Kmat;   
end
y_at_pi_with_transition_matrix = Kmat*[1;1;1;1]

y_at_pi_with_transition_matrix = 4×1

9.754957068626140 2.363232167775413 -0.641295161941960 -2.287853242583156

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

% Compute eigenvalues
eig(Kmat)

ans = 
  3.282252311329557 + 1.080784504109311i
  3.282252311329557 - 1.080784504109311i
 -1.030044815844306 + 0.000000000000000i
 -0.039890272044380 + 0.000000000000000i

function Mat_K = K(t,dt)
  n = 4;
  At = A(t);
  Atpdth = A(t+dt/2);
  Atpdt = A(t+dt);
  E = Atpdth*(eye(n)+dt/2*At);
  F = Atpdth*(eye(n)+(-1/2+1/sqrt(2))*dt*At+(1-1/sqrt(2))*dt*E);
  G = Atpdt*(eye(n)-dt/sqrt(2)*E+(1+1/sqrt(2))*dt*F);
  K = eye(n)+dt/6*(At+2*(1-1/sqrt(2))*E+2*(1+1/sqrt(2))*F+G);
  Mat_K = K;
end

function Mat_A = A(t)
   ct = cos(2*t);
   st = sin(2*t);
   % Define inertia matrix [M]
   M11 = 3 + 1 + ct;
   M12 = -st;
   M21 = -st;
   M22 = 3 + 1 - ct;
   M  = [M11 M12;M21 M22];
   % Define damping matrix [D]
   D11 =  (1 - ct) - (1 + ct) - 2 * st;
   D12 =  st + 0.4*st - 2*(1 + ct) - 2*ct;
   D21 =  st + 0.4*st + 2*(1 - ct) - 2*ct;
   D22 =  (1 + ct) - (1 - ct) + 2 * st;
   D = [D11 D12; D21 D22];
   % Define stiffness matrix [K]
   K11 = - (1 - ct) + st;
   K12 = -st + 0.5*(1 + ct);
   K21 = -st - 0.5*(1 - ct);
   K22 = - (1 + ct) - st;
   K = [K11 K12; K21 K22];
   % Compute the terms -M \ K and -M \ C 
   MK =-M \K; % Solves M * y = K for y
   MD =-M \D; % Solves M * y = D for y
   % Compute A(t)
   A = [ 0, 1, 0, 0;
         MK(1,1), MD(1,1), MK(1,2), MD(1,2);
         0, 0, 0, 1;
         MK(2,1), MD(2,1), MK(2,2), MD(2,2)];
   Mat_A = A;
end

[Nikoo]

Thank you very much, @Torsten, for all your guidance.

[Nikoo]

Hi @Torsten I intend to implement the following diagram, which corresponds to a 2-degree-of-freedom system with periodic coefficients, but the boundaries do not fully match the reference.

clc
clear all ;

% Define parameters 
global I_thetaoverI_b I_psioverI_b h gamma H_u M_b M_u V Omega omega_rad_s 

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 = 1 * 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)) ;

%%%%%%%%% Aerodynamic coefficients calculations
M_u = 0.25*sqrt(2) - 0.25*log(1 + sqrt(2));
M_b = 0.0625*sqrt(2) - 0.1875*log(1 + sqrt(2));
H_u = 0.5*log(1 + sqrt(2));

% Frequency ranges:
f_pitch= 0.06:0.01:0.5 ; % Frequency [Cycle/s or Hz]
f_yaw= 0.06:0.01:0.5; % Frequency [Cycle/s or Hz]

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

% Modify the loop to iterate over frequency points

for kk1 = 1:length(f_pitch)
    for kk2 = 1:length(f_yaw)

% Angular frequencies for the current frequency points
w_omega_pitch = 2 * pi * f_pitch(kk1);
w_omega_yaw = 2 * pi * f_yaw(kk2);

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

% Time settings

t0 = 0;
tf = pi / omega_rad_s^2;
N = 20;
dt = (tf - t0) / N;
t = linspace(t0, tf, N+1);

MOmat = MO(0, dt, K_theta, K_psi);

for i = 1:N-1
  MOmat = MO(i*dt,dt, K_theta, K_psi)*MOmat;   
end

y_at_pi_with_transition_matrix = MOmat*[1;1;1;1];

% Compute eigenvalues
eig(MOmat);

eigenvalues= eig(MOmat);
EigVals{kk1,kk2} = eigenvalues;

% Flutter 
for k = 1:length(eigenvalues)
    if any(abs(eigenvalues(k)) > 1)
       unstable = [unstable; K_psi, K_theta];
    end
    
% 1/Ω *(Divergence) proximity check
      if any(abs(abs(eigenvalues(k)) - 2/Omega) < 0.99887)
         Rdivergence_map = [Rdivergence_map; K_psi, K_theta];
      end 
end

 
            end 
        end


% Plot the Flutter and divergence maps

figure;
hold on;

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

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

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

hold off;


function Mat_MO = MO(t,dt, K_theta, K_psi)
  n = 4;

  At = A(t, dt, K_theta, K_psi);
  Atpdth = A(t+dt/2, dt,  K_theta, K_psi);
  Atpdt = A(t+dt, dt,  K_theta, K_psi);

  E = Atpdth*(eye(n)+dt/2*At);
  F = Atpdth*(eye(n)+(-1/2+1/sqrt(2))*dt*At+(1-1/sqrt(2))*dt*E);
  G = Atpdt*(eye(n)-dt/sqrt(2)*E+(1+1/sqrt(2))*dt*F);

  MO = eye(n)+dt/6*(At+2*(1-1/sqrt(2))*E+2*(1+1/sqrt(2))*F+G);
  Mat_MO = MO;
end
   
function Mat_A = A(t, dt, K_theta, K_psi)

global I_thetaoverI_b I_psioverI_b h gamma H_u M_b M_u V Omega omega_rad_s

    ct = cos(2*omega_rad_s^2*t);
    st = sin(2*omega_rad_s^2*t);

    % Define inertia matrix [M]
    M11 = I_thetaoverI_b + 1 + ct;
    M12 = -st;
    M21 = -st;
    M22 = I_psioverI_b + 1 - ct;

    M  = [M11 M12; M21 M22];
    
    % Define damping matrix [D]
    D11 =  h^2*gamma*H_u*(1 - ct) - gamma*M_b*(1 + ct) - (2 + 2*h*gamma*M_u)*st;
    D12 =  h^2*gamma*H_u*st + gamma*M_b*st - 2*(1 + ct) - 2*h*gamma*M_u*ct;
    D21 =  h^2*gamma*H_u*st + gamma*M_b*st + 2*(1 - ct) - 2*h*gamma*M_u*ct;
    D22 =  h^2*gamma*H_u*(1 + ct) - gamma*M_b*(1 - ct) + (2 + 2*h*gamma*M_u)*st;

    D = [D11 D12; D21 D22];

    % Define stiffness matrix [K]
    K11 = K_theta - h*gamma*H_u*(1 - ct) + gamma*M_u*st;
    K12 = -h*gamma*H_u*st + gamma*M_u*(1 + ct);
    K21 = -h*gamma*H_u*st - gamma*M_u*(1 - ct);
    K22 = K_psi - h*gamma*H_u*(1 + ct) - gamma*M_u*st;

    K = [K11 K12; K21 K22];

   % Compute the terms -M \ K and -M \ C 
   MK =-M \K; % Solves M * y = K for y
   MD =-M \D; % Solves M * y = D for y

   % Compute A(t)
   A = [ 0, 1, 0, 0;
         MK(1,1), MD(1,1), MK(1,2), MD(1,2);
         0, 0, 0, 1;
         MK(2,1), MD(2,1), MK(2,2), MD(2,2)];

   Mat_A = A;

end

The system has been solved, and the eigenvalues were found using the Floquet method. I have some doubts about applying the Runge-Kutta computational method.

Here is the link to the only reference I used:

https://ntrs.nasa.gov/citations/19740019387

Could you please review my code and confirm if there are no issues with the coding?

A 1/rev divergence also appears where one eigenvalue would be expected to be near 1/rev.(which is 2/Ω for this problem with a period of π)