您所在位置：主题：大家帮我看看为什么这个程序运行不下去啊 谢谢大家了
yy0451392yy
[专家分：0] 发布于 2010-04-27 17:33:00
x和y的初值都取0.7
program main
implicit real*8(a-h,o-z)
parameter(pai=3.1415926535897932384626433832795D0,ge=9.81)
! parameter(f0=0.12,beta=3,b=0.9,omig0=25,omig=2.98)
parameter(u=0.04, ksi=0.12,delta=0.16)
parameter(nmax=4,neqmax=21)
dimension x(nmax)
open(1,file='G:\zuhe3WATJB_GENGX\zuhe3WATJB_GENGX\Jeffcott.dat')
call input(kc,x)
call solve(kc,x)
close(1)
stop
end
subroutine input(kc,x)
implicit real*8(a-h,o-z)
parameter(pai=3.1415926535897932384626433832795D0,ge=9.81)
parameter(hrk0=2.*pai/100.,hmin=1.e-7,hmax=2.)
parameter(nmax=4,neqmax=21)
common/sys/ex,ey,ez,w,dw0,wstat,wend,t
common/system/wm1,wm2,wm3,xl,xl1,xl2,ei,wjox1,wjoy1,wjox2,wjoy2
common/eps/epsrk,epsi,epsw
dimension x(nmax),di(20)
read(1,*) kc
read(1,*) ex,ey,ez,rpm,drpm0,rpmstat,rpmend,t
read(1,*) epsrk,epsi,epsw
read(1,*) (x(i),i=1,4)
write(2,*) kc
write(2,*) ex,ey,ez,rpm,drpm0,rpmstat,rpmend,t
write(2,*) epsrk,epsi,epsw
write(2,*) (x(i),i=1,4)
t=t*2.*pai
return
end
!----------------------------------------------------------------------
subroutine outorbit(kpfile,x,tt)
implicit real*8(a-h,o-z)
external diffun
parameter(pai=3.1415926535897932384626433832795D0,ge=9.81,
& delta=0.16)
parameter(hrk0=2.*pai/100.,hmin=1.e-7,hmax=2.)
parameter(nmax=4,neqmax=21)
common/sys/ex,ey,ez,w,dw0,wstat,wend,t
common/system/wm1,wm2,wm3,xl,xl1,xl2,ei,wjox1,wjoy1,wjox2,wjoy2
common/eps/epsrk,epsi,epsw
dimension xc(neqmax),x(nmax)
dimension dy(neqmax),ym(3,neqmax),ysave(neqmax),y1(neqmax),
& y2(neqmax),y3(neqmax),dyn(neqmax),dy1(neqmax),yc(neqmax)
dimension ffx(10000),ffy(10000)
character kpfile*6
open(11,file=kpfile)
xc(1)=0.
do i=1,4
xc(i+1)=x(i)
enddo
hrk=hrk0
errrk=1.
do ii=1,5
ym(1,ii)=1.
enddo
k=0
! write(11,'(5f13.6)') xc(1),(xc(i)/delta,i=2,3),(xc(i)/delta/
! & omigs,i=4,5)
write(11,'(5f13.6)') xc(1),(xc(i),i=2,3),(xc(i),i=4,5)
call atomrk(5,diffun,hrk,-1,1,hmin,hmax,epsrk,errrk,xc,
& dy,ym,ysave,dyn,dy1,y1,y2,y3,yc)
40 call atomrk(5,diffun,hrk,1,1,hmin,hmax,epsrk,errrk,xc,
& dy,ym,ysave,dyn,dy1,y1,y2,y3,yc)
write(*,*)'xc(1)=',xc(1)
if(xc(1).lt.hrk0*(k+1)) goto 40
hsave=hrk
hrk=hrk0*(k+1)-ym(2,1)
call atomrk(5,diffun,hrk,0,-1,hmin,hmax,epsrk,errrk,xc,
& dy,ym,ysave,dyn,dy1,y1,y2,y3,yc)
k=k+1
hrk=hsave
! write(11,'(5f13.6)') xc(1),(xc(i)/delta,i=2,3),(xc(i)/delta/
! & omigs,i=4,5)
write(11,'(5f13.6)') xc(1),(xc(i),i=2,3),(xc(i),i=4,5)
if(xc(1).lt.tt) goto 40
close(11)
return
end
!----------------------------------------------------------------------
subroutine solve(kc,x)
implicit real*8(a-h,o-z)
parameter(nmax=4,neqmax=21)
parameter(pai=3.1415926535897932384626433832795D0,ge=9.81)
parameter(c=0.000432)
common/sys/ex,ey,ez,w,dw0,wstat,wend,t
common/system/wm1,wm2,wm3,xl,xl1,xl2,ei,wjox1,wjoy1,wjox2,wjoy2
common/eps/epsrk,epsi,epsw
dimension x(nmax),hw(nmax),xmin(nmax),xmax(nmax)
dimension ahz(nmax+2,nmax+2),ah(nmax+2,nmax+2),js(nmax+2)
dimension fai(nmax,nmax),wr(nmax),wi(nmax),fm(nmax),
& scalee(nmax),intt(nmax)
if(kc.eq.0) then
call outorbit('p1.dat',x,500.*t)
elseif(kc.eq.2) then
elseif(kc.eq.3) then
elseif(kc.eq.4) then
endif
return
end
!----------------------------------------------------------------------
subroutine diffun(l,xc,xd)
implicit real*8(a-h,o-z)
parameter(nmax=5,neqmax=5)
parameter(pai=3.1415926535897932384626433832795D0,ge=9.81)
parameter(f0=0.12,beta=3,b=0.1,omig0=25,omig=2.98,u=0.25,
& u0=0.25,ksi=0.12,delta=0.16,a=0)
dimension xc(neqmax),xd(neqmax)
gra=ge/omig0/omig0/delta
gamma=a*delta*delta
e=sqrt(xc(2)*xc(2)+xc(3)*xc(3))
if(e.le.1)then
hh=0
else
hh=1
end if
ffx=-hh*(e-1)*(xc(2)-(f0+b*delta*omig0*sqrt(xc(4)*xc(4)+
& xc(5)*xc(5)))*xc(3))*beta/e
ffy=-hh*(e-1)*((f0+sqrt(xc(4)*xc(4)+xc(5)*xc(5))*b*delta*omig0)*
& xc(2)+xc(3))* beta/e
xd(1)=1
xd(2)=xc(4)
xd(3)=xc(5)
xd(4)=ffx+u0*omig*omig*cos(omig*xc(1))-2*ksi*xc(4)-xc(2)-gamma*
& (xc(2)*xc(2)+xc(3)*xc(3))*xc(2)
xd(5)=ffy+u0*omig*omig*sin(omig*xc(1))-gra-2*ksi*xc(5)-xc(3)-
& gamma*(xc(2)*xc(2)+xc(3)*xc(3))*xc(3)
return
end
!----------------------------------------------------------------------
!------------------------atomrk solve----------------------------------
subroutine atomrk(n,diffun,h,l,jstart,hmin,hmax,eps,errbe,
& y,dy,ym,ysave,dyn,dy1,y1,y2,y3,yc)
implicit real*8(a-h,o-z)
parameter(nmax=4,neqmax=21)
external diffun
dimension y(neqmax),dy(neqmax),ym(3,neqmax),ysave(neqmax),
& y1(neqmax),y2(neqmax),y3(neqmax),dyn(neqmax),dy1(neqmax),
& yc(neqmax)
eps=abs(eps)
if(jstart.lt.0) goto 20
if(l)100,50,50
100 continue
do 10 i=1,n
10 ysave(i)=y(i)
call diffun(n,y,dy)
return
20 do 30 i=1,n
ysave(i)=ym(2,i)
ym(1,i)=ym(3,i)
30 continue
call diffun(n,ysave,dy)
50 call rk(n,diffun,h,ysave,dy,dy1,yc,y1)
if(l.eq.0) goto 90
130 hhalf=0.5*h
call rk(n,diffun,hhalf,ysave,dy,dy1,yc,y2)
call diffun(n,y2,dyn)
call rk(n,diffun,hhalf,y2,dyn,dy1,yc,y3)
errmax=0.
do 60 i=1,n
if(ym(1,i).lt.abs(y1(i))) ym(1,i)=abs(y1(i))
if(ym(1,i).lt.abs(y2(i))) ym(1,i)=abs(y2(i))
if(ym(1,i).lt.abs(y3(i))) ym(1,i)=abs(y3(i))
err1=abs(y3(i)-y1(i))/(abs(eps)*ym(1,i))
if(errmax.lt.err1) errmax=err1
y(i)=(32.*y3(i)-y1(i))/31.
60 continue
if(errmax.lt.1.0) goto 150
do 70 i=1,n
70 y1(i)=y2(i)
h=0.5*h
if(abs(h).gt.abs(hmin)) goto 130
write(*,11)
11 format(3x,'the step length is less then (hmin) !')
if(eps.lt.0.) goto 140
eps=-eps
goto 130
150 if(32.*errmax.lt.errbe) h=2.*h
if(abs(h).gt.abs(hmax)) then
h=hmax
write(*,12)
12 format(3x,'the step length is great then (hmax) !')
endif
errbe=errmax
140 do 160 i=1,n
ym(2,i)=ysave(i)
ym(3,i)=ym(1,i)
160 continue
goto 100
90 do 110 i=1,n
110 y(i)=y1(i)
goto 140
end
subroutine rk(n,diffun,h,y,dy,dy1,yc,y1)
implicit real*8(a-h,o-z)
parameter(nmax=4,neqmax=21)
external diffun
dimension y(neqmax),dy(neqmax),dy1(neqmax),yc(neqmax),
& y1(neqmax),a(4)
a(1)=h/2.
a(2)=a(1)
a(3)=h
a(4)=h
do 10 i=1,n
dy1(i)=dy(i)
y1(i)=y(i)
10 continue
do 30 k=1,3
do 20 i=1,n
yc(i)=y(i)+a(k)*dy1(i)
y1(i)=y1(i)+a(k+1)*dy1(i)/3.
20 continue
call diffun(n,yc,dy1)
30 continue
do 40 i=1,n
40 y1(i)=y1(i)+a(1)*dy1(i)/3.
return
end
program main
implicit real*8(a-h,o-z)
parameter(pai=3.1415926535897932384626433832795D0,ge=9.81)
! parameter(f0=0.12,beta=3,b=0.9,omig0=25,omig=2.98)
parameter(u=0.04, ksi=0.12,delta=0.16)
parameter(nmax=4,neqmax=21)
dimension x(nmax)
open(1,file='G:\zuhe3WATJB_GENGX\zuhe3WATJB_GENGX\Jeffcott.dat')
call input(kc,x)
call solve(kc,x)
close(1)
stop
end
subroutine input(kc,x)
implicit real*8(a-h,o-z)
parameter(pai=3.1415926535897932384626433832795D0,ge=9.81)
parameter(hrk0=2.*pai/100.,hmin=1.e-7,hmax=2.)
parameter(nmax=4,neqmax=21)
common/sys/ex,ey,ez,w,dw0,wstat,wend,t
common/system/wm1,wm2,wm3,xl,xl1,xl2,ei,wjox1,wjoy1,wjox2,wjoy2
common/eps/epsrk,epsi,epsw
dimension x(nmax),di(20)
read(1,*) kc
read(1,*) ex,ey,ez,rpm,drpm0,rpmstat,rpmend,t
read(1,*) epsrk,epsi,epsw
read(1,*) (x(i),i=1,4)
write(2,*) kc
write(2,*) ex,ey,ez,rpm,drpm0,rpmstat,rpmend,t
write(2,*) epsrk,epsi,epsw
write(2,*) (x(i),i=1,4)
t=t*2.*pai
return
end
!----------------------------------------------------------------------
subroutine outorbit(kpfile,x,tt)
implicit real*8(a-h,o-z)
external diffun
parameter(pai=3.1415926535897932384626433832795D0,ge=9.81,
& delta=0.16)
parameter(hrk0=2.*pai/100.,hmin=1.e-7,hmax=2.)
parameter(nmax=4,neqmax=21)
common/sys/ex,ey,ez,w,dw0,wstat,wend,t
common/system/wm1,wm2,wm3,xl,xl1,xl2,ei,wjox1,wjoy1,wjox2,wjoy2
common/eps/epsrk,epsi,epsw
dimension xc(neqmax),x(nmax)
dimension dy(neqmax),ym(3,neqmax),ysave(neqmax),y1(neqmax),
& y2(neqmax),y3(neqmax),dyn(neqmax),dy1(neqmax),yc(neqmax)
dimension ffx(10000),ffy(10000)
character kpfile*6
open(11,file=kpfile)
xc(1)=0.
do i=1,4
xc(i+1)=x(i)
enddo
hrk=hrk0
errrk=1.
do ii=1,5
ym(1,ii)=1.
enddo
k=0
! write(11,'(5f13.6)') xc(1),(xc(i)/delta,i=2,3),(xc(i)/delta/
! & omigs,i=4,5)
write(11,'(5f13.6)') xc(1),(xc(i),i=2,3),(xc(i),i=4,5)
call atomrk(5,diffun,hrk,-1,1,hmin,hmax,epsrk,errrk,xc,
& dy,ym,ysave,dyn,dy1,y1,y2,y3,yc)
40 call atomrk(5,diffun,hrk,1,1,hmin,hmax,epsrk,errrk,xc,
& dy,ym,ysave,dyn,dy1,y1,y2,y3,yc)
write(*,*)'xc(1)=',xc(1)
if(xc(1).lt.hrk0*(k+1)) goto 40
hsave=hrk
hrk=hrk0*(k+1)-ym(2,1)
call atomrk(5,diffun,hrk,0,-1,hmin,hmax,epsrk,errrk,xc,
& dy,ym,ysave,dyn,dy1,y1,y2,y3,yc)
k=k+1
hrk=hsave
! write(11,'(5f13.6)') xc(1),(xc(i)/delta,i=2,3),(xc(i)/delta/
! & omigs,i=4,5)
write(11,'(5f13.6)') xc(1),(xc(i),i=2,3),(xc(i),i=4,5)
if(xc(1).lt.tt) goto 40
close(11)
return
end
!----------------------------------------------------------------------
subroutine solve(kc,x)
implicit real*8(a-h,o-z)
parameter(nmax=4,neqmax=21)
parameter(pai=3.1415926535897932384626433832795D0,ge=9.81)
parameter(c=0.000432)
common/sys/ex,ey,ez,w,dw0,wstat,wend,t
common/system/wm1,wm2,wm3,xl,xl1,xl2,ei,wjox1,wjoy1,wjox2,wjoy2
common/eps/epsrk,epsi,epsw
dimension x(nmax),hw(nmax),xmin(nmax),xmax(nmax)
dimension ahz(nmax+2,nmax+2),ah(nmax+2,nmax+2),js(nmax+2)
dimension fai(nmax,nmax),wr(nmax),wi(nmax),fm(nmax),
& scalee(nmax),intt(nmax)
if(kc.eq.0) then
call outorbit('p1.dat',x,500.*t)
elseif(kc.eq.2) then
elseif(kc.eq.3) then
elseif(kc.eq.4) then
endif
return
end
!----------------------------------------------------------------------
subroutine diffun(l,xc,xd)
implicit real*8(a-h,o-z)
parameter(nmax=5,neqmax=5)
parameter(pai=3.1415926535897932384626433832795D0,ge=9.81)
parameter(f0=0.12,beta=3,b=0.1,omig0=25,omig=2.98,u=0.25,
& u0=0.25,ksi=0.12,delta=0.16,a=0)
dimension xc(neqmax),xd(neqmax)
gra=ge/omig0/omig0/delta
gamma=a*delta*delta
e=sqrt(xc(2)*xc(2)+xc(3)*xc(3))
if(e.le.1)then
hh=0
else
hh=1
end if
ffx=-hh*(e-1)*(xc(2)-(f0+b*delta*omig0*sqrt(xc(4)*xc(4)+
& xc(5)*xc(5)))*xc(3))*beta/e
ffy=-hh*(e-1)*((f0+sqrt(xc(4)*xc(4)+xc(5)*xc(5))*b*delta*omig0)*
& xc(2)+xc(3))* beta/e
xd(1)=1
xd(2)=xc(4)
xd(3)=xc(5)
xd(4)=ffx+u0*omig*omig*cos(omig*xc(1))-2*ksi*xc(4)-xc(2)-gamma*
& (xc(2)*xc(2)+xc(3)*xc(3))*xc(2)
xd(5)=ffy+u0*omig*omig*sin(omig*xc(1))-gra-2*ksi*xc(5)-xc(3)-
& gamma*(xc(2)*xc(2)+xc(3)*xc(3))*xc(3)
return
end
!----------------------------------------------------------------------
!------------------------atomrk solve----------------------------------
subroutine atomrk(n,diffun,h,l,jstart,hmin,hmax,eps,errbe,
& y,dy,ym,ysave,dyn,dy1,y1,y2,y3,yc)
implicit real*8(a-h,o-z)
parameter(nmax=4,neqmax=21)
external diffun
dimension y(neqmax),dy(neqmax),ym(3,neqmax),ysave(neqmax),
& y1(neqmax),y2(neqmax),y3(neqmax),dyn(neqmax),dy1(neqmax),
& yc(neqmax)
eps=abs(eps)
if(jstart.lt.0) goto 20
if(l)100,50,50
100 continue
do 10 i=1,n
10 ysave(i)=y(i)
call diffun(n,y,dy)
return
20 do 30 i=1,n
ysave(i)=ym(2,i)
ym(1,i)=ym(3,i)
30 continue
call diffun(n,ysave,dy)
50 call rk(n,diffun,h,ysave,dy,dy1,yc,y1)
if(l.eq.0) goto 90
130 hhalf=0.5*h
call rk(n,diffun,hhalf,ysave,dy,dy1,yc,y2)
call diffun(n,y2,dyn)
call rk(n,diffun,hhalf,y2,dyn,dy1,yc,y3)
errmax=0.
do 60 i=1,n
if(ym(1,i).lt.abs(y1(i))) ym(1,i)=abs(y1(i))
if(ym(1,i).lt.abs(y2(i))) ym(1,i)=abs(y2(i))
if(ym(1,i).lt.abs(y3(i))) ym(1,i)=abs(y3(i))
err1=abs(y3(i)-y1(i))/(abs(eps)*ym(1,i))
if(errmax.lt.err1) errmax=err1
y(i)=(32.*y3(i)-y1(i))/31.
60 continue
if(errmax.lt.1.0) goto 150
do 70 i=1,n
70 y1(i)=y2(i)
h=0.5*h
if(abs(h).gt.abs(hmin)) goto 130
write(*,11)
11 format(3x,'the step length is less then (hmin) !')
if(eps.lt.0.) goto 140
eps=-eps
goto 130
150 if(32.*errmax.lt.errbe) h=2.*h
if(abs(h).gt.abs(hmax)) then
h=hmax
write(*,12)
12 format(3x,'the step length is great then (hmax) !')
endif
errbe=errmax
140 do 160 i=1,n
ym(2,i)=ysave(i)
ym(3,i)=ym(1,i)
160 continue
goto 100
90 do 110 i=1,n
110 y(i)=y1(i)
goto 140
end
subroutine rk(n,diffun,h,y,dy,dy1,yc,y1)
implicit real*8(a-h,o-z)
parameter(nmax=4,neqmax=21)
external diffun
dimension y(neqmax),dy(neqmax),dy1(neqmax),yc(neqmax),
& y1(neqmax),a(4)
a(1)=h/2.
a(2)=a(1)
a(3)=h
a(4)=h
do 10 i=1,n
dy1(i)=dy(i)
y1(i)=y(i)
10 continue
do 30 k=1,3
do 20 i=1,n
yc(i)=y(i)+a(k)*dy1(i)
y1(i)=y1(i)+a(k+1)*dy1(i)/3.
20 continue
call diffun(n,yc,dy1)
30 continue
do 40 i=1,n
40 y1(i)=y1(i)+a(1)*dy1(i)/3.
return
end