主题:龙格库塔程序运行问题
下面的程序是用龙格库塔解方程组的,但是运行时出现了:'stepsize not significant in rkdumb.'
请问谁能帮忙修改一下?
!RK4
subroutine rk4(y,dydx,n,x,h,yout,derivs1)
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
external derivs1
hh=h*0.5
h6=h/6.
xh=x+hh
do i=1,n
yt(i)=y(i)+hh*dyt(i)
enddo
call derivs1(xh,yt,dym)
do i=1,n
yt(i)=y(i)+h*dym(i)
dym(i)=dyt(i)+dym(i)
enddo
call derivs1(x+h,yt,dyt)
do i=1,n
yout(i)=y(i)+h6*(dydx(i)+dyt(i)+2.*dym(i))
enddo
end subroutine rk4
!rkdumb
subroutine rkdumb(vstart,nvar,x1,x2,nstep,derivs1)
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 derivs1
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 derivs1(x,v,dv)
call rk4(v,dv,nvar,x,h,v,derivs1)
if(x+h==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
end subroutine rkdumb
program d14r2
implicit double precision (a-h,o-z)
!drive for routine rkdumb
parameter(nvar=2,nstep=200)
dimension vstart(nvar)
common /path/x(200),y(10,200)
external derivs1
x1=0.0
vstart(1)=0.740741
vstart(2)=0.0
x2=6.283186
call rkdumb(vstart,nvar,x1,x2,nstep,derivs1)
write(*,'(/1x,t9,a,t17,a/)')'x','integrated'
do i=1,nstep
write(*,'(1x,f12.8,2x,2f12.8)') x(i),y(1,i),y(2,i)
enddo
end program
subroutine derivs1(x,y,dydx)
dimension y(*),dydx(*)
real x
dydx(1)=f(x,y(1))
dydx(2)=f1(x,y(1))
return
end subroutine derivs1
function f(x,y1)
real f,x,y1
f=1.5*y1*sin(x)/(1+0.5*cos(x))-1.5*sin(x)/(1+0.5*cos(x))**3
end function f
!函数过程f1
function f1(x,y1)
real f1,x,y1
f1=y1
end function f1
请问谁能帮忙修改一下?
!RK4
subroutine rk4(y,dydx,n,x,h,yout,derivs1)
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
external derivs1
hh=h*0.5
h6=h/6.
xh=x+hh
do i=1,n
yt(i)=y(i)+hh*dyt(i)
enddo
call derivs1(xh,yt,dym)
do i=1,n
yt(i)=y(i)+h*dym(i)
dym(i)=dyt(i)+dym(i)
enddo
call derivs1(x+h,yt,dyt)
do i=1,n
yout(i)=y(i)+h6*(dydx(i)+dyt(i)+2.*dym(i))
enddo
end subroutine rk4
!rkdumb
subroutine rkdumb(vstart,nvar,x1,x2,nstep,derivs1)
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 derivs1
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 derivs1(x,v,dv)
call rk4(v,dv,nvar,x,h,v,derivs1)
if(x+h==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
end subroutine rkdumb
program d14r2
implicit double precision (a-h,o-z)
!drive for routine rkdumb
parameter(nvar=2,nstep=200)
dimension vstart(nvar)
common /path/x(200),y(10,200)
external derivs1
x1=0.0
vstart(1)=0.740741
vstart(2)=0.0
x2=6.283186
call rkdumb(vstart,nvar,x1,x2,nstep,derivs1)
write(*,'(/1x,t9,a,t17,a/)')'x','integrated'
do i=1,nstep
write(*,'(1x,f12.8,2x,2f12.8)') x(i),y(1,i),y(2,i)
enddo
end program
subroutine derivs1(x,y,dydx)
dimension y(*),dydx(*)
real x
dydx(1)=f(x,y(1))
dydx(2)=f1(x,y(1))
return
end subroutine derivs1
function f(x,y1)
real f,x,y1
f=1.5*y1*sin(x)/(1+0.5*cos(x))-1.5*sin(x)/(1+0.5*cos(x))**3
end function f
!函数过程f1
function f1(x,y1)
real f1,x,y1
f1=y1
end function f1