florenz_xyz.m

一个小小的c语言编程

function florenz_xyz(htext_speed,str_speed,time,htext_sigma,str_sigma,sigma,htext_b,str_b,b,htext_r,str_r,r,X0,Y0,Z0,t0,tf)
global X1 X2 X3 t %解微分方程的函数
cla,set(htext_sigma,'string',[str_sigma,sprintf('%1.4g\',sigma)]);
cla,set(htext_b,'string',[str_b,sprintf('%1.4g\',b)]);
cla,set(htext_r,'string',[str_r,sprintf('%1.10f\',r)]);
cla,set(htext_speed,'string',[str_speed,sprintf('%1.4g\',time)]);

tspan=[t0,tf];
f0=[X0;Y0;Z0];
[t,X]=ode45('DYDt2',tspan,f0,[],sigma,b,r);
X1=X(:,1);X2=X(:,2);X3=X(:,3);
x1=-20;x2=20;y1=-30;y2=30;z1=0;z2=60;

singularity1_x=0;singularity1_y=0;singularity1_z=0;
singularity2_x=(b*(r-1))^(1/2);singularity2_y=(b*(r-1))^(1/2);singularity2_z=r-1;
singularity3_x=-(b*(r-1))^(1/2);singularity3_y=-(b*(r-1))^(1/2);singularity3_z=r-1;

figure(1)
plot3(X1,X2,X3,'W'),hold on

if r<=1
    plot3(singularity1_x,singularity1_y,singularity1_z,'k.','MarkerSize',20)       %画奇点
    text(singularity1_x,singularity1_y,singularity1_z,' O')
      else 
   plot3(singularity1_x,singularity1_y,singularity1_z,'k.','MarkerSize',20)       %画奇点
   text(singularity1_x,singularity1_y,singularity1_z,' O')
   plot3(singularity2_x,singularity2_y,singularity2_z,'k.','MarkerSize',20)       %画奇点
   text(singularity2_x,singularity2_y,singularity2_z,' C+')
   plot3(singularity3_x,singularity3_y,singularity3_z,'k.','MarkerSize',20)       %画奇点
   text(singularity3_x,singularity3_y,singularity3_z,' C-')
end
temp_time=time/400;
mycomet3(X1,X2,X3,0.02,temp_time);hold on,
title('Lorenz phase'),xlabel('x'),ylabel('y'),zlabel('z'),
%legend('z',4)    %选4可使彗星图不被刷掉                 
hold off