龙格库塔法求解简单的微分方程组出错了

fanghuikeer

我应用'[Visual Fortran常用数值算法集][何光渝]'这本书里面的四阶龙格库塔法求解一个简单的常微分方程组,但是得到的数值解与经典解完全不一样,请高手给予指点,小妹不胜感激!
   方程组如下:dy1/dt=y1
              dy2/dt=-y2
   初值：t=0时 y1(0)=1,y2(0)=1
　我编的Fortran程序如下：
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!              Runge-Kutta算法           !!!!!!!
!!!!                     　　　　　　　　　!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      　　SUBROUTINE RKT(y,dydx,n,x,h,yout,F)
DIMENSION Y(n),dydx(n),yout(n),yt(n),dyt(n),dym(n)
real hh,h6,xh,x,h
integer i,n

hh=h/2
h6=h/6
xh=x+hh

do i=1,n
yt(i)=y(i)+hh*dydx(i)
end do
call F(xh,yt,dyt)

do i=1,n
      yt(i)=y(i)+hh*dyt(i)
end do
call F(xh,yt,dym)
do i=1,n
yt(i)=y(i)+h*dym(i)
dym(i)=dyt(i)+dym(i)
end do
call F(x+h,yt,dyt)
do i=1,n
yout=y(i)+h6*(dydx(i)+dyt(i)+2*dym(i))
end do
end SUBROUTINE RKT
!!!!!!!!!!!!!!!!!!!!!!!
!!!!      main   !!!!            
!!!!               !!!!
!!!!!!!!!!!!!!!!!!!!!!!
program d14r1
      external F
      parameter(n=2)
dimension y(n),dydx(n),yout(n)
x=0.0
y(1)=1
y(2)=1

dydx(1)=y(1)
dydx(2)=-y(2)

do i=1,11
h=0.01*(i-1)
call RKT(y,dydx,n,x,h,yout,F)
write(*,*)'t=',h
write(*,*)(yout(j),j=1,2)
end do
end program
subroutine F(x,y,dydx)
dimension y(2),dydx(2)
dydx(1)=y(1)
dydx(2)=-y(2)
end subroutine F

[ 本帖最后由 xinyuxf 于 2007-9-20 20:28 编辑 ]

风花雪月

先调用fortran在带库中的rk法程序比较一下结果是否一致
如果不一致则可能是你程序输入错误