李亚普诺夫指数计算方法

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