Lorenz系统的Poincare截面计算方法
看近来对自治系统的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
!*************************************************************************************** 频闪法,是什么方法？:lol 频闪法就是根据系统激励频率进行取点的方法！ octopussheng可以推荐一些关于频闪法的资料么，谢谢！:handshake :lol
回复 楼主 的帖子
oct能否将这个程序转换成matlab的 呵呵，转成matlab计算太慢了，而且现在确实也没时间了，呵呵！因为一样可以用的
回复 4楼 的帖子
频闪法你可以看看闻邦椿的《非线性振动理论中的解析方法及工程应用》这本书里面对频闪法讲的还是比较详细的，而且附录还有C语言做的一个例子！不过我C学的很差，没有去调试，我已经在非线性板块贴出这些程序了，你可以找一下！
回复 6楼 的帖子
对于我这种只会matlab的人来说，就比较难了哈回复 8楼 的帖子
呵呵，无水可以移植一下哈！基本思想你也是知道的嘛，哈！
回复 9楼 的帖子
呵呵有这个想法，但是现在都没有时间去做这些东西了 呵呵，其实fortran还是很方便的，而且计算速度比matlab快，我推荐使用！
回复 11楼 的帖子
fortran是比较快，但是我现在都没有时间学习这个哈 哈，matlab不是也可以直接调用fortran的嘛！可以省很多功夫哦！回复 13楼 的帖子
不会折腾这些东西了还是按照古老的办法解决算了 呵呵，无水所说的古老办法是什么哦？