Poincare截面的截取讨论

410610288

我也正在找poincare截面图的资料，有哪位高手能指点一下，质子束在直线加速器螺线管产生的周期性聚焦磁场通道中运动，当质子束满足K-V分布，其包络方程。求poincare截面。谢谢！！！

张某某zyl

这个对于lornez这样的系统也能用吗？

学渣硬要撑学霸

张某某zyl 发表于 2015-1-12 16:45
这个对于lornez这样的系统也能用吗？

洛伦兹系统可以有庞加莱截面。。。你取Y=0 就可以得出截面图了

qimiaohao

无水1324 发表于 2007-7-2 15:49
是的
理论清楚了，程序应该不是很大的问题

想问你两个问题！
1，如果我的这个信号周期不知道怎么画庞加莱截面呢？FFT分析出来的频谱并非单一分量！
2，当时主频率知道的时候，因为是离散的问题，没办法实现精确的取一个周期！能帮我看看我的程序吗？

%plot the graph of the fig4.2(Book from simple to complex) of the mutual
%syn
%which include the realization Poincare map and spectre
%do the FFT first. FFT determine the main fre and then we can calculate the
%period of the signal and strobe the phase trajectory by this period.
%complete the transformation of parameter
%qm2015.4.9
clear all;
clc;
close all;
lamda1=0.1;
lamda2=0.1;
BR=0;
BD=0.04;
w1=1;
w2=0.98;%p=w2/w1
f0=w1/(2*pi);%采样频率
T=2*pi/w1;%采样周期
nn=5000;%一个周期的点数
NN=200;%计算多少格周期
NS=round(0.5*nn*NN);%从哪个点开始算庞加莱图
dt=T/nn;
[t,y]=ode23('zjzdfun_mut',[0:T/nn:NN*T],[0.01,0.01,0.01,0.01],[],lamda1,lamda2,BR,BD,w1,w2);
Ts=T/nn;
fs=1/Ts;
x1=y(NS:1:nn*NN,3)';
N=length(x1);
x2=y(NS:1:nn*NN,4)';
X1=fft(x1,N);%进行fft变换
X2=fft(x2,N);%进行fft变换
mag1=abs(X1)/(N/2);%求幅值
mag2=abs(X2)/(N/2);%求幅值
ww=(0:length(X1)-1)*fs*2*pi/length(X1);%进行对应的频率转换
figure(2)
plot(ww,mag1,ww,mag2);%做频谱图
[pk1,pki1]=max(mag1);
wx1=ww(pki1)
[pk2,pki2]=max(mag2);
wx2=ww(pki2)
axis([0.9,1.1,0,1]);%设定最后作图的坐标
xlabel('频率(Hz)');
ylabel('幅值S');
grid;
T_1=2*pi/wx1;%1信号的周期
N_T=round(T_1/dt);
figure(1)
subplot(2,1,1)
plot(t,y(:,3),t,y(:,4));%y3是x1的值，y4是x2的值
title('realiazation');
xlabel('t');ylabel('x1x2');
subplot(2,2,3)
plot(y(1:end,3),y(1:end,4));
axis([-1 1 -1 1])
xlabel('x1');ylabel('x2');%横坐标是xdot
title('相图');
subplot(2,2,4)
axis([-1 1 -1 1])
hold on
for i=NS:N_T:NN*nn;
   plot(y(i,3),y(i,4),'r.');
   hold on;
end
xlabel('x1');ylabel('x2');
title('庞加莱');



最开始的T是采样周期，fs是采样频率。最后的wx1是ode23计算出来信号的信号频率，T_1是信号的周期，我把周期除以dt再取整转化成相应的点！但是这样总有一个连续时间到离散时间的误差，而且这个误差是会累积的，这样做出来的庞加莱图就不精确了。本来应该是一个点的，却变成了一个弧形！

qimiaohao

无水1324 发表于 2007-7-2 15:49
是的
理论清楚了，程序应该不是很大的问题

我把那个zjzdfun_mut的函数也放上来吧！方便您查看，vdp方程！

%y1--y(1) y2--y(2) x1--y(3) x2--y(4)
function ydot=vdbfun(t,y,flag,lamda1,lamda2,BR,BD,w1,w2);
ydot=[(lamda1-y(3).^2)*y(1)-w1.^2*y(3)+BR*(y(4)-y(3))+BD*(y(2)-y(1))
      (lamda2-y(4).^2)*y(2)-w2.^2*y(4)+BR*(y(3)-y(4))+BD*(y(1)-y(2))
       y(1)
       y(2)      
      ];