[Fortran] 纯文本查看 复制代码
program main
INTEGER*4 MTIME, MXPARM,NP, IDO, MC, istep
PARAMETER (MXPARM=50,MTIME=5000000000, NP=50)
REAL*8 Y(NP),YPRIME(NP),Y2(NP),YPRIME2(NP)
REAL*8 PARAM(MXPARM)
real*8 x(-300:300),xe(0:2000,-300:300)
real*8 vex(-300:300),ve(0:2000)
real*8 t(0:2000),a(0:3000)
real*8 t_stop, t_start,del_t
REAL*8 TOL,T1,T2,t11,t21
external DIVPRK,DSET,FCN1,fcn2
t_start=0d0
t_stop=5d0
del_t=0.1d0
TOL=1.0D-3
CCCCC初值:y(1)=x0=-12,y(2)=v0=3
y(1)=-12
y(2)=3
open(unit=3725,file=
1'C:\USERS\YANGBO\DESKTOP\test_tang\test_xiabishi
1\test123.dat',
1status='unknown')
CALL DSET(MXPARM, 0.0D0, PARAM, 1)
PARAM(10)=1.0D0
PARAM(4)=MTIME+1
istep=0
IDO=1
T1=T_START
mt=0
istep=istep+1
T2=T1+DEL_T
CALL DIVPRK(IDO,NP,FCN1,T1,T2,TOL,PARAM,Y)
if (y(1)>16) then
y2(1)=y(1) !新的初值,求出来的y(1)变成新的x初值
y2(2)=y(2)
CALL DSET(MXPARM, 0.0D0, PARAM, 1)
PARAM(10)=1.0D0
PARAM(4)=MTIME+1
istep=0
IDO=1
T11=t2 !y(1)>16时,t的初值变为停止时的t2
mt=0
end if
110 CONTINUE
istep=istep+1
T21=T11+DEL_T
CALL DIVPRK(IDO,NP,FCN2,T11,T21,TOL,PARAM,Y2)
write(3725,1345)t2,y(2),y(1),y2(2),y2(1)
1345 format(1x,4(x,d13.6))
IF(T2<=T_STOP) GOTO 110
print *,'end'
stop
end
subroutine FCN1 (NP, T2, Y, YPRIME)
INTEGER*4 NP,npv21,npz,npx
REAL*8 Y(NP),YPRIME(NP)
a=2
YPRIME(1) = Y(2)
YPRIME(2) = a
return
end subroutine
subroutine FCN2 (NP, T21, Y2, YPRIME2)
INTEGER*4 NP,npv21,npz,npx
REAL*8 Y2(NP),YPRIME2(NP)
a=10
YPRIME2(1) = Y2(2)
YPRIME2(2) = a
return
end subroutine