Моделирование пружинного маятника в MATLAB

clear, clc
global k m r
k = 200; % жесткость пружины
m = 0.5; % масса
r = 0.8; % коэффициент сопротивления
tk = [0 20]; % интервал времени наблюдения
x0 = [0.2, 0]; % начальные условия (координата и скорость)
 
%% исходное решение при типовых значениях параметров
[t,x] = ode45(@dsys,tk,x0); % решаем ДУ
figure(1)
plot(t, x(:,1))
grid on
title('График колебаний маятника')
xlabel('Время, с')
ylabel('Координата х, м')
xlim([0 8]) % для наглядности
 
%% 1. Определить зависимость времени затухания колебаний маятника от жёсткости пружины.
k_list = 10:50:500; % набор значений жесткости
tau = zeros(size(k_list)); % заготовка из нулей
for i = 1:length(k_list) 
    k = k_list(i); 
    [t,x] = ode45(@dsys,tk,x0);
    % определяем время затухания по уровню 5%:
    A = abs( x(:,1) ); % значения по модулю
    mA = max(A); % максимальное значение
    % находим один последний элемент, который больше 5% от максимума:
    id = find( A>=0.05*mA, 1, 'last' );
    tau(i) = t(id); % сохраняем в массив
end
figure(2)
plot(k_list, tau,'.-')
grid on
title({'Зависимость времени затухания колебаний', 'от жёсткости пружины'})
xlabel('Жесткость, k') 
ylabel('Время затухания \tau, с')
ylim([0 5]) % для наглядности
 
%% 2. Определить зависимость времени затухания колебаний маятника от его массы.
k = 200; % возвращаем типовое значение
m_list = 0.05:0.1:2; % набор значений массы
tau = zeros(size(m_list));
for i = 1:length(m_list)
    m = m_list(i);
    [t,x] = ode45(@dsys,tk,x0);
    % определяем время затухания по уровню 5%:
    A = abs( x(:,1) );
    mA = max(A);
    id = find( A>=0.05*mA, 1, 'last' );
    tau(i) = t(id);
end
figure(3)
plot(m_list, tau,'.-')
grid on
title({'Зависимость времени затухания колебаний', 'от его массы'})
xlabel('Масса, m') 
ylabel('Время затухания \tau, с')
 
%% 3. Определить зависимость времени затухания колебаний маятника от начального отклонения.
m = 0.5;
x0_list = 0.05:0.1:2; % набор значений начального отклонения
tau = zeros(size(x0_list));
for i = 1:length(x0_list)
    [t,x] = ode45(@dsys,tk,[x0_list(i), 0]);
    % определяем время затухания по уровню 5%:
    A = abs( x(:,1) );
    mA = max(A);
    id = find( A>=0.05*mA, 1, 'last' );
    tau(i) = t(id);
end
figure(4)
plot(x0_list, tau,'.-')
grid on
title({'Зависимость времени затухания колебаний', 'от начального отклонения'})
xlabel('Начальное отклонение х0, м') 
ylabel('Время затухания \tau, с')
ylim([0 5])
 
%% 4. Определить зависимость периода колебания маятника от коэффициента сопротивления.
r_list = 0.1:0.5:50; % набор значений коэффициента сопротивления
T = zeros(size(r_list));
for i = 1:length(r_list)
    r = r_list(i);
    [t,x] = ode45(@dsys,tk,x0);
    % определяем период колебаний:
    A = x(:,1); % значения координаты
    % находим 2 самых крупных пика, ставим ограничение на высоту пиков
    % 10^-5, чтобы не мешали шумы
    [pks,locs] = findpeaks(A,t,'NPeaks',2,'SortStr','descend','MinPeakHeight',10^-5);
    if length(locs)==2 % если два максимума
        T(i) = locs(2)-locs(1); % сохраняем время между ними
    else % если пик один, значит колебания кончились
        T(i:end) = nan;
        disp('Колебания отсутствуют, когда коэффициент сопротивления более:')
        disp(r) % выводим значение
        break % выходим из цикла
    end
end
figure(5)
plot(r_list, T,'.-')
grid on
title({'Зависимость периода колебаний', 'от коэффициента сопротивления'})
xlabel('Коэффициент сопротивления, r') 
ylabel('Период Т, с')
ylim([0, 0.5])
 
 
%% 5. Нарисовать фазовую диаграмму для исходных значений.
r = 0.8; % коэффициент сопротивления (как было в начале)
[t,x] = ode45(@dsys,tk,x0); 
figure(6)
plot(x(:,1), x(:,2))
grid on
title('Фазовая диаграмма')
xlabel('Координата х, м') 
ylabel('Скорость v, м/с')
 
 
 
% % % фнукиця системы % % %
function dx = dsys(t,x)
global k m r
v = x(2); % скорость (v=x')
beta = r/(2*m); % коэффициент затухания
w0 = sqrt(k/m); % собственная частота
dx(1,1) = v; % x' = v
dx(2,1) = -2*beta*v - w0^2*x(1); % x'' = -2*beta*x' - w0^2*x
end