lorenz系统的poincare截面计算方法.txt

Lorenz系统的Poincare截面计算方法。用于分析转子点映射特性

Lorenz系统的Poincare截面计算方法
看近来对自治系统的Poincare截面问题颇多，现将我所做的Lorenz系统的Poincare截面计算代码分享如下，由于计算量很大，所以我采用Fortran编程，计算结果和方法我在之前的帖子已有说明，请参考！

————该方法只适用于自治系统，对于非自治系统，请参考频闪法！
[code]! ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
!    主程序         
  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 
!***************************************************************************************[/code]