program lylorenz.m

计算lyapunov指数的程序源码,利用matlab编程

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