renyiweifen.m

Duffing振子的Matlab程序

clc
clear all
global mu 
global r
global w
global w1
global fs
global noise
global signal
global yy
global last
global x1
global tspan
%%%%%%%%%%%%%%%%%%%%%  初始化 %%%%%%%%%%%%%%%%%%%%%%%%
last=500;                             %仿真时间／单位（Ｓ） 
a=0;                                  %噪声均值
b=0;                                  %噪声方差
w1=1;                                 %信号频率
w=1;                                  %参考信号频率
k=50;                                 %k倍采样
fs=k*w/(2*pi);                        %采样频率
Ts=1/fs;                              %采样间隔
T=2*pi/w1;                            %信号周期
r=0.8200;                              %内部策动力
signal=0.05;                           %信号振幅
mu=0.5;                               %阻尼系数
dd=0.001;
tspan=[0:Ts:last];
n=length(tspan);
%%%%%%%%%%%%%% 产生均值为a，方差为b的高斯随机序列 %%%%%%%%%%%%                     
n0=randn(1,n);
noise=a+sqrt(b)*n0; 
%%%%%%%%%%%%%%%%%%%%%% CW信号生成 %%%%%%%%%%%%%%%%%%%%%
arriveTim=200;               % CW信号到达时间
nof=20;                           % CW信号持续的周期数
endTim=arriveTim+nof*T;           %CW信号结束时间
arrive=arriveTim/Ts;
endc=endTim/Ts;

for j=1:arrive-1;
    yy(j)=0;
end

for j=arrive:endc;
    yy(j)=signal*cos(w1*j*Ts);
end

for j=endc+1:n;
    yy(j)=0;
end
%%%%%%%%%%%%%%%%%%%% 输入信号合成 %%%%%%%%%%%%%%%%%%%%%%%
for i=1:n;
    input(i)=yy(i)+noise(i);
end
%%%%%%%%%%%%%%%%%%%% Duffing方程求解 %%%%%%%%%%%%%%%%%%%%%%%
x0=[0;0];                           %% 注释 %%
[t,y]=ode45('ad1',tspan,x0);        %% 'ad1'输入为正弦信号 %%
x1=y(:,1);                          %%  'ad2'输入为CW信号  %%
x2=y(:,2);
%%%%%%%%%%%%%%%%%%%% 输出矩阵求解 %%%%%%%%%%%%%%%%%%%%%%
z0=[1 0 0 0 1 0];
[t,z]=ode45('ad5',tspan,z0);
z1=z(:,1);
z2=z(:,2);
z3=z(:,3); 
z4=z(:,4);
z5=z(:,5);
z6=z(:,6);
%%%%%%%%%%%%%%%%%%%%% 最大Lyapunov特征指数求解 %%%%%%%%%%%%%%%%%
for i=1:n-1;
    Y=[z1(i) z2(i);z4(i) z5(i)];
    Y=Y'*Y;
    [a,R]=eig(Y);
    lce1(i)=log(abs(R(2,2)))/(i*Ts);
end
%%%%%%%%%%%%%%%%%%%%%% 画图 %%%%%%%%%%%%%%%%%%%%%%%%%%%
subplot(3,1,1)
plot(input)
xlabel('Time')
ylabel('Input')
grid on
subplot(3,1,2)
plot(x1)
xlabel('Time')
ylabel('Output')
grid on
subplot(3,1,3)
plot(lce1)
xlabel('Time')
ylabel('The largest LCE')
grid on

save dataofz z