plotphase.m

基于Matlab对约瑟夫森结（Josephson Junction）RCSJ模型的交直流I-V特性及非线性混沌现象进行数值模拟。通过计算机数值模拟得到该模型的非线性微分方程数值解

%
% 画时序图、相图、庞加莱截面及功率谱
% 
% ****** Author:J.F.Peng([e-mail]) ******
%
clear
clc
global Bc ia w id; 
Bc=0.500; 
ia=1.32;
w=0.660; 
id=0.000; 
%
% 解微分方程
% 
[t,Y]=ode45('Josephson_Junction',[0,800],[0;0;0]);
%
% 画φ、dφ/dt、ψ（dψ/dt=w）的三维图
%
figure
plot3(Y(:,1),Y(:,2),Y(:,3))
xlabel('Phase  φ'); 
ylabel('Voltage  u=dφ/dt'); 
zlabel('Phase of Aternative Current ψ'); 
%
% 画归一化电压的时序图(Time Dependence Voltage)
%
miny=min(Y(:,2));
maxy=max(Y(:,2));
figure
subplot(2,2,4)
plot(t(700:end),Y(700:end,2))
axis([700 800 miny maxy])
xlabel('Time  t'); 
ylabel('Voltage  u=dφ/dt'); 
title('Time Dependence Voltage')
%
% 画相图(Phase Diagram)
%
T=2*pi/w;
[t,Y]=ode45('Josephson_Junction',[0:T/100:1000*T],[0;0;0]);
n=length(t);
m=round(n/60);
subplot(2,2,1)
%plot(Y(m:n,1),Y(m:n,2),'k')
%plot(sin(Y(m:n,1)),Y(m:n,2),'k')
%
% 如果相图窗口被涂黑可选择用下个
%
plot(Y(50*m:n,1),Y(50*m:n,2),'k')
%plot(sin(Y(50*m:n,1)),Y(50*m:n,2),'k')
axis tight
%axis([560 600 -2 2])
axis([-5600 -5550 -2 1])
title('Phase Diagram')
xlabel('Phase φ')
%xlabel('sin(φ)')
ylabel('Voltage u=dφ/dt')
%
% 画庞加莱截面(Poincare Section)
%
% 通过观察Poincare截面上截点的情况可以判断是否发生混沌：
% 当Poincare截面上有且只有一个不动点或少数离散点时，运动是周期的；
% 当Poincare截面上是一封闭曲线时，运动是准周期的；
% 当Poincare截面上是一些成片的具有分形结构的密集点时，运动便是混沌。
%
Z1=[];
Z2=[];
%
% 每隔一个外加交流电周期(T=2π/ω)取一个点
%
for i=m:100:length(Y)
    Z1=[Z1,Y(i,1)];
    Z2=[Z2,Y(i,2)];
end
subplot(2,2,2)
plot(Z1,Z2,'b.','markersize',6)
minx=min(Y(m:n,1));
maxx=max(Y(m:n,1));
miny=min(Y(m:n,2));
maxy=max(Y(m:n,2));
axis([minx maxx miny maxy])
title('Poincare Section with（T=2π/ω）')
xlabel('Phase φ')
ylabel('Voltage u=dφ/dt')
%
% 画功率谱(Power Spectrum)
%
subplot(2,2,3)
%
% Fs 采样频率
% CYn 自相关函数
%  
Fs=1;
t=0:1/Fs:1000;
[T,Y]=ode45('Josephson_Junction',t,[0;0;0]);
n=length(t);
m=round(n/6);
Yn=Y(m:n,2);
nfft=1024;
CYn=xcorr(Yn,'unbiased'); 
CYk=fft(CYn,nfft);
Pyy=abs(CYk);
index=0:round(nfft/2-1);
%
% 功率谱以f=Fs/2为对称轴，左右对称，故画图频率范围为0~Fs/2
%
k=index*Fs/nfft;
plot_Pyy=1000*Pyy(index+1);
plot(k,10*log10(plot_Pyy));
title('Power Spectrum')
xlabel('Frequency  f')
ylabel('power  |Yk|^2') 
axis tight