四阶龙格库塔方程组代码，为什么不出结果？

c275633094

如下所示，运算结果不对，请问哪里有问题？

[Fortran] 纯文本查看 复制代码

module counts
 integer,parameter::nstep=100
 integer,parameter::nvar=3 !Funcion dimentions
 end module counts

 module d1
 use counts
 real::xx(nstep),y(nvar,nstep)
 end module d1


 subroutine derivs(x,y,dydx)
 implicit none
 integer::n !方程维数
real::dydx(*),x,y(*)
 dydx(1)=y(2)*y(3)
 dydx(2)=-y(1)*y(3)
 dydx(3)=-0.51*y(1)*y(2)
 end subroutine derivs


 subroutine rk4(y,dydx,n,x,h,yout,derivs)
 integer,parameter::nmax=3
 real::h,x,dydx(n),y(n),yout(n)
 external::derivs
 integer::i
 real::h6,hh,xh,dym(nmax),dyt(nmax),yt(nmax)
 hh=h*0.5
 h6=h/6.
 xh=x+hh

 do i=1,n !First step,k1=dydx(i)
 yt(i)=y(i)+hh*dyt(i)
 enddo

 call derivs(xh,yt,dyt) !Second step,k2=dyt
 do i=1,n
 yt(i)=y(i)+hh*dyt(i)
 enddo

 call derivs(xh,yt,dym) !Third step,k3=dym
 do i=1,n
 yt(i)=y(i)+h*dym(i)
 dym(i)=dyt(i)+dym(i) !k3=k3+k2
 enddo

 call derivs(x+h,yt,dyt)!Fourth step,k4=dyt
 do i=1,n
 yout(i)=y(i)+h6*(dydx(i)+dyt(i)+2.*dym(i))
 end do !Every step calculate dydx from x and above y(i)

 return
 end subroutine rk4



 subroutine rkdumb(vstart,nvar,x1,x2,nstep,derivs)
 integer::nstep,nvar
 integer,parameter::nmax=10,nstpmx=100
 real::x2,x1,vstart(nvar),xx(nstpmx),y(nmax,nstpmx)
 external derivs
 integer::i,k
 real::h,x,dv(nmax),v(nmax)
 do i=1,nvar !方程维数
v(i)=vstart(i)
 y(i,1)=v(i)
 enddo
 xx(1)=x1
 x=x1
 h=(x2-x1)/nstep !步长
do k=1,nstep
 call derivs(x,v,dv)
 call rk4(v,dv,nvar,x,h,v,derivs)
 if(x+h.eq.x) pause'stepsize not significant in rkdumb.'
 x=x+h
 xx(k+1)=x
 do i=1,nvar
 y(i,k+1)=v(i)
 enddo
 enddo
 return
 end subroutine rkdumb



 program solve319
 use counts
 use d1
 implicit none
 integer::i
 real::x2,x1,vstart(nvar)
 external::derivs,rk4,rkdumb
 vstart(1)=0
 vstart(2)=1
 vstart(3)=1
 x1=0
 x2=12
 print*,vstart(1)
 open(unit=11,file='20150327out.dat')
 do i=1,nstep
 write(unit=11,fmt='(3f12.6)')y(1,i),y(2,i),y(3,i)

 call rkdumb(vstart,nvar,x1,x2,nstep,derivs)
 enddo
 end 

necrohan

这种应该有比较通用的算法程序吧，直接找个现成的比较看看。

糖盒love玲珑

我刚好最近写了一个RK4变步长的Fortran程序, 不知道你有兴趣要没?

llsvvv

我有有兴趣要

chenk

糖盒love玲珑 发表于 2015-12-28 18:55
我刚好最近写了一个RK4变步长的Fortran程序, 不知道你有兴趣要没?

请问可求解2阶微分方程组吗？我现在需要用RK4解2阶微分方程组。。需要解决算法。。。求助