李亚普诺夫指数计算方法

求取lyapunov指数的小数据量方法，采用混合编程lylorenz.f90的程序如下：

program lylorenz

parameter(n=3,m=12,st=100)

integer::i,j,k

real y(m),z(n),s(n),yc(m),h,y1(m),a,b,r,f(m),k1,k2,k3 y(1)=10.

y(2)=1.

y(3)=0.

a=10.

b=8./3.

r=28.

t=0.

h=0.01

!!!!!initial conditions

do i=n+1,m

y(i)=0.

end do

do i=1,n

y((n+1)*i)=1.

s(i)=0

end do

open(1,file='lorenz1.dat')

open(2,file='ly lorenz.dat')

do 100 k=1,st !!!!!!!!st iterations

call rgkt(m,h,t,y,f,yc,y1)

!!!!normarize vector 123

z=0.

do i=1,n

do j=1,n

z(i)=z(i)+y(n*j+i)**2

enddo

if(z(i)>0.)z(i)=sqrt(z(i))

do j=1,n

y(n*j+i)=y(n*j+i)/z(i)

enddo

end do

!!!!generate gsr coefficient

k1=0.

k2=0.

k3=0.

do i=1,n

k1=k1+y(3*i+1)*y(3*i+2)

k2=k2+y(3*i+3)*y(3*i+2)

k3=k3+y(3*i+1)*y(3*i+3)

end do

!!!!conduct new vector2 and 3

do i=1,n

y(3*i+2)=y(3*i+2)-k1*y(3*i+1)

y(3*i+3)=y(3*i+3)-k2*y(3*i+2)-k3*y(3*i+1) end do

!!!generate new vector2 and 3,normarize them do i=2,n

z(i)=0.

do j=2,n

z(i)=z(i)+y(n*j+i)**2

enddo

if(z(i)>0.)z(i)=sqrt(z(i))

do j=2,n

y(n*j+i)=y(n*j+i)/z(i)

end do

end do

!!!!!!!update lyapunov exponent

do i=1,n

if(z(i)>0)s(i)=s(i)+log(z(i))

enddo

100 continue

do i=1,n

s(i)=s(i)/(1.*st*h*1000.)

write(2,*)s(i)

enddo

end

!!!!!!!!!!!!!!!!!!!!!!!!!

subroutine rgkt(m,h,t,y,f,yc,y1)

real y(m),f(m),y1(m),yc(m),a,b,r

integer::i,j

do j=1,1000

call df(m,t,y,f)

t=t+h/2.0

do i=1,m

yc(i)=y(i)+h*f(i)/2.0

y1(i)=y(i)+h*f(i)/6.0

end do

call df(m,t,yc,f)

do i=1,m

yc(i)=y(i)+h*f(i)/2.0

y1(i)=y1(i)+h*f(i)/3.0

end do

call df(m,t,yc,f)

t=t+h/2.0

do i=1,m

yc(i)=y(i)+h*f(i)

y1(i)=y1(i)+h*f(i)/3.0

end do

call df(m,t,yc,f)

do i=1,m

y(i)=y1(i)+h*f(i)/6.0

end do

if(j>500)write(1,*)t,y(1),y(2),y(3)

end do

return

end

!!!!!!!!!!!!!!!!!!!!!!!!

subroutine df(m,t,y,f)

real y(m),a,b,r,f(m)

common a,b,r

a=10.

b=8./3.

r=28.

f(1)=a*(y(2)-y(1))

f(2)=y(1)*(r-y(3))-y(2)

f(3)=y(1)*y(2)-b*y(3)

do i=0,2

f(4+i)=a*y(7+i)-y(4+i)

f(7+i)=y(4+i)*(r-y(3))-y(7+i)-y(1)*y(10+i)

f(10+i)=y(2)*y(4+i)-b*y(10+i)+y(1)*y(7+i) enddo

return

end