马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。
您需要 登录 才可以下载或查看，没有账号？我要加入
x
看近来对自治系统的Poincare截面问题颇多，现将我所做的Lorenz系统的Poincare截面计算代码分享如下，由于计算量很大，所以我采用Fortran编程，计算结果和方法我在之前的帖子已有说明，请参考！
————该方法只适用于自治系统，对于非自治系统，请参考频闪法！
- ! ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
- ! 主程序
- program Poincare_section_Lorenz_fit_FRK
- implicit none
- external f
- dimension y(3),d(3),b(3),z1(2000000,3),z2(300000,3),z(2,4)
- dimension yy1(3000000,2),yy2(3000000,2),yy3(3000000,2)
- double precision y,d,t,h,b,z,z1,z2,yy1,yy2,yy3
- double precision xa,ya,za,tol_e
- integer i,j,k,l,m,n,p
- !**************************************************************************
- ! 积分求解过程
- y(1)=0
- y(2)=1
- y(3)=0
- t=0.0
- h=0.01
- m=3
- l=2000000
- ! p=200000
- do 110 i=1,l
- n=2
- call rkt2(t,y,m,h,n,z,f,d,b)
- do 15 j=1,m
- z1(i,j)=z(2,j)
- 15 continue
- t=t+h
- print *,i
- 110 continue
- ! do 115 i=1800000,l
- ! z2(i-1799999,1)=z1(i,1)
- ! z2(i-1799999,2)=z1(i,2)
- ! z2(i-1799999,3)=z1(i,3)
- ! 115 continue
- open(1,file='Phase_Lorenz.dat',status='new')
- do 120 i=1800000,l
- write(1,*)(z1(i,j),j=1,3)
- 120 continue
- close(1)
- !***************************************************************
- !!!!!! 求均值
- tol_e=1.0d-1
- xa=0.0
- DO 60 i=1,l
- xa=xa+z1(i,1)/l
- 60 continue
- ya=0.0
- DO 70 i=1,l
- ya=ya+z1(i,2)/l
- 70 continue
- za=0.0
- DO 80 i=1,l
- za=za+z1(i,3)/l
- 80 continue
- !
- !
- open(2,file='Section_Lorenz_xy.dat')
- !判断，如果y(3)-za<0.01，则将相应的y(1),y(2)值赋给yy1和yy2，然后输出
- do 130 k=1,l
- if(abs(z1(k,3)-za).lt.tol_e)then
- yy1(k,1)=z1(k,1)
- yy1(k,2)=z1(k,2)
- end if
- write(2,*)(yy1(k,j),j=1,2)
- 130 continue
- close(2)
- open(3,file='Section_Lorenz_yz.dat')
- !判断，如果y(2)<25，则将相应的y(1),y(2)值赋给yy1和yy2，然后输出
- do 140 k=1,l
- if(abs(z1(k,2)-ya).lt.tol_e)then
- yy2(k,1)=z1(k,1)
- yy2(k,2)=z1(k,3)
- end if
- write(3,*)(yy2(k,j),j=1,2)
- 140 continue
- close(3)
- open(4,file='Section_Lorenz_xz.dat')
- !判断，如果y(1)<25，则将相应的y(1),y(2)值赋给yy1和yy2，然后输出
- do 150 k=1,l
- if(abs(z1(k,1)-xa).lt.tol_e)then
- yy3(k,1)=z1(k,2)
- yy3(k,2)=z1(k,3)
- end if
- write(4,*)(yy3(k,j),j=1,2)
- 150 continue
- close(4)
- end
- ! ------------------------------------------------------------------------
- ! *****************定义Lorenz系统微分方程****************************
- subroutine f(t,y,m,d)
- dimension y(m),d(m)
- real*8 y,d,t
- real*8 a,b,c
- a=10
- b=8/3
- c=28
- d(1)=-a*(y(1)-y(2))
- d(2)=-y(1)*y(3)+c*y(1)-y(2)
- d(3)=y(1)*y(2)-b*y(3)
- return
- end
- !***************************************************************************************
- ! 四阶定步长Rounge－Kutta法积分子程序
- subroutine rkt2(t,y,m,h,n,z,f,d,b)
- dimension y(m),d(m),z(m,n),a(4),b(m)
- double precision y,d,z,a,b,t,h,x,tt
- !!!!!!!!!!!!!!!!!!!!!参数说明!!!!!!!!!!!!!!!!!!!!
- ! t——双精度实型变量，输入参数。积分起始点
- ! y——双精度实型一维数组，长度为m，输入参数。m个未知函数在起始点t处的初值
- ! 为yi(t)(i=1,2,...,m)，返回时其值在z(i,1)(i=1,2,...,m)中
- ! m——整型变量，输入参数。方程组中方程的个数，也是未知函数的个数
- ! h——双精度实型变量，输入参数。积分步长
- ! n——整型变量，输入参数。积分的步数（包括起始点这一步）
- ! z——双精度实型二维数组，体积为m*n，输出参数。返回m个未知函数在n个积分点上的函数值
- ! 即：z(i,j)=yij, i=1,2,...,m; j=1,2,...,n
- ! 其中，z(i,1)(i=1,2,...,m)为初值
- ! f——子程序名，输入参数。用于计算方程组中各方程右端函数值fi(t,y1,y2,...ym)
- ! (i=1,2,...,m)。在主程序中必须使用外部语句对相应的实参数进行说明。
- ! 子程序由用户自己编译，语句形式为：
- ! subroutine f(t,y,m,d)
- ! 其中，t——双精度实型变量，自变量值；
- ! y——双精度实型一维数组，长m，存放m个未知函数值
- ! d——双精度实型一维数组，长m，返回m个方程右端函数值
- ! d——双精度实型一维数组，长m。本程序的工作数组，用于存放m个方程的右端函数值
- ! b——双精度实型一维数组，长m。本程序的工作数组。
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- a(1)=h/2.0
- a(2)=a(1)
- a(3)=h
- a(4)=h
- do 5 i=1,m
- 5 z(i,1)=y(i)
- x=t
- do 100 j=2,n
- call f(t,y,m,d)
- do 10 i=1,m
- 10 b(i)=y(i)
- do 30 k=1,3
- do 20 i=1,m
- y(i)=z(i,j-1)+a(k)*d(i)
- b(i)=b(i)+a(k+1)*d(i)/3.0
- 20 continue
- tt=t+a(k)
- call f(tt,y,m,d)
- 30 continue
- do 40 i=1,m
- 40 y(i)=b(i)+h*d(i)/6.0
- do 50 i=1,m
- 50 z(i,j)=y(i)
- t=t+h
- 100 continue
- t=x
- return
- end
- !***************************************************************************************
|
发表于 2008-1-9 15:51
提升卡
变色卡
显身卡
发表于 2008-1-10 19:29
楼主