poincare_section请教与讨论

410610288

我最近在做混沌束晕方面的研究，可是混沌方面的基础比较薄弱，Matlab用的也不是很熟悉，现在遇到以下问题、请指教
function Z=f1(t,Z)  
k=0;  %k=3时为混沌态
x=Z(1);
y=Z(2);
Z=[y,k/x+1/(x^3)-kz(t)*x]; %kx(t)为周期聚焦场
%%%
function kzs=kz(s)   %上面用到的子函数程序
yitaa=1/6;  %1/6为占空比
k0=3.79;  %k0=3.79
s1=rem(abs(s),1);
if s1>yitaa/2&&s1<1-yitaa/2
    kzs=0;
else
    kzs=k0;
end
%%%因为直接调用ode45出现错误，主要可能是不知道kz(t)怎么调用进去（高手若知道，请指教），所以自己编了一下龙格库塔法的程序
function [T,Z]=rks4(F,a,b,Za,M)
    h=(b-a)/M;
    T=zeros(1,M+1);
    Z=zeros(M+1,length(Za));
    T=a:h:b;
    Z(1,:)=Za;
for j=1:M
    k1=h*feval(F,T(j),Z(j,:));
    k2=h*feval(F,T(j)+h/2,Z(j,:)+k1/2);
    k3=h*feval(F,T(j)+h/2,Z(j,:)+k2/2);
    k4=h*feval(F,T(j)+h,Z(j,:)+k3);
    Z(j+1,:)=Z(j,:)+(k1+2*k2+2*k3+k4)/6;
end
%%%下面是主程序
clear
clc
for i=1.1:15.1  %取15个初值点
[T,Z]=rks4('f1',0,500,[i;0],50000);
plot(Z(40001:100:end,1),Z(40001:100:end,2),'.');
hold on;
end
title('shuyun系统的 Poincare 映像');
xlabel('rb'),ylabel('drb/ds');
%%%
出的图与看的文章（Nonlinear Resonances and Chaotic Behavior in a Periodically Focused Intense Charged_Particle Beam，有兴趣可查阅一下）上的不一样，图片暂无法上传，请指教，谢谢！

410610288

期待高手予以援手，谢谢！

410610288

贴一下图，第一个是原文的图

第二个是我做出来的图

k=3时的混沌状态的图更不一致了，这里就不贴出来了


另外，上述程序中的kz(s)函数也可用k0/2*(1+square(2*pi*(t+1/12),100/6))更为简便

meiyongyuandeze

感觉做庞加莱截面对积分的精度要求很高，最好还是用变步长或者是高阶的龙格库塔方法，或者这是时间用ode45吧！

kangarooli

你的图这不已经挺像的了吗，至于有差别的地方可能和步长与精度有关吧，在不知道作者怎么取值的情况下想做出一模一样的图还是有困难的

410610288

回复 4 # meiyongyuandeze 的帖子

感觉做庞加莱截面对积分的精度要求很高，最好还是用变步长或者是高阶的龙格库塔方法，或者这是时间用ode45吧！

变步长是指时间上的积分区间呢，还是初值点的选取？
积分区间的话具体怎么实现呢？请指教

410610288

回复 5 # kangarooli 的帖子

姑且可以说挺像的，可我的图最外圈有“重影”。更何况那是k=0时的情况，k=3时的图根本就出现不了原图的那种效果，是不是这种情况下的poincare截面图不是这样做的，我的程序有问题吧？

meiyongyuandeze

回复 6 # 410610288 的帖子

是按照积分的精度来自动选择下个积分时间的步长。

410610288

回复 8 # meiyongyuandeze 的帖子

具体怎么实现呢？可以帮忙改一下程序吗，谢谢

410610288

自己终于改出来了，谢谢大家的帮助！

kangarooli

如果再能与大家分享下经验就更好了，呵呵