How can I evaluate the transition matrix using the fourth order Runge-Kutta ODEs
148 views (last 30 days)
Show older comments
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
Answers (2)
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
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" ?
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]
% Compute eigenvalues
eig(Knum)
% 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,:).'
6 Comments
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!