program Benettin
!四阶runge-kutta方法求解微分方程,并利用Benettin方法计算其最大李雅谱诺夫指数
!使用fortran65调试成功
!2013年7月18
implicit none
integer,parameter::p=500000,n=2
integer,parameter::pi=3.141592653
integer:: i,j
real(kind=8)::y(0:p,n),n_y(0:p,n)
real(kind=8)::h,d0,dis(p),lamda,le
real(kind=8)::x(n),q(n),n_x(n),n_q(n)!临时参量
real::Gauss,sigma,t
h=0.001
t=0
!!!!!!!!!!!!!!!!!!!!!!!!!!初值
y(0,1)=-0.15
y(0,2)=0.1
!y(0,3)=1
n_y(0,1)=-0.15001 !初始邻域点
n_y(0,2)=0.10001
!n_y(0,3)=1.0001
!!!!!!!!!!!!!!!!!!!!!!!!!!!!初始距离
do j=1,n
d0=d0+(n_y(0,j)-y(0,j))**2
enddo
d0=sqrt(d0)
! write(*,*) d0
open(unit=20,file="lya.txt")
do sigma=0,0.01,0.0001
open(unit=10,file="zt.txt")
do i=0,P-1
t=i*h
dis(i+1)=0.0
do j=1,n
x(j)=y(i,j)
n_x(j)=n_y(i,j)
enddo
read(10,"(f10.6)")Gauss
call runge_Kutta(h,n,x,q,Gauss,sigma,t)
call runge_Kutta(h,n,n_x,n_q,Gauss,sigma,t)
do j=1,n
y(i+1,j)=q(j)
n_y(i+1,j)=n_q(j)
enddo
do j=1,n
dis(i+1)=dis(i+1)+(n_y(i+1,j)-y(i+1,j))**2
enddo
dis(i+1)=dsqrt(dis(i+1))
lamda=d0/dis(i+1)
!取下一个邻域内的点
do j=1,n
n_y(i+1,j)=lamda*n_y(i+1,j)+(1.0-lamda)*y(i+1,j)
enddo
enddo
do i=1,p
if(i>p/5) then
! write(*,*) dis(i)/d0
le=le+dlog(dis(i)/d0)
endif
enddo
le=le/(P*h)! le=le/(P*h*4/5)
write(*,"(1x,a20,f18.8)") "最大Lyapuunov指数是:",le
write(20,"(f10.6)")le
close(10)
end do
close(20)
stop
end
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!四阶runge-kutta方法求解微分方程 !
subroutine runge_Kutta(h,M,y1,y,Gauss,sigma,t) !
implicit none !
integer::M,j
real::Gauss,sigma,t !
real(kind=8)::y1(M),y(M),c(4,M),z(M),k(M),l(M),h !
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !
do j=1,M !
k(j)=y1(j) !
enddo !
call fun(k,M,l,Gauss,sigma,t) !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !
do j=1,M !
c(1,j)=l(j) !
k(j)=c(1,j)*h/2+y1(j) !
enddo !
call fun(k,M,l,Gauss,sigma,t) !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !
do j=1,M !
c(2,j)=l(j) !
k(j)=c(2,j)*h/2+y1(j) !
enddo !
call fun(k,M,l,Gauss,sigma,t) !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !
do j=1,M !
c(3,j)=l(j) !
k(j)=c(3,j)*h+y1(j) !
enddo !
call fun(k,M,l,Gauss,sigma,t) !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !
do j=1,M !
c(4,j)=l(j) !
enddo !
!
do j=1,M !
y(j)=y1(j)+h*(c(1,j)+2*c(2,j)+2*c(3,j)+c(4,j))/6 !
enddo !
!
return !
end !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!系统方程
subroutine fun(x,m,y,Gauss,sigma,t)
implicit none
integer,parameter::pi=3.141592653
integer::m
real(kind=8)::x(m),y(m)
real::k,alpha,beta,g,A,omega
real::Gauss,sigma,t
k=-1.875
alpha=-0.0403
beta=2.6896
g=0.0147
A=0.4
omega=0.8
y(1)=x(2)
y(2)=-x(1)-k*x(1)**3-alpha*x(2)-beta*x(2)**3-g+A*omega**2*sin(omega*t)-sigma*Gauss
return
end
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!