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.