yb6231290 发表于 2018-4-9 17:09:49

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

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

1345format(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 发表于 2018-4-9 21:25:02

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

yb6231290 发表于 2018-4-10 08:19:43

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

代码是完整的啊

li913 发表于 2018-4-10 10:15:27

DIVPRK,DSET 没有。

yb6231290 发表于 2018-4-10 12:57:34

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

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

li913 发表于 2018-4-10 19:14:09

本帖最后由 li913 于 2018-4-10 19:20 编辑

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


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


aoeo2jam 发表于 2019-5-26 15:54:31

简单问题整这么复杂？
页: [1]
查看完整版本: 第二次改条件求常微分方程报错