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

148 views (last 30 days)
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
  2 Comments
Torsten
Torsten on 13 Nov 2024 at 20:41
Edited: Torsten on 13 Nov 2024 at 20:41
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
Nikoo on 13 Nov 2024 at 21:04
k1 = Function_system(t(i), y(:, i));
k2 = Function_system(t(i) + ts/2, y(:, i) + ts/2 * k1);
All should have 2

Sign in to comment.

Answers (2)

Torsten
Torsten on 13 Nov 2024 at 21:11
Edited: Torsten on 13 Nov 2024 at 21:17
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
  13 Comments
Nikoo
Nikoo on 18 Nov 2024 at 21:40
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
Torsten on 18 Nov 2024 at 22:33
Edited: Torsten on 18 Nov 2024 at 22:34
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" ?

Sign in to comment.


Torsten
Torsten on 22 Nov 2024 at 16:04
Edited: Torsten on 22 Nov 2024 at 16:15
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>
  6 Comments
Nikoo
Nikoo on 29 Nov 2024 at 4:05
Edited: Nikoo on 6 Dec 2024 at 19:52
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 π)

Sign in to comment.

Categories

Find more on Programming in Help Center and File Exchange

Tags

Products


Release

R2023a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!