testlorenz.m

来自「lorenz吸引子的画法」· M 代码 · 共 30 行

clear;
clc;
X0=[8;9;10];
X=X0;
z=X0;
Tf=60;T0=0;h=0.002;t=T0;
N=round((Tf-T0)/h)-1;
m=round((Tf/2-T0)/h)-1;
 for i=1:m
           K1=lorenz(X);
    K2=lorenz(X'+K1*h/2);
    K3=lorenz(X'+K2*h/2);
    K4=lorenz(X'+K3*h);
    X=X+h*(K1'+2*K2'+2*K3'+K4')/6;
    z(:,i)=X;
    t=[t,t(i)+h];
 end
for i=m:N
           K1=lorenz1(X);
    K2=lorenz1(X'+K1*h/2);
    K3=lorenz1(X'+K2*h/2);
    K4=lorenz1(X'+K3*h);
    X=X+h*(K1'+2*K2'+2*K3'+K4')/6;
    z(:,i)=X;
    t=[t,t(i)+h];
end
figure(1),plot3(z(1,:),z(2,:),z(3,:)),grid on;
figure(2),subplot(3,1,1),plot(t(1:length(z(1,:))),z(1,:)),grid on;
subplot(3,1,2),plot(t(1:length(z(2,:))),z(2,:)),grid on;
subplot(3,1,3),plot(t(1:length(z(3,:))),z(3,:)),grid on;