%% Programme MATLAB
clc; clear; close all;
%% Paramètres du système
m1 = 1.0; m2 = 1.0; m3 = 1.0; % Masses (kg)
k01 = 10; k12 = 10; k23 = 10; k34 = 10; k13 = 5; % Raideurs (N/m)
c12 = 0.5; c23 = 0.5; % Coefficients d'amortissement (Ns/m)
%% Matrices du système
% Matrice de masse
M = diag([m1, m2, m3]);
% Matrice de rigidité
K = -[ (k01 + k12 + k13), -k12, -k13;
-k12, (k12 + k23), -k23;
-k13, -k23, (k23 + k34 + k13) ];
% Matrice d'amortissement
C = [ c12, -c12, 0 ;
-c12, (c12 + c23), -c23;
0, -c23, c23 ];
%% Calcul des fréquences propres (non amorties)
[V,D] = eig(K,M);
omega_n = sqrt(diag(D)); % Fréquences propres en rad/s
f_n = omega_n/(2*pi); % Fréquences propres en Hz
% Tri par ordre croissant
[omega_n, idx] = sort(omega_n);
f_n = f_n(idx);
V = V(:,idx);
%% Affichage des fréquences propres dans la fenêtre de commande
disp('=== Fréquences propres du système ===');
for i = 1:length(omega_n)
fprintf('Mode %d: ?_%d = %.4f rad/s (f_%d = %.4f Hz)\n',...
i, i, omega_n(i), i, f_n(i));
end
disp(' ');
%% Calcul des fréquences de résonance (amorties)
% On prend en compte l'amortissement pour trouver les pics de résonance
Z = zeros(3);
I = eye(3);
A = [ Z, I;
-M\K, -M\C ];
% Calcul des valeurs propres du système amorti
lambda = eig(A);
lambda = lambda(imag(lambda) > 1e-6 ); % On garde seulement les pôles à partie
imaginaire positive
if isempty(lambda)
warning('Aucune fréquence amortie trouvée - utilisation des fréquences propres');
omega_d = omega_n;
f_d = f_n;
else
omega_d = abs(imag(lambda));
f_d = omega_d/(2*pi);
[omega_d, idx] = unique(sort(omega_d));
f_d = f_d(idx);
end
if length(f_d) < 2
f_d(end+1) = f_d(end)*1.1;
end
%% Affichage des fréquences
disp('=== Fréquences propres ===');
disp(table(omega_n, f_n, 'RowNames', {'Mode 1','Mode 2','Mode 3'}));
disp('=== Fréquences de résonance ===');
if length(f_d) >= 3
disp(table(omega_d(1:3), f_d(1:3), 'RowNames', {'Mode 1','Mode 2','Mode 3'}));
else
disp(table(omega_d, f_d, 'RowNames', compose('Mode %d',1:length(f_d))));
end
disp('=== Fréquences de résonance (amorties) ===');
for i = 1:length(omega_d)
fprintf('Pic de résonance %d: ?_%d = %.4f rad/s (f_%d = %.4f Hz)\n',...
i, i, omega_d(i), i, f_d(i));
end
disp(' ');
%% Figures pour les fréquences
figure;
subplot(2,1,1);
bar(omega_n);
title('Fréquences propres non amorties (?_n)');
xlabel('Mode');
ylabel('Fréquence (rad/s)');
xticks(1:3);
grid on;
subplot(2,1,2);
bar(omega_d);
title('Fréquences de résonance amorties (?_d)');
xlabel('Mode');
ylabel('Fréquence (rad/s)');
xticks(1:3);
grid on;
%% Conditions initiales
x0 = [0; 0.1; 0]; % Positions initiales (m)
v0 = [0; 0; 0]; % Vitesses initiales (m/s)
%% Force externe avec vérification
if length(f_d) >= 2
freq_force = f_d(2); % On excite près de la 2ème fréquence de
résonance
else
freq_force = f_n(2);
warning('Utilisation fréquence propre comme fréquence de forçage');
end
F_amplitude = 1.0; % Amplitude de la force (N)
F = @(t) [sin(2*pi*freq_force*t);
cos(2*pi*freq_force*t);
0.5*sin(2*pi*freq_force*t)];
%% Paramètres de simulation
t_start = 0;
t_end = 20; % Temps final (s)
h = 0.01; % Pas de temps (s)
steps = ceil((t_end - t_start)/h);
%% Conversion en système du premier ordre
B = [ Z;
M\eye(3) ];
%% Initialisation
t = linspace(t_start, t_end, steps)';
X_euler = zeros(2*3, steps);
X_rk4 = zeros(2*3, steps);
X_euler(:,1) = [x0; v0];
X_rk4(:,1) = [x0; v0];
%% Méthode d'Euler explicite
for i = 1:steps-1
X_euler(:,i+1) = X_euler(:,i) + h*(A*X_euler(:,i) + B*F(t(i)));
end
%% Méthode de Runge-Kutta d'ordre 4
for i = 1:steps-1
k1 = A*X_rk4(:,i) + B*F(t(i));
k2 = A*(X_rk4(:,i) + h/2*k1) + B*F(t(i)+h/2);
k3 = A*(X_rk4(:,i) + h/2*k2) + B*F(t(i)+h/2);
k4 = A*(X_rk4(:,i) + h*k3) + B*F(t(i)+h);
X_rk4(:,i+1) = X_rk4(:,i) + h/6*(k1 + 2*k2 + 2*k3 + k4);
end
%% Extraction des résultats
x_euler = X_euler(1:3,:);
v_euler = X_euler(4:6,:);
x_rk4 = X_rk4(1:3,:);
v_rk4 = X_rk4(4:6,:);
%% Calcul des accélérations
% Pour Euler
a_euler = zeros(3,steps);
% Pour Runge-Kutta
a_rk4 = zeros(3,steps);
for i = 1:steps
a_euler(:,i) = M\(F(t(i)) - K*x_euler(:,i) - C*v_euler(:,i));
a_rk4(:,i) = M\(F(t(i)) - K*x_rk4(:,i) - C*v_rk4(:,i));
end
%% Affichage des résultats dynamiques complets
figure('Name','Comparaison Euler vs RK4 - Positions','Position',[100 100 1200 800]);
for mass = 1:3
subplot(3,1,mass);
plot(t, x_euler(mass,:), 'b-', t, x_rk4(mass,:), 'r--');
title(['Position de la masse ' num2str(mass)]);
legend('Euler', 'RK4', 'Location','best');
xlabel('Temps (s)');
ylabel('Position (m)');
grid on;
end
figure('Name','Comparaison Euler vs RK4 - Vitesses','Position',[100 100 1200 800]);
for mass = 1:3
subplot(3,1,mass);
plot(t, v_euler(mass,:), 'b-', t, v_rk4(mass,:), 'r--');
title(['Vitesse de la masse ' num2str(mass)]);
legend('Euler', 'RK4', 'Location','best');
xlabel('Temps (s)');
ylabel('Vitesse (m/s)');
grid on;
end
figure('Name','Comparaison Euler vs RK4 - Accélérations','Position',[100 100 1200 800]);
for mass = 1:3
subplot(3,1,mass);
plot(t, a_euler(mass,:), 'b-', t, a_rk4(mass,:), 'r--');
title(['Accélération de la masse ' num2str(mass)]);
legend('Euler', 'RK4', 'Location','best');
xlabel('Temps (s)');
ylabel('Accélération (m/s²)');
grid on;
end
%% Calcul des erreurs entre les méthodes
erreur_pos = abs(x_euler - x_rk4);
erreur_vit = abs(v_euler - v_rk4);
erreur_acc = abs(a_euler - a_rk4);
%% Affichage des erreurs maximales
disp('=== Erreurs maximales entre Euler et RK4 ===');
for mass = 1:3
fprintf('Masse %d:\n', mass);
fprintf(' Position: %.4e m\n', max(erreur_pos(mass,:)));
fprintf(' Vitesse: %.4e m/s\n', max(erreur_vit(mass,:)));
fprintf(' Accélération: %.4e m/s²\n\n', max(erreur_acc(mass,:)));
end
%% Affichage des RMS (Root Mean Square) des erreurs
disp('=== RMS des erreurs entre Euler et RK4 ===');
for mass = 1:3
fprintf('Masse %d:\n', mass);
fprintf(' Position: %.4e m\n', rms(erreur_pos(mass,:)));
fprintf(' Vitesse: %.4e m/s\n', rms(erreur_vit(mass,:)));
fprintf(' Accélération: %.4e m/s²\n\n', rms(erreur_acc(mass,:)));
end
%% Animation comparative (optionnelle)
figure;
for i = 1:10:steps
clf;
subplot(1,2,1);
plot([0 x_euler(1,i) x_euler(2,i) x_euler(3,i)], [0 0 0 0], 'bo-');
title('Méthode Euler');
xlim([-0.5 0.5]); ylim([-0.1 0.1]);
subplot(1,2,2);
plot([0 x_rk4(1,i) x_rk4(2,i) x_rk4(3,i)], [0 0 0 0], 'ro-');
title('Méthode Runge-Kutta');
xlim([-0.5 0.5]); ylim([-0.1 0.1]);
sgtitle(sprintf('Temps: %.2f s', t(i)));
drawnow;
end
%% Analyse spectrale (FFT) pour comparer les fréquences
figure('Name','Analyse spectrale des positions','Position',[100 100 1200 800]);
for mass = 1:3
subplot(3,1,mass);
% Calcul de la FFT
Fs = 1/h; % Fréquence d'échantillonnage
L = length(t);
f = Fs*(0:(L/2))/L;
Y_euler = fft(x_euler(mass,:));
P2_euler = abs(Y_euler/L);
P1_euler = P2_euler(1:L/2+1);
P1_euler(2:end-1) = 2*P1_euler(2:end-1);
Y_rk4 = fft(x_rk4(mass,:));
P2_rk4 = abs(Y_rk4/L);
P1_rk4 = P2_rk4(1:L/2+1);
P1_rk4(2:end-1) = 2*P1_rk4(2:end-1);
plot(f, P1_euler, 'b-', f, P1_rk4, 'r--');
title(['Spectre de la masse ' num2str(mass)]);
legend('Euler', 'RK4', 'Location','best');
xlabel('Fréquence (Hz)');
ylabel('Amplitude');
xlim([0 2]);
grid on;
end
%% Affichage des performances computationnelles
disp('=== Performances computationnelles ===');
tic;
% Test Euler
X_test = [x0; v0];
for i = 1:steps-1
X_test = X_test + h*(A*X_test + B*F(t(i)));
end
t_euler = toc;
tic;
% Test RK4
X_test = [x0; v0];
for i = 1:steps-1
k1 = A*X_test + B*F(t(i));
k2 = A*(X_test + h/2*k1) + B*F(t(i)+h/2);
k3 = A*(X_test + h/2*k2) + B*F(t(i)+h/2);
k4 = A*(X_test + h*k3) + B*F(t(i)+h);
X_test = X_test + h/6*(k1 + 2*k2 + 2*k3 + k4);
end
t_rk4 = toc;
fprintf('Temps de calcul Euler: %.4f s\n', t_euler);
fprintf('Temps de calcul RK4: %.4f s\n', t_rk4);
fprintf('Ratio RK4/Euler: %.2f\n', t_rk4/t_euler);
%% Calcul de l'énergie totale (cinétique + potentielle)
% Énergie potentielle élastique
Epe_euler = zeros(1,steps);
Epe_rk4 = zeros(1,steps);
for i = 1:steps
Epe_euler(i) = 0.5 * x_euler(:,i)' * K * x_euler(:,i);
Epe_rk4(i) = 0.5 * x_rk4(:,i)' * K * x_rk4(:,i);
end
% Énergie cinétique
Eke_euler = 0.5 * sum(v_euler.^2 .* diag(M), 1);
Eke_rk4 = 0.5 * sum(v_rk4.^2 .* diag(M), 1);
% Énergie totale
Etot_euler = Epe_euler + Eke_euler;
Etot_rk4 = Epe_rk4 + Eke_rk4;
%% Affichage des résultats énergétiques
disp('=== Énergie du système ===');
fprintf('Énergie initiale (Euler): %.4f J\n', Etot_euler(1));
fprintf('Énergie initiale (RK4): %.4f J\n', Etot_rk4(1));
fprintf('Énergie moyenne (Euler): %.4f J\n', mean(Etot_euler));
fprintf('Énergie moyenne (RK4): %.4f J\n', mean(Etot_rk4));
fprintf('Perte d''énergie relative (Euler): %.2f%%\n',...
100*(Etot_euler(1)-Etot_euler(end))/Etot_euler(1));
fprintf('Perte d''énergie relative (RK4): %.2f%%\n',...
100*(Etot_rk4(1)-Etot_rk4(end))/Etot_rk4(1));
disp(' ');
%% Figures pour les résultats temporels
figure;
% Positions
subplot(3,1,1);
plot(t, x_euler(1,:), 'b-', t, x_rk4(1,:), 'r--');
title(['Position de la masse 1 - Excitation à ' num2str(freq_force) ' Hz']);
legend('Euler', 'RK4');
xlabel('Temps (s)');
ylabel('Position (m)');
grid on;
subplot(3,1,2);
plot(t, x_euler(2,:), 'b-', t, x_rk4(2,:), 'r--');
title('Position de la masse 2');
legend('Euler', 'RK4');
xlabel('Temps (s)');
ylabel('Position (m)');
grid on;
subplot(3,1,3);
plot(t, x_euler(3,:), 'b-', t, x_rk4(3,:), 'r--');
title('Position de la masse 3');
legend('Euler', 'RK4');
xlabel('Temps (s)');
ylabel('Position (m)');
grid on;
% Énergie
figure;
plot(t, Etot_euler, 'b-', t, Etot_rk4, 'r--');
title('Évolution de l''énergie totale');
legend('Euler', 'RK4');
xlabel('Temps (s)');
ylabel('Énergie (J)');
grid on;
% Positions - Comparaison Euler vs RK4
figure;
subplot(3,1,1);
plot(t, x_euler(1,:), 'r', t, x_rk(1,:), 'r--', 'LineWidth', 1.5);
hold on;
plot(t, x_euler(2,:), 'g', t, x_rk(2,:), 'g--', 'LineWidth', 1.5);
plot(t, x_euler(3,:), 'b', t, x_rk(3,:), 'b--', 'LineWidth', 1.5);
title('Positions des masses - Comparaison Euler et RK4');
xlabel('Temps (s)');
ylabel('Position (m)');
legend('m_1 Euler', 'm_1 RK4', 'm_2 Euler', 'm_2 RK4', 'm_3 Euler', 'm_3 RK4');
grid on;
% Vitesses - Comparaison Euler vs RK4
subplot(3,1,2);
plot(t, v_euler(1,:), 'r', t, v_rk(1,:), 'r--', 'LineWidth', 1.5);
hold on;
plot(t, v_euler(2,:), 'g', t, v_rk(2,:), 'g--', 'LineWidth', 1.5);
plot(t, v_euler(3,:), 'b', t, v_rk(3,:), 'b--', 'LineWidth', 1.5);
title('Vitesses des masses - Comparaison Euler et RK4');
xlabel('Temps (s)');
ylabel('Vitesse (m/s)');
grid on;
% Accélérations - Comparaison Euler vs RK4
subplot(3,1,3);
plot(t, a_euler(1,:), 'r', t, a_rk(1,:), 'r--', 'LineWidth', 1.5);
hold on;
plot(t, a_euler(2,:), 'g', t, a_rk(2,:), 'g--', 'LineWidth', 1.5);
plot(t, a_euler(3,:), 'b', t, a_rk(3,:), 'b--', 'LineWidth', 1.5);
title('Accélérations des masses - Comparaison Euler et RK4');
xlabel('Temps (s)');
ylabel('Accélération (m/s^2)');
grid on;
%% Comparaison des méthodes
erreur_pos = max(abs(x_euler - x_rk4), [], 'all');
fprintf('Erreur maximale entre Euler et RK4: %.4e m\n', erreur_pos);
