Fortran Coder

查看: 129|回复: 5

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

[复制链接]

7

帖子

4

主题

0

精华

入门

F 币
40 元
贡献
27 点
发表于 2018-4-9 17:09:49 | 显示全部楼层 |阅读模式
10F 币
当位移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



回复

使用道具 举报

254

帖子

1

主题

0

精华

宗师

F 币
1430 元
贡献
971 点
发表于 2018-4-9 21:25:02 | 显示全部楼层
1、给出错误提示;
2、代码不全,没法调试;
3、最好所有文件打包上传。
回复

使用道具 举报

7

帖子

4

主题

0

精华

入门

F 币
40 元
贡献
27 点
 楼主| 发表于 2018-4-10 08:19:43 | 显示全部楼层
li913 发表于 2018-4-9 21:25
1、给出错误提示;
2、代码不全,没法调试;
3、最好所有文件打包上传。

代码是完整的啊
回复

使用道具 举报

254

帖子

1

主题

0

精华

宗师

F 币
1430 元
贡献
971 点
发表于 2018-4-10 10:15:27 | 显示全部楼层
DIVPRK,DSET 没有。
回复

使用道具 举报

7

帖子

4

主题

0

精华

入门

F 币
40 元
贡献
27 点
 楼主| 发表于 2018-4-10 12:57:34 | 显示全部楼层
li913 发表于 2018-4-10 10:15
DIVPRK,DSET 没有。

这个是imsl函数,估计你没装
回复

使用道具 举报

254

帖子

1

主题

0

精华

宗师

F 币
1430 元
贡献
971 点
发表于 2018-4-10 19:14:09 | 显示全部楼层
本帖最后由 li913 于 2018-4-10 19:20 编辑

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


最重要的,给出你的错误提示!


回复

使用道具 举报

您需要登录后才可以回帖 登录 | 极速注册

本版积分规则

QQ|捐赠本站|Archiver|关于我们 About Us|群聊|Fcode

GMT+8, 2018-7-21 00:22

Powered by Discuz! X3.2

© 2001-2017 Comsenz Inc.

快速回复 返回顶部 返回列表