Lorenz_attractor

% 洛伦兹的蝴蝶
%% 洛伦兹方程参数与计算
r=28;s=10;b=8/3;
% 初始数据点数，演示时间
N=16;ti=0:0.005:108;
% ode设置
opts=odeset('RelTol',1e8,'AbsTol',1e10,'NormControl','on',...
    'MaxStep',1e-2,'InitialStep',1e-3);
for ii=1:N
[~,u] = ode45(@(t,u)Lorenz(t,u,r,s,b),ti,[0 ii/N 0],opts);
x(:,ii)=u(:,1);y(:,ii)=u(:,2);z(:,ii)=u(:,3);
end
%% 绘图
f=figure;f.MenuBar='none';
dpi=1.25;scr=[1920 1080];res=[1280 960];
fig=gcf;fig.Color='k';
fig.Position = [(scr-res)/2,res]/dpi;% 居中显示
cj=jet;% 颜色集
for ii=1:N
    a(ii)=animatedline('MaximumNumPoints',1250, ...
        'Color',cj(ii*16,:),'LineWidth',1.5);
end
axis off;axis equal
ax=gca;ax.Clipping='off';
ax.CameraViewAngle = 22.5;%视野范围
ax.CameraTarget = [0 0 25];%目标位置
%% 动画
mv = 1;% 是否录制视频，1录制，其它不录制
if mv == 1
v = VideoWriter('D:\Programs\MATLAB\美丽的蝴蝶.mp4','MPEG-4');
v.Quality = 100;v.FrameRate = 60;
open(v);
end
for n=1:length(ti)
    for ii=1:N
        addpoints(a(ii),x(n,ii),y(n,ii),z(n,ii));
    end
    ax.CameraPosition = [100*cos(ti(n)/2) 100*sin(ti(n)/2) 25];%相机位置
    drawnow% limitrate
    if mv==1;F=getframe(fig);writeVideo(v,F);end
end
if mv == 1;close(v);end
%% 洛伦兹方程
function du=Lorenz(~,u,r,s,b)
x=u(1);y=u(2);z=u(3);
dx=s*(y-x);
dy=r*x-y-x*z;
dz=-b*z+x*y;
du=[dx dy dz]';
end