lorenz.f90

fortran程序

!Lorenz System : 
!     dX/dt=-aX+aY
!     dY/dt=-XZ+rX-Y
!     dZ/dt=XY-bZ
!Where a, r, b is parameter coulb be changed in the follow section.
!We call on Fourth Runge-Kutta integration method. Which is the 
! class Runge-Kutta method as follow:
!     Y(i,n+1)=Y(i,n)+(h/6)*(K1+2*K2+2*K3+K4)
!     K1(i)=f(t(n),Y(1,n),Y(2,n),...,Y(m,n))
!     K2(i)=f(t(n)+h/2,Y(1,n)+h/2*K1(1),Y(2,n)+h/2*K1(2),...,Y(m,n)+h/2*K1(m))
!     K3(i)=f(t(n)+h/2,Y(1,n)+h/2*K2(1),Y(2,n)+h/2*K2(2),...,Y(m,n)+h/2*K2(m))
!     K4(i)=f(t(n)+h,Y(1,n)+h*K3(1),Y(2,n)+h*K3(2),...,Y(m,n)+h*K3(m))
!In this subroutine, we define h=dt as the step,i=1,2,...,
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!Writed by Dr. Wang.
!2006.03.09
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!Subroutine Lorenz()
implicit none
integer,parameter::Nstep=10000
integer::n,i,j

double precision::p,r,b
double precision::h
double precision::KX1,KX2,KX3,KX4
double precision::KY1,KY2,KY3,KY4
double precision::KZ1,KZ2,KZ3,KZ4
double precision::X(0:Nstep),Y(0:Nstep),Z(0:Nstep)

X(0)=1.0
Y(0)=-0.3
Z(0)=1.2

p=10.0
r=28
b=8.0/3.0


h=0.01

open(unit=1,file='outputmy.txt')

do n=0,Nstep-1
   KX1=-p*X(n)+p*Y(n)
   KY1=-X(n)*Z(n)+r*X(n)-Y(n)
   KZ1=X(n)*Y(n)-b*Z(n)

   KX2=-p*(X(n)+h/2.0*KX1-Y(n)-h/2.0*KY1)
   KY2=-(X(n)+h/2.0*KX1)*(Z(n)+h/2.0*KZ1)+r*(X(n)+h/2.0*KX1)-(Y(n)+h/2.0*KY1)
   KZ2=(X(n)+h/2.0*KX1)*(Y(n)+h/2.0*KY1)-b*(Z(n)+h/2.0*KZ1)

   KX3=-p*(X(n)+h/2.0*KX2-(Y(n)+h/2.0*KY2))
   KY3=-(X(n)+h/2.0*KX2)*(Z(n)+h/2.0*KZ2)+r*(X(n)+h/2.0*KX2)-(Y(n)+h/2.0*KY2)
   KZ3=(X(n)+h/2.0*KX2)*(Y(n)+h/2.0*KY2)-b*(Z(n)+h/2.0*KZ2)

   KX4=-p*((X(n)+h*KX3)-(Y(n)+h*KY3))
   KY4=-(X(n)+h*KX3)*(Z(n)+h*KZ3)+r*(X(n)+h*KX3)-(Y(n)+h*KY3)
   KZ4=(X(n)+h*KX3)*(Y(n)+h*KY3)-b*(Z(n)+h*KZ3)

   X(n+1)=X(n)+h/6.0*(KX1+2.0*KX2+2.0*KX3+KX4)
   Y(n+1)=Y(n)+h/6.0*(KY1+2.0*KY2+2.0*KY3+KY4)
   Z(n+1)=Z(n)+h/6.0*(KZ1+2.0*KZ2+2.0*KZ3+KZ4)

   write(1,'(i,f20.10,f20.10,f20.10)') n+1,X(n+1),Y(n+1),Z(n+1)

enddo

end