Lorenz系统混沌同步

Lorenz系统PC同步法

Lorenz系统

$$\{\begin{array}{l}\dot{x}=\sigma (y-x)\\ \dot{y}=x(\rho -z)-y\\ \dot{z}=xy-\beta z\end{array}$$

同步过程描述

Drive System -> Response System
采用一套动力学方程描述
D系统的x1分量完全输入R系统，作为R系统的x分量
D系统和R系统的初值不同
采用4阶Lunge-Kutta积分，精度0.001

代码

% Date : 2023/10/24
% Author: Zhou
% Ref: Pecora L M, Carroll T L. Synchronization of chaotic Systems[J]. Chaos: An Interdisciplinary Journal of Nonlinear Science, 2015, 25(9): 097611.
clc;clear;close all;
%% Time define
% time step
dt=0.001;
maxtime = 20;
maxtime_pts=round(maxtime/dt);
t_eval=linspace(0,maxtime,maxtime_pts);
%% Lorenz Generator
% Drive System
xt0 = [17+rand(1);18+rand(1);20+rand(1)];
[~,XD] = LKOde45(@ODE_Lorenz,[0 maxtime],xt0,dt);

% AWGN
XDv  = XD(1,:);
SNRdB = 100;
XDv = awgn(XDv,SNRdB,'measured');

% Response System
xt0 = [5+rand(1);4+rand(1);6+rand(1)];
[t,XR] = LKOde45Response(@ODE_Lorenz2,[0 maxtime],xt0,XDv,dt);

%% visualization
%  Error
labels = {'x','y','z'};
figure;set(gcf,'Color','w','Position',[100 100 400 300]); 
for i = 1:3
    subplot(3,1,i);plot(t,XD(i,:) - XR(i,:),'k');
    ylabel(labels(i))
end
xlabel('time');

% Time Domain
figure;set(gcf,'Color','w','Position',[100 100 400 300]); 
subplot(423);plot(t,XD(1,:));ylabel('xR');title('Drive System');
subplot(425);plot(t,XD(2,:));ylabel('yR');
subplot(427);plot(t,XD(3,:));ylabel('zR');
xlabel('time');
subplot(424);plot(t,XR(1,:));ylabel('xD');title('Response System');
subplot(426);plot(t,XR(2,:));ylabel('yD');
subplot(428);plot(t,XR(3,:));ylabel('zD');
xlabel('time');
subplot(4,2,[1,2]);plot(t,XD(2,:));ylabel('x+n');title('Drive Signal');

% D vs R
figure;set(gcf,'Color','w','Position',[100 100 400 120]); 
subplot(131);plot(XD(1,:),XR(1,:));title('xR VS xD');
subplot(132);plot(XD(2,:),XR(2,:));title('yR VS yD');
subplot(133);plot(XD(3,:),XR(3,:));title('zR VS zD');

% hyperplane
figure;set(gcf,'Color','w','Position',[100 100 400 300]); 
plot3(XD(1,:),XD(2,:),XR(2,:));
xlabel('x_D');ylabel('y_D');zlabel('y_R');
box on;
grid on;

% 4 order classic form of Runge_Kutta 
% To solve the normal differential functions y'=f(x,y),y(x0)=y0
% dyfun----vector function f(x,y)
% tspan----solve span [x0,xn]
% x0----initial value    h----step size  
% t ----return point     y---return result vector
function [t,y]=LKOde45(dyfun,tspan,x0,h)
t = tspan(1):h:tspan(2);% time
y = zeros(length(x0),length(t));% output
y(:,1) = x0(:);
for n=1:(length(t)-1)
    k1 = dyfun(t(n),y(:,n));
    k2 = dyfun(t(n)+h/2,y(:,n)+h*k1/2);
    k3 = dyfun(t(n)+h/2,y(:,n)+h*k2/2);
    k4 = dyfun(t(n)+h,y(:,n)+h*k3);
    y(:,n+1) = y(:,n)+h*(k1+2*k2+2*k3+k4)/6;
end

function df1dt = ODE_Lorenz(~,X)
df1dt = zeros(3,1);
sigma = 10;
beta = 8/3;
rho = 28;
df1dt(1) = -sigma*X(1) + sigma*X(2); 
df1dt(2) = rho*X(1) - X(2) - X(1)*X(3);
df1dt(3) = -beta*X(3) + X(1)*X(2);
end

% 4 order classic form of Runge_Kutta 
% To solve the normal differential functions y'=f(x,y),y(x0)=y0
% dyfun----vector function f(x,y)
% tspan----solve span [x0,xn]
% x0----initial value    h----step size  
% t ----return point     y---return result vector
function [t,y]=LKOde45Response(dyfun,tspan,x0,X1v,h)
t = tspan(1):h:tspan(2);% time
y = zeros(length(x0),length(t));% output
y(:,1) = x0(:);
for n=1:(length(t)-1)
    X1 = X1v(n);
    k1 = dyfun(t(n),y(:,n),X1);
    k2 = dyfun(t(n)+h/2,y(:,n)+h*k1/2,X1);
    k3 = dyfun(t(n)+h/2,y(:,n)+h*k2/2,X1);
    k4 = dyfun(t(n)+h,y(:,n)+h*k3,X1);
    y(:,n+1) = y(:,n)+h*(k1+2*k2+2*k3+k4)/6;
end

function df1dt = ODE_Lorenz2(~,X,X1)
df1dt = zeros(3,1);
sigma = 10;
beta = 8/3;
rho = 28;
df1dt(1) = -sigma*X(1) + sigma*X(2); % x2 = f_1(x2,y2,z2)
df1dt(2) = rho*X1 - X(2) - X1*X(3);  % y2 = f_2(x1,y2,z2)
df1dt(3) = -beta*X(3) + X1*X(2);     % z2 = f_3(x1,y2,z2)
end

[1]: Pecora L M, Carroll T L. Synchronization of chaotic Systems[J]. Chaos: An Interdisciplinary Journal of Nonlinear Science, 2015, 25(9): 097611.