unified.m

分数阶统一混沌系统的离散算法-This the program for fractional order Unified system, a simple but fast method in ti

%This is the program for fractional order Unified system, a simple but fast
%method in time domain.
clc
clear
clf
a=1; %parameter a
q1=0.99;q2=0.99;q3=0.99; %the fractional order
x10=3;x20=1;x30=5;
h=0.01;N=1000;
x1(N+1)=[0];x2(N+1)=[0];x3(N+1)=[0];
xx1(N+1)=[0];xx2(N+1)=[0];xx3(N+1)=[0];
M1=0;M2=0;M3=0;N1=0;N2=0;N3=0;

xx1(1)=x10+h^q1*(25*a+10)*(x20-x10)/(q1*gamma(q1));
xx2(1)=x20+h^q2*((28-35*a-x30)*x10+(29*a-1)*x20)/(q2*gamma(q2));
xx3(1)=x30+h^q3*(x10*x20-x30*(8+a)/3)/(q3*gamma(q3));

x1(1)=x10+h^q1*(25*a+10)*(xx2(1)-xx1(1)+q1*(x20-x10))/gamma(q1+2);
x2(1)=x20+h^q2*((28-35*a-xx3(1))*xx1(1)+(29*a-1)*xx2(1)+q2*((28-35*a-x30)*x10+(29*a-1)*x20))/gamma(q2+2);
x3(1)=x30+h^q3*(xx1(1)*xx2(1)-xx3(1)*(8+a)/3+q3*(x10*x20-x30*(8+a)/3))/gamma(q3+2);

for n=1:N
    M1=(n^(q1+1)-(n-q1)*(n+1)^q1)*(x20-x10);
    M2=(n^(q2+1)-(n-q2)*(n+1)^q2)*((28-35*a-x30)*x10+(29*a-1)*x20);
    M3=(n^(q3+1)-(n-q3)*(n+1)^q3)*(x10*x20-x30*(8+a)/3);
    
    N1=((n+1)^q1-n^q1)*(x20-x10);
    N2=((n+1)^q2-n^q2)*((28-35*a-x30)*x10+(29*a-1)*x20);
    N3=((n+1)^q3-n^q3)*(x10*x20-x30*(8+a)/3);
    
    for j=1:n
        M1=M1+((n-j+2)^(q1+1)+(n-j)^(q1+1)-2*(n-j+1)^(q1+1))*(x2(j)-x1(j));
        M2=M2+((n-j+2)^(q2+1)+(n-j)^(q2+1)-2*(n-j+1)^(q2+1))*((28-35*a-x3(j))*x1(j)+(29*a-1)*x2(j));
        M3=M3+((n-j+2)^(q3+1)+(n-j)^(q3+1)-2*(n-j+1)^(q3+1))*(x1(j)*x2(j)-x3(j)*(8+a)/3);
        
        N1=N1+((n-j+1)^q1-(n-j)^q1)*(x2(j)-x1(j));
        N2=N2+((n-j+1)^q2-(n-j)^q2)*((28-35*a-x3(j))*x1(j)+(29*a-1)*x2(j));
        N3=N3+((n-j+1)^q3-(n-j)^q3)*(x1(j)*x2(j)-x3(j)*(8+a)/3);        
    end
    
    xx1(n+1)=x10+(25*a+10)*h^q1*N1/(q1*gamma(q1));
    xx2(n+1)=x20+h^q2*N2/(q2*gamma(q2));
    xx3(n+1)=x30+h^q3*N3/(q3*gamma(q3));
    
    x1(n+1)=x10+(25*a+10)*h^q1*(xx2(n+1)-xx1(n+1)+M1)/gamma(q1+2);
    x2(n+1)=x20+h^q2*((28-35*a-xx3(n+1))*xx1(n+1)+(29*a-1)*xx2(n+1)+M2)/gamma(q2+2);
    x3(n+1)=x30+h^q3*(xx1(n+1)*xx2(n+1)-xx3(n+1)*(8+a)/3+M3)/gamma(q3+2);   
    
end
x1=x1';x2=x2';x3=x3';

plot3(x3(100:end),x1(100:end),x2(100:end),'k');
xlabel('\itz','fontsize',18);
ylabel('\itx','fontsize',18);
zlabel('\ity','fontsize',18);
figure
plot(x1(100:end),x2(100:end),'k');
xlabel('\itx','fontsize',18);
ylabel('\ity','fontsize',18);
figure
plot(x1(100:end),x3(100:end),'k');
xlabel('\itx','fontsize',18);
ylabel('\itz','fontsize',18);
figure
plot(x2(100:end),x3(100:end),'k');
xlabel('\ity','fontsize',18);
ylabel('\itz','fontsize',18);