Duffing方程胞映射计算程序

clear;ux1=1.5;dx1=-1.5;ux2=0.8;dx2=-1.5;nx1=400;nx2=200;globalFomigaF=0.345;omiga=1.2;fori=1:nx1*nx2b=mod(i,nx1);a=ceil(i/nx1);ifb==0b=nx1;endCELL{i}=[(b-1)*(ux1-dx1)/nx1+dx1b*(ux1-dx1)/nx1+dx1(b-1)*(ux1-dx1)/nx1+dx1+0.5*(ux1-dx1)/nx1;(a-1)*(ux2-dx2)/nx2+dx2a*(ux2-dx2)/nx2+dx2(a-1)*(ux2-dx2)/nx2+dx2+0.5*(ux2-dx2)/nx2];endNg=1;g=zeros(1,nx1*nx2+1);g(nx1*nx2+1)=1;s(nx1*nx2+1)=0;c(nx1*nx2+1)=nx1*nx2+1;forz=1:nx1*nx2j=0;b=z;whileg(b)==0g(b)=-1;c(b)=getNextCellNum(CELL,b,ux1,dx1,ux2,dx2,nx1,nx2);b=c(b);j=j+1;end%ifg(b)==-1%Ng=Ng+1;%forI=0:j%g(z)=Ng;%z=c(z);%end%endifg(b)==-1Ng=Ng+1;i=0;zTemp=z;whileb~=zTempzTemp=c(zTemp);i=i+1;endforI=1;ig(z)=Ng;s(z)=i-I;z=c(z);endforI=i+1:jg(z)=Ng;s(z)=0;z=c(z);endendifg(b)0forI=0:j-1g(z)=g(b);s(z)=s(b)+j-I;z=c(z);endendendNgholdonplot(0,0,'x')colormatrix=hsv(Ng);fori=1:nx1*nx2plot(CELL{i}(1,3),CELL{i}(2,3),'.','color',colormatrix(g(i),:),'LineWidth',2)endfunctionm=getNextCellNum(CELL,n,ux1,dx1,ux2,dx2,nx1,nx2)X1=CELL{n}(1,3);X2=CELL{n}(2,3);globalFomigaperiod=2*pi/omiga;step=2*pi/omiga/100;dt=step;%m=240;a1=18.86*1e3;a2=0.13*1e6;a3=0.004*1e9;b0=33.052*1e3;b1=32.303;c0=1385.4;%c1=524.28;X3=0;x0=[X1X2X3]';options=odeset('RelTol',1e-4,'AbsTol',[1e-51e-51e-5]);[t,x]=ode45('duffing',[0:step:period],x0,options);X11=x(end,1);X22=x(end,2);if(X11dx1||X11=ux1)||(X22dx2||X22=ux2)m=nx1*nx2+1;elseb=floor((X11-dx1)/(ux1-dx1)*nx1)+1;a=floor((X22-dx2)/(ux2-dx2)*nx2)+1;m=(a-1)*nx1+b;endend