龙格库塔法求解简单的微分方程组出错了
我应用'[何光渝]'这本书里面的四阶龙格库塔法求解一个简单的常微分方程组,但是得到的数值解与经典解完全不一样,请高手给予指点,小妹不胜感激!方程组如下: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法程序比较一下结果是否一致
如果不一致则可能是你程序输入错误
页:
[1]