rk4.f90

来自「这是一个可以用来计算洛仑兹方程解的程序」· F90 代码 · 共 55 行

subroutine rk4(t0,x0,y0,z0,h,n,x1)
real*8 k_1,k_2,k_3,k_4,D_1,D_2,D_3,D_4,xx,yy,zz
integer n
h2=h/2
h6=h/6
xx=x0
yy=y0
zz=z0
call f(t0,xx,yy,zz,n,D_1)
     k_1=D_1
     t0=t0+h2
	 if(n.EQ.1)then
	    x0=xx+h2*k_1
        call f(t0,x0,yy,zz,1,D_2)
		     k_2=D_2
             x0=xx+h2*k_2
			 call f(t0,x0,yy,zz,1,D_3)
			      k_3=D_3
                  x0=xx+h*k_3
                  call f(t0,x0,yy,zz,1,D_4)
				       k_4=D_4
					   x1=xx+h6*(k_1+2*(k_2+k_3)+k_4)
!---------------------------------------------------
	  else if(n.EQ.2)then 
	    x0=xx
	    y0=y0+h2*k_1
	    z0=zz
        call f(t0,x0,y0,z0,2,D_2)
		     k_2=D_2
             y0=yy+h2*k_2
			 call f(t0,x0,y0,z0,2,D_3)
			      k_3=D_3
                  y0=yy+h*k_3
                  call f(t0,x0,y0,z0,2,D_4)
				       k_4=D_4
					   x1=yy+h6*(k_1+2*(k_2+k_3)+k_4)
!----------------------------------------------------
	  else if(n.EQ.3)then
	    x0=xx
	    z0=z0+h2*k_1
	    y0=yy
        call f(t0,x0,y0,z0,3,D_2)
		     k_2=D_2
             z0=zz+h2*k_2
			 call f(t0,x0,y0,z0,3,D_3)
			      k_3=D_3
                  z0=zz+h*k_3
                  call f(t0,x0,y0,z0,3,D_4)
				       k_4=D_4
					   x1=zz+h6*(k_1+2*(k_2+k_3)+k_4)
	 end if
!-----------------------------------------------------

end