Моделирование системы Лоренца, нахождение корреляции

clear all;
clc;
del_t=.0001; n=10^6;
sig=10; b=8/3;r=28;
x(1)=1.442; y(1)=2.582; z(1)=17.444;
t=(1:n)*del_t;
for i=2:n
x(i)=x(i-1)+sig*(y(i-1)-x(i-1))*del_t;
y(i)=y(i-1)+((r-z(i-1))*x(i-1)-y(i-1))*del_t;
z(i)=z(i-1)+(x(i-1)*y(i-1)-b*z(i-1))*del_t;
end
m1=mean(x); D1=var(x); cko1=sqrt(var(x));
subplot(3,1,1), plot(t,x)
grid on
text(2,15,['Expected value = ' num2str(m1) ',']);
text(35,15,['Variance = ' num2str(D1) ',']);
text(60,15,['Standard deviation = ' num2str(cko1) '.']);
title('Dynamic component X')
xlabel('Time')
ylabel('Amplitude')
m2=mean(y); D2=var(y); cko2=sqrt(var(y));
subplot(3,1,2), plot(t,y)
grid on
text(2,30,['Expected value = ' num2str(m2) ',']);
text(35,30,['Variance = ' num2str(D2) ',']);
text(60,30,['Standard deviation = ' num2str(cko2) '.']);
title('Dynamic component Y')
xlabel('Time')
ylabel('Amplitude')
m3=mean(z); D3=var(z); cko3=sqrt(var(z));
subplot(3,1,3), plot(t,z)
grid on
text(2,50,['Expected value = ' num2str(m3) ',']);
text(35,50,['Variance = ' num2str(D3) ',']);
text(60,50,['Standard deviation = ' num2str(cko3) '.']);
title('Dynamic component Z')
xlabel('Time')
ylabel('Amplitude')
figure
plot3 (x, y, z)
grid on