lorenz.f90

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

parameter n=3
real y(n),f(n),y1(n),yc(n),a,b,c,t,h
integer::j,i
a=10.
b=28.
c=8./3.
t=0.
h=0.01
y(1)=10.00000
y(2)=1.
y(3)=2.
open(1,file='lorenz1.dat')
do i=1,5000
call rgkt(n,h,t,y,f,yc,y1)
end do 
do i=1,10000
call rgkt(n,h,t,y,f,yc,y1)
write(1,*)y(1),y(2),y(3)
enddo 
end 
!~!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

subroutine df(n,t,y,f)
	real y(n),f(n),a,b,c
	a=10.
b=28.
c=8./3.
	f(1)=a*(y(2)-y(1))
	f(2)=y(1)*(b-y(3))-y(2)
	f(3)=y(1)*y(2)-c*y(3)
	return
end 

!!!!!!!!!!!!!!!!!!!!!
subroutine rgkt(n,h,t,y,f,yc,y1)
real y(n),f(n),y1(n),yc(n),a,b,c
integer::i
	call df(n,t,y,f)
	t=t+h/2.0
	do i=1,n
		yc(i)=y(i)+h*f(i)/2.0
		y1(i)=y(i)+h*f(i)/6.0		
		end do			
	call df(n,t,yc,f)
	do i=1,n
		yc(i)=y(i)+h*f(i)/2.0
		y1(i)=y1(i)+h*f(i)/3.0		
	end do	
	call df(n,t,yc,f)
	t=t+h/2.0
	do i=1,n
		yc(i)=y(i)+h*f(i)
		y1(i)=y1(i)+h*f(i)/3.0
	end do
	call df(n,t,yc,f)
	do i=1,n
		y(i)=y1(i)+h*f(i)/6.0
	end do

returnend