ly lorenz.f90

附件中是我用Fortran写的lorenz混沌吸引子的lyapunov指数谱产生程序。包括三部分内容：如何产生lorenz吸引子

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 the first vector 
	z=0.   
	do j=1,n
	z(1)=z(1)+y(n*j+1)**2
	enddo
	if(z(1)>0.)z(1)=sqrt(z(1))
	do j=1,n
    y(n*j+1)=y(n*j+1)/z(1)
	enddo
  
   !!!!generate gsr coefficient
   k1=y(4)*y(5)+y(7)*y(8)+y(10)*y(11)
   k2=y(5)*y(6)+y(8)*y(9)+y(11)*y(12)
   k3=y(4)*y(6)+y(7)*y(9)+y(10)*y(12)
!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+3)=y(3*i+3)-k2*y(3*i+2)-k3*y(3*i+1)
		y(3*i+2)=y(3*i+2)-k1*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.)

enddo
print*,s(1),s(2),s(3)
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

	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