% 基于反馈控制的分数阶混沌系统同步
c=5;
% q1=0.98;q2=0.98;q3=0.98;
% q1=1;q2=1;q3=1;
q1=0.8;q2=0.8;q3=0.8;
h=0.01; % iteration step size
N=4000; % iteration number
kk=15;
x0=-6.7781;y0=-8.3059;z0=9.2133; % initial values of the system values;
x(N+1)=0;y(N+1)=0;z(N+1)=0;
x1(N+1)=0;y1(N+1)=0;z1(N+1)=0;
X0=1;Y0=2;Z0=3; % initial values of the sYstem values;
X(N+1)=0;Y(N+1)=0;Z(N+1)=0;
X1(N+1)=0;Y1(N+1)=0;Z1(N+1)=0;
x1(1)=x0+h^q1*10*(y0-x0)/(gamma(q1)*q1);
y1(1)=y0+h^q2*((24-4*c)*x0-x0*z0+c*y0)/(gamma(q2)*q2);
z1(1)=z0+h^q3*(x0*y0-8*z0/3)/(gamma(q3)*q3);
x(1)=x0+h^q1*(10*(y1(1)-x1(1))+q1*10*(y0-x0))/gamma(q1+2);
y(1)=y0+h^q2*((24-4*c)*x1(1)-x1(1)*z1(1)+c*y1(1)+q2*((24-4*c)*x0-x0*z0+c*y0))/gamma(q2+2);
z(1)=z0+h^q3*(x1(1)*y1(1)-8*z1(1)/3+q3*(x0*y0-8*z0/3))/gamma(q3+2);
X1(1)=X0+h^q1*(10*(Y0-X0)-kk*(X0-x0))/(gamma(q1)*q1);
Y1(1)=Y0+h^q2*((24-4*c)*X0-X0*Z0+c*Y0-kk*(Y0-y0))/(gamma(q2)*q2);
Z1(1)=Z0+h^q3*(X0*Y0-8*Z0/3-kk*(Z0-z0))/(gamma(q3)*q3);
X(1)=X0+h^q1*(10*(Y1(1)-X1(1)-kk*(X1(1)-x1(1)))+q1*(10*(Y0-X0)-kk*(X0-x0)))/gamma(q1+2);
Y(1)=Y0+h^q2*((24-4*c)*X1(1)-X1(1)*Z1(1)+c*Y1(1)-kk*(Y1(1)-y1(1))+q2*((24-4*c)*X0-X0*Z0+c*Y0-kk*(Y0-y0)))/gamma(q2+2);
Z(1)=Z0+h^q3*(X1(1)*Y1(1)-8*Z1(1)/3-kk*(Z1(1)-z1(1))+q3*(X0*Y0-8*Z0/3-kk*(Z0-z0)))/gamma(q3+2);
for n=1:N
% kk = 15;
n
% system 1
m1=(n^(q1+1)-(n-q1)*(n+1)^q1)*10*(y0-x0);
m2=(n^(q2+1)-(n-q2)*(n+1)^q2)*((24-4*c)*x0-x0*z0+c*y0);
m3=(n^(q3+1)-(n-q3)*(n+1)^q3)*(x0*y0-8*z0/3);
n1=((n+1)^q1-n^q1)*10*(y0-x0);
n2=((n+1)^q2-n^q2)*((24-4*c)*x0-x0*z0+c*y0);
n3=((n+1)^q3-n^q3)*(x0*y0-8*z0/3);
% system 2
M1=(n^(q1+1)-(n-q1)*(n+1)^q1)*(10*(Y0-X0)-kk*(X0-x0));
M2=(n^(q2+1)-(n-q2)*(n+1)^q2)*((24-4*c)*X0-X0*Z0+c*Y0-kk*(Y0-y0));
M3=(n^(q3+1)-(n-q3)*(n+1)^q3)*(X0*Y0-8*Z0/3-kk*(Z0-z0));
N1=((n+1)^q1-n^q1)*(10*(Y0-X0)-kk*(X0-x0));
N2=((n+1)^q2-n^q2)*((24-4*c)*X0-X0*Z0+c*Y0-kk*(Y0-y0));
N3=((n+1)^q3-n^q3)*(X0*Y0-8*Z0/3-kk*(Z0-z0));
for j=1:n
m1=m1+((n-j+2)^(q1+1)+(n-j)^(q1+1)-2*(n-j+1)^(q1+1))*(10*(y(j)-x(j)));
m2=m2+((n-j+2)^(q2+1)+(n-j)^(q2+1)-2*(n-j+1)^(q2+1))*((24-4*c)*x(j)-x(j)*z(j)+c*y(j));
m3=m3+((n-j+2)^(q3+1)+(n-j)^(q3+1)-2*(n-j+1)^(q3+1))*(x(j)*y(j)-8*z(j)/3);
n1=n1+((n-j+1)^q1-(n-j)^q1)*(10*(y(j)-x(j)));
n2=n2+((n-j+1)^q2-(n-j)^q2)*((24-4*c)*x(j)-x(j)*z(j)+c*y(j));
n3=n3+((n-j+1)^q3-(n-j)^q3)*(x(j)*y(j)-8*z(j)/3);
M1=M1+((n-j+2)^(q1+1)+(n-j)^(q1+1)-2*(n-j+1)^(q1+1))*(10*(Y(j)-X(j))-kk*(X(j)-x(j)));
M2=M2+((n-j+2)^(q2+1)+(n-j)^(q2+1)-2*(n-j+1)^(q2+1))*((24-4*c)*X(j)-X(j)*Z(j)+c*Y(j)-kk*(Y(j)-y(j)));
M3=M3+((n-j+2)^(q3+1)+(n-j)^(q3+1)-2*(n-j+1)^(q3+1))*(X(j)*Y(j)-8*Z(j)/3-kk*(Z(j)-z(j)));
N1=N1+((n-j+1)^q1-(n-j)^q1)*(10*(Y(j)-X(j))-kk*(X(j)-x(j)));
N2=N2+((n-j+1)^q2-(n-j)^q2)*((24-4*c)*X(j)-X(j)*Z(j)+c*Y(j)-kk*(Y(j)-y(j)));
N3=N3+((n-j+1)^q3-(n-j)^q3)*(X(j)*Y(j)-8*Z(j)/3-kk*(Z(j)-z(j)));
end
x1(n+1)=x0+h^q1*n1/(gamma(q1)*q1);
y1(n+1)=y0+h^q2*n2/(gamma(q2)*q2);
z1(n+1)=z0+h^q3*n3/(gamma(q3)*q3);
x(n+1)=x0+h^q1*(10*(y1(n+1)-x1(n+1))+m1)/gamma(q1+2);
y(n+1)=y0+h^q2*((24-4*c)*x1(n+1)-x1(n+1)*z1(n+1)+c*y1(n+1)+m2)/gamma(q2+2);
z(n+1)=z0+h^q3*(x1(n+1)*y1(n+1)-8*z1(n+1)/3+m3)/gamma(q3+2);
X1(n+1)=X0+h^q1*N1/(gamma(q1)*q1);
Y1(n+1)=Y0+h^q2*N2/(gamma(q2)*q2);
Z1(n+1)=Z0+h^q3*N3/(gamma(q3)*q3);
X(n+1)=X0+h^q1*(10*(Y1(n+1)-X1(n+1))-kk*(X1(n+1)-x1(n+1))+M1)/gamma(q1+2);
Y(n+1)=Y0+h^q2*((24-4*c)*X1(n+1)-X1(n+1)*Z1(n+1)+c*Y1(n+1)-kk*(Y1(n+1)-y1(n+1))+M2)/gamma(q2+2);
Z(n+1)=Z0+h^q3*(X1(n+1)*Y1(n+1)-8*Z1(n+1)/3-kk*(Z1(n+1)-z1(n+1))+M3)/gamma(q3+2);
end
figure(1)
plot(x, y)
figure(2)
plot(X,Y)
figure(3)
t1 = 0:1:N;
plot(t1,X-x,'r', t1,Y-y,'b',t1,Z-z,'g')
