Скелетная линия колебательной системы с 3 степенями свободы

function dxdt = diffsys(t, x)
% обозначения:
% x(1) -> x(t) (искомая функция)
% x(2) -> dx/dt (производная)
% dxdt(1,1) -> dx/dt (первая производная)
% dxdt(2,1) -> d2x/dt2 (вторая производная)
dxdt(1,1) = x(2); 
dxdt(2,1) = x(1)-0.5*x(1).^3;
end

clear; clc; close all;
N=2^10;
T=2*pi;
T=32*T;
nfft=N;
t=0:T/N:T;
fs=1/(t(2)-t(1));
i_ot=0;
for i0=0.01:0.01:10
x0 = [i0 0];
[T0 X] = ode45('diffsys',t,x0);
X=X';
x=X(2,:);
S=fft(x,nfft)/nfft;
mag=2*abs(S);
mag(1)=mag(1)/2;
freg=fs/2*linspace(-1,1,nfft);
f_mag=fftshift(mag);
d_f=abs(freg(2)-freg(1))/2;
freg=freg-d_f;
s_freg=freg(nfft/2:end);
s_freg(2)=0;
s_mag=f_mag(nfft/2:end);
[pks locks]=findpeaks(s_mag);
i_ot=i_ot+1;
otX(i_ot)=s_freg(locks(1));
otY(1,i_ot)=pks(1);
otY(2,i_ot)=pks(2);
otY(3,i_ot)=pks(3);
end
plot (otX.*2*pi, otY([1 2 3],:))
legend('1я гармоника','2я гармоника','3я гармоника',3)

clear; clc
 
dsys = @(t,x) [x(2);
    x(1)-0.5*x(1)^3];
Tk = [0 100];
nfft = 2^10;
ti = linspace(0,Tk(end),nfft);
fs=1/(ti(2)-ti(1));
freq=fs/2*linspace(-1,1,nfft);
d_f=abs(freq(2)-freq(1))/2;
freq=freq-d_f;
s_freq=freq(nfft/2:end);
s_freq(2)=0;
 
x0 = 0.1:0.1:1;
 
[F X0] = meshgrid(s_freq, x0);
S = zeros(size(F));
for i = 1:length(x0)
    [t X] = ode45(dsys,Tk,[x0(i),0]);
    x = X(:,2);
    xi = spline(t,x, ti);
    s=fft(xi,nfft)/nfft;
    mag=2*abs(s);
    mag(1)=mag(1)/2;
    f_mag=fftshift(mag);
    s_mag=f_mag(nfft/2:end);
    S(i,:) = s_mag;
end
 
surf(F,X0,S,'FaceAlpha',0.5,'MeshStyle','row')
 
view(30,30)
xlim([0 1])
xlabel('частота')
ylabel('нач. условия')
title(['Длительность Т=',num2str(Tk(end))])

clear; clc
 
dsys = @(t,x) [x(2); 
                -sin(x(1))];
Tk = [0 64*pi];
count=0;
iN=0.01;
iK=5;
iH=0.01;
for i0=iN:iH:iK
    count=count+1;
x0 = [i0, 0];
[t X] = ode45(dsys,Tk,x0);
x = X(:,2);
 
nfft = 2^10;
ti = linspace(0,t(end),nfft);
xi = spline(t,x, ti);
 
 
S=fft(xi,nfft)/nfft;
mag=2*abs(S);
mag(1)=mag(1)/2;
f_mag=fftshift(mag);
s_mag=f_mag(nfft/2:end);
 
fs=1/(ti(2)-ti(1));
freq=fs/2*linspace(-1,1,nfft);
d_f=abs(freq(2)-freq(1))/2;
freq=freq-d_f;
s_freq=freq(nfft/2:end);
s_freq(2)=0;
[pks locks]=findpeaks(s_mag);
xFin(count)=s_freq(locks(1));
yFin(1,count)=pks(1);
yFin(2,count)=pks(2);
yFin(3,count)=pks(3);
end
 
subplot(211)
plot(t,x)
title(['x0=',num2str(iN),':',num2str(iH),':',num2str(iK),' , T=',num2str(t(end)/pi),'\pi'])
 
subplot(212)
plot(xFin*2*pi,yFin)

[pks locks]=findpeaks(s_mag);
xFin(count)=s_freq(locks(1));
yFin(1,count)=pks(1);
yFin(2,count)=pks(2);
yFin(3,count)=pks(3);