主题:龙格库塔程序运行错误
下面的程序是龙格库塔程序,但运行时出了错误,结果全部都是初值。我觉得可能是外部函数有问题。要求解的方程式:dy/dx=1.5*y*cosx/(1+0.5cosx)-1.5*sinx/(1+0.5cosx)**3,初值是:x=2.29336289,y=0,x的终值是3.98982311.积分步长是h=0.03141593。请哪位帮忙指点一下?
subroutine rk4(y,dydx,n,x,h,yout,derivs)
parameter(nmax=10)
dimension y(n),dydx(n),yout(n),yt(nmax),dyt(nmax),&
dym(nmax)
real hh,h6,xh,x,h
integer i,n
hh=h*0.5
h6=h/6.
xh=x+hh
do i=1,n
yt(i)=y(i)+hh*dyt(i)
enddo
call derivs(xh,yt,dym)
do i=1,n
yt(i)=y(i)+h*dym(i)
dym(i)=dyt(i)+dym(i)
enddo
call derivs(x+h,yt,dyt)
do i=1,n
yout(i)=y(i)+h6*(dydx(i)+dyt(i)+2.*dym(i))
enddo
end subroutine rk4
subroutine rkdumb(vstart,nvar,x1,x2,nstep,derivs)
parameter(nmax=10)
common/path/ xx(200) ,y(10,200)
!uses rk4
dimension vstart(nvar),v(nmax),dv(nmax)
real x,x1,x2,h
integer i,nstep,k
external derivs
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==x)&
pause 'stepsize not signficant in rkdump.'
x=x+h
xx(k+1)=x
do i=1,nvar
y(i,k+1)=v(i)
enddo
enddo
end subroutine rkdumb
program d14r2
!drive for routine rkdumb
parameter(nvar=1,nstep=54)
dimension vstart(nvar)
common /path/x(200),y(10,200)
external derivs1
x1=2.29336298
vstart(1)=0.
x2=3.98982311
call rkdumb(vstart,nvar,x1,x2,nstep,derivs1)
write(*,'(/1x,t9,a,t17,a,t31,a/)')'x','integrated'
do i=1,nstep
write(*,'(1x,f10.8,2x,2f12.6)') x(i),y(1,i)
enddo
end program
subroutine derivs1(x,y,dydx)
dimension y(*),dydx(*)
dydx(1)=-3*a*sin(x)/(1+a*cos(x))**3+3*y(1)*a*cos(x)/(1+a*cos(x))
return
end subroutine derivs1
subroutine rk4(y,dydx,n,x,h,yout,derivs)
parameter(nmax=10)
dimension y(n),dydx(n),yout(n),yt(nmax),dyt(nmax),&
dym(nmax)
real hh,h6,xh,x,h
integer i,n
hh=h*0.5
h6=h/6.
xh=x+hh
do i=1,n
yt(i)=y(i)+hh*dyt(i)
enddo
call derivs(xh,yt,dym)
do i=1,n
yt(i)=y(i)+h*dym(i)
dym(i)=dyt(i)+dym(i)
enddo
call derivs(x+h,yt,dyt)
do i=1,n
yout(i)=y(i)+h6*(dydx(i)+dyt(i)+2.*dym(i))
enddo
end subroutine rk4
subroutine rkdumb(vstart,nvar,x1,x2,nstep,derivs)
parameter(nmax=10)
common/path/ xx(200) ,y(10,200)
!uses rk4
dimension vstart(nvar),v(nmax),dv(nmax)
real x,x1,x2,h
integer i,nstep,k
external derivs
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==x)&
pause 'stepsize not signficant in rkdump.'
x=x+h
xx(k+1)=x
do i=1,nvar
y(i,k+1)=v(i)
enddo
enddo
end subroutine rkdumb
program d14r2
!drive for routine rkdumb
parameter(nvar=1,nstep=54)
dimension vstart(nvar)
common /path/x(200),y(10,200)
external derivs1
x1=2.29336298
vstart(1)=0.
x2=3.98982311
call rkdumb(vstart,nvar,x1,x2,nstep,derivs1)
write(*,'(/1x,t9,a,t17,a,t31,a/)')'x','integrated'
do i=1,nstep
write(*,'(1x,f10.8,2x,2f12.6)') x(i),y(1,i)
enddo
end program
subroutine derivs1(x,y,dydx)
dimension y(*),dydx(*)
dydx(1)=-3*a*sin(x)/(1+a*cos(x))**3+3*y(1)*a*cos(x)/(1+a*cos(x))
return
end subroutine derivs1