第二次改条件求常微分方程报错

yb6231290

当位移x的初值x0=-12,v0=3时，进行运动学求解，采用ivprk求解微分方程，此时加速度为2，但我想要当得到的位移xt>16时，停止计算，将停止时的时间t2变成新的t0，停止时的速度vt和xt变成新的速度和位移的初值v0,x0,然后再进行第二次求微分方程，此时加速度为10。我的程序已经列出，但是总是报错，求大神指导。

[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

li913

1、给出错误提示；
2、代码不全，没法调试；
3、最好所有文件打包上传。

yb6231290

li913 发表于 2018-4-9 21:25
1、给出错误提示；
2、代码不全，没法调试；
3、最好所有文件打包上传。

代码是完整的啊

li913

DIVPRK,DSET 没有。

yb6231290

li913 发表于 2018-4-10 10:15
DIVPRK,DSET 没有。

这个是imsl函数，估计你没装

li913

FAQ之 各类函数库的使用方法(如IMSL,MKL)
http://fcode.cn/guide-58-1.html


最重要的，给出你的错误提示！

aoeo2jam

简单问题整这么复杂？