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

Question

Nikoo 2024-11-13，20:29

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2166119-how-can-i-evaluate-the-transition-matrix-using-the-fourth-order-runge-kutta-odes

编辑： Torsten about 1 hour 前

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 个评论
显示无隐藏无

Torsten 2024-11-13，20:41

编辑：Torsten 2024-11-13，20:41

在 MATLAB Online 中打开

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 2024-11-13，21:04

在 MATLAB Online 中打开

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 2024-11-13，21:11

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2166119-how-can-i-evaluate-the-transition-matrix-using-the-fourth-order-runge-kutta-odes#answer_1544939

编辑：Torsten 2024-11-13，21:17

在 MATLAB Online 中打开

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 个评论
显示 11更早的评论隐藏 11更早的评论

Nikoo 2024-11-14，4:03

在 MATLAB Online 中打开

% 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?

Torsten 2024-11-14，10:24

编辑：Torsten 2024-11-14，10:47

在 MATLAB Online 中打开

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)

Nikoo 2024-11-17，23:22

在 MATLAB Online 中打开

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 2024-11-18，2:34

编辑：Torsten 2024-11-18，2:37

在 MATLAB Online 中打开

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 2024-11-18，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 2024-11-18，22:33

编辑：Torsten 2024-11-18，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" ?

请先登录，再进行评论。

Answer 2

Torsten about 1 hour 前

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2166119-how-can-i-evaluate-the-transition-matrix-using-the-fourth-order-runge-kutta-odes#answer_1548598

编辑：Torsten about 1 hour 前

在 MATLAB Online 中打开

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>

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

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

2 个评论
显示无隐藏无

回答（2 个）

13 个评论
显示 11更早的评论隐藏 11更早的评论

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

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

2 个评论 显示 无隐藏 无

回答（2 个）

13 个评论 显示 11更早的评论隐藏 11更早的评论

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

2 个评论
显示无隐藏无

13 个评论
显示 11更早的评论隐藏 11更早的评论

0 个评论
显示 -2更早的评论隐藏 -2更早的评论