rk.f90

这是一个可以用来计算洛仑兹方程解的程序

subroutine rk(t0,x0,y0,z0,h,n,x)
!real k_1,k_2,k_3,k_4,D
!integer M
h2=h/2
h6=h/6
if(n.EQ.1)then 
k_1=theta*(y0-x0)
!call F(t0,x0,y0,z0,n,D_1)
!k_1=D_1
k_2=theta*(y0-x0-h2*k_1)
!call F(t0+h2,x0+h2*k_1,y0+h2*k_1,z0+h2*k_1,n,D_2)
!k_2=D_2
k_3=theta*(y0-x0-h2*k_2)
!call F(t0+h2,x0+h2*k_2,y0+h2*k_2,z0+h2*k_2,n,D_3)
!k_3=D_3
!t1=t0+h
k_4=theta*(y0-x0-h*k_3)
!call F(t1,x0+h*k_3,y0+h*k_3,z0+h*k_3,n,D_4)
!k_4=D_4
x=x0+h6*(k_1+2*(k_2+k_3)+k_4)
end if
!---------------------------------------------------
if(n.EQ.2)then 
k_1=x0*(r-z0)-y0
k_2=x0*(r-z0)-y0-h2*k_1
k_3=x0*(r-z0)-y0-h2*k_2
!t1=t0+h
k_4=x0*(r-z0)-y0-h*k_3
x=y0+h6*(k_1+2*(k_2+k_3)+k_4)
end if
!----------------------------------------------------
if(n.EQ.3)then
k_1=x0*y0-b*z0
k_2=x0*y0-b*(z0+h2*k_1)
k_3=x0*y0-b*(z0+h2*k_2)
!t1=t0+h
k_4=x0*y0-b*(z0+h*k_3)
x=z0+h6*(k_1+2*(k_2+k_3)+k_4)
end if
!k_1=D_1(t0,x0,y0,z0)
!k_2=D_1(t0+h2,x0+h2*k_1,y0+h2*k_1,z0+h2*k_1)
!k_3=D_1(t0+h2,x0+h2*k_2,y0+h2*k_2,z0+h2*k_2)
!t1=t0+h
!k_4=D_1(t1,y0+h*k_3)
!y=y0+h6*(k_1+2*(k_2+k_3)+k_4)
end