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

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

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

[Fortran] 纯文本查看 复制代码

01        program main
02         INTEGER*4 MTIME, MXPARM,NP, IDO, MC, istep
03        PARAMETER (MXPARM=50,MTIME=5000000000, NP=50)
04        REAL*8 Y(NP),YPRIME(NP),Y2(NP),YPRIME2(NP)
05      REAL*8 PARAM(MXPARM)
06        real*8 x(-300:300),xe(0:2000,-300:300)
07        real*8 vex(-300:300),ve(0:2000)
08        real*8 t(0:2000),a(0:3000)
09        real*8 t_stop, t_start,del_t
10        REAL*8 TOL,T1,T2,t11,t21
11      external DIVPRK,DSET,FCN1,fcn2
12 
13        t_start=0d0
14        t_stop=5d0
15        del_t=0.1d0
16        TOL=1.0D-3
17CCCCC初值:y(1)=x0=-12,y(2)=v0=3
18       y(1)=-12
19        y(2)=3
20 
21      open(unit=3725,file=
22        1'C:\USERS\YANGBO\DESKTOP\test_tang\test_xiabishi
23     1\test123.dat',
24     1status='unknown')
25 
26      CALL DSET(MXPARM, 0.0D0, PARAM, 1)
27        PARAM(10)=1.0D0
28        PARAM(4)=MTIME+1
29        istep=0
30        IDO=1
31        T1=T_START
32        mt=0
33 
34        istep=istep+1
35        T2=T1+DEL_T
36        CALL DIVPRK(IDO,NP,FCN1,T1,T2,TOL,PARAM,Y)   
37 
38        if (y(1)>16) then
39         
40        y2(1)=y(1)   !新的初值,求出来的y(1)变成新的x初值
41        y2(2)=y(2)
42         
43        CALL DSET(MXPARM, 0.0D0, PARAM, 1)
44        PARAM(10)=1.0D0
45        PARAM(4)=MTIME+1
46        istep=0
47        IDO=1
48        T11=t2    !y(1)>16时，t的初值变为停止时的t2
49        mt=0
50        end if
51110        CONTINUE
52 
53      istep=istep+1
54        T21=T11+DEL_T
55         
56      CALL DIVPRK(IDO,NP,FCN2,T11,T21,TOL,PARAM,Y2)
57         
58        write(3725,1345)t2,y(2),y(1),y2(2),y2(1)
59 
601345  format(1x,4(x,d13.6))
61 
62         
63        IF(T2<=T_STOP) GOTO 110
64         
65        print *,'end'
66        stop
67        end
68 
69 
70 
71 
72      subroutine FCN1 (NP, T2, Y, YPRIME)
73      INTEGER*4 NP,npv21,npz,npx
74        REAL*8 Y(NP),YPRIME(NP)
75        a=2
76      YPRIME(1) = Y(2)
77      YPRIME(2) = a
78      return
79      end subroutine
80 
81 
82        subroutine FCN2 (NP, T21, Y2, YPRIME2)
83      INTEGER*4 NP,npv21,npz,npx
84        REAL*8 Y2(NP),YPRIME2(NP)
85        a=10
86      YPRIME2(1) = Y2(2)
87      YPRIME2(2) = a
88      return
89      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函数，估计你没装

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

简单问题整这么复杂？

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