Fortran Coder

查看: 15404|回复: 5

[微积分] 积分发散怎么处理?

[复制链接]

4

帖子

3

主题

0

精华

入门

F 币
55 元
贡献
34 点
发表于 2020-9-27 16:05:28 | 显示全部楼层 |阅读模式
20F 币
本帖最后由 xuweijie 于 2020-11-4 09:09 编辑

我的积分Fortran算不了,MAthematic算的太慢,希望大神指导,可能需要一个积分插件用Mathematic算,Fortran如下:      PROGRAM MAIN
        EXTERNAL D
        DOUBLE PRECISION D,e,D1
        DO e=0.1,1.0,0.01
        D1=D(e)
      WRITE(*,*)e,D1
        END DO

        END

        FUNCTION D(e)
        EXTERNAL F5
        DOUBLE PRECISION D,F5,G5,e
        CALL FLAGS3(F5,G5,e)
C        WRITE(*,*)G5
        D=G5

        END

        FUNCTION F5(b,e)
        EXTERNAL IS,IC
        DOUBLE PRECISION F5,IS,IC,b,e
        F5=(1/3.1415926)*(COS(b*(e-IC(b))))*EXP(-1*b*IS(b))
        RETURN

        END

        FUNCTION IS(b)
        EXTERNAL F4
        DOUBLE PRECISION IS,F4,G4,b
        CALL FLAGS2(F4,G4,b)
C        WRITE(*,*)G4
        IS=G4

        END

        FUNCTION F4(V1,b)
        EXTERNAL N
        DOUBLE PRECISION F4,N,V1,b
        F4=(SIN(b*V1))*N(V1)
        RETURN

        END

        FUNCTION IC(b)
        EXTERNAL F3
        DOUBLE PRECISION IC,F3,G3,b
        CALL FLAGS2(F3,G3,b)
C        WRITE(*,*)G3
        IC=G3
       
        END

        FUNCTION F3(V1,b)
        EXTERNAL N
        DOUBLE PRECISION F3,b,V1,N
        F3=(COS(b*V1))*N(V1)
        RETURN

        END


        FUNCTION N(V1)
        EXTERNAL F
        DOUBLE PRECISION N,V1,F,A,B,G1,G
        A=0.0
        B=V1
        EPS=0.005
        CALL FLAGS1(A,B,F,EPS,G1)
        CALL FLAGS(F,G)
C        WRITE(*,*)G1
        N=G1-G
        RETURN

        END

       
      FUNCTION F(V)
        DOUBLE PRECISION F,V
        F=(3/(2*3.1415926*V))*LOG((COSH(SQRT(1/(2*V)))**2-
     *        COS(SQRT(1/(2*V)))**2)/(2*V))
        RETURN

        END

      SUBROUTINE FLAGS(F,G)
        DIMENSION T(5),C(5)
        DOUBLE PRECISION F,G,T,C,X
        DATA C/0.6790941054,1.638487956,2.769426772,
     *       4.31594400,7.104896230/
        DATA T/0.26355990,1.41340290,3.59642600,
     *       7.08580990,12.64080000/

       
        G=0.0D0
        DO 10 I=1,5
          X=T(I)
          G=G+F(X)*C(I)
10    CONTINUE

        END

        SUBROUTINE FLAGS1(A,B,F,EPS,G)            
        DIMENSION T(5),C(5)
        DOUBLE PRECISION A,B,F,G,T,C,S,P,H,AA,BB,W,X,Q
        DATA T/-0.061798459,-0.5384693101,0.0,
     *       0.5384693101,0.9061798459/
        DATA C/0.2369268851,0.4786286705,0.5688888889,
     *       0.4786286705,0.2369268851/
        M=1
        S=(B-A)*0.001
        P=0.0
10    H=(B-A)/M
        G=0.0
        DO 30 I=1,M
         AA=A+(I-1)*H
         BB=A+I*H
         W=0.0
         DO 20 J=1,5
          X=((BB-AA)*T(J)+(BB+AA))/2.0
          W=W+F(X)*C(J)
20         CONTINUE
       G=G+W
30    CONTINUE
      G=G*H/2.0
        Q=ABS(G-P)/(1.0+ABS(G))
        IF ((Q.GE.EPS).AND.(ABS(H).GT.ABS(S)))THEN
          P=G
          M=M+1
          GOTO 10
        END IF
        RETURN
       
        END

      SUBROUTINE FLAGS2(F,G,b)
        DIMENSION T(5),C(5)
        DOUBLE PRECISION F,G,T,C,X,b
        DATA C/0.6790941054,1.638487956,2.769426772,
     *       4.31594400,7.104896230/
        DATA T/0.26355990,1.41340290,3.59642600,
     *       7.08580990,12.64080000/

       
        G=0.0D0
        DO 10 I=1,5
          X=T(I)
          G=G+F(X,b)*C(I)
10    CONTINUE

        END

       
      SUBROUTINE FLAGS3(F,G,e)
        DIMENSION T(5),C(5)
        DOUBLE PRECISION F,G,T,C,X,e
        DATA C/0.6790941054,1.638487956,2.769426772,
     *       4.31594400,7.104896230/
        DATA T/0.26355990,1.41340290,3.59642600,
     *       7.08580990,12.64080000/

       
        G=0.0D0
        DO 10 I=1,5
          X=T(I)
          G=G+F(X,e)*C(I)
10    CONTINUE

        END




830

帖子

2

主题

0

精华

大宗师

F 币
3905 元
贡献
2323 点
发表于 2020-9-28 10:06:29 | 显示全部楼层
只要是数值积分,fortran都可以,难易程度而已。给出积分方程。
回复

使用道具 举报

250

帖子

2

主题

0

精华

宗师

F 币
1731 元
贡献
872 点

规矩勋章

发表于 2020-9-30 07:33:55 | 显示全部楼层
怎么判断的积分发散?
回复

使用道具 举报

4

帖子

3

主题

0

精华

入门

F 币
55 元
贡献
34 点
 楼主| 发表于 2020-9-30 08:51:28 | 显示全部楼层
necrohan 发表于 2020-9-30 07:33
怎么判断的积分发散?

算出的数非常大
回复

使用道具 举报

2

帖子

0

主题

0

精华

入门

F 币
142 元
贡献
23 点

规矩勋章

发表于 2020-10-14 00:23:19 | 显示全部楼层
首先你要判断究竟是因为什么原因发散。
如果是由于方程本身不满足李雅普诺夫稳定性导致的,那么这种发散是必然而无解的,发散速度取决于李雅普诺夫系数。
对于刚性过程,如果使用非刚性积分器来计算必然发散,那么需要使用满足适用于刚性条件的积分器来进行计算,
如果是由于截断引起的,那么应该采用更高阶的算法,或者使用更小的步长来计算。
如果是因为辛结构破坏引入的,那么应该采取保辛的算法。
回复

使用道具 举报

213

帖子

2

主题

0

精华

宗师

F 币
2126 元
贡献
875 点

规矩勋章

发表于 2020-10-15 00:43:11 | 显示全部楼层
函数 Ic(b) 在v' 是零时,无穷大。
回复

使用道具 举报

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

本版积分规则

捐赠本站|Archiver|关于我们 About Us|小黑屋|Fcode ( 京ICP备18005632-2号 )

GMT+8, 2024-10-10 15:52

Powered by Tencent X3.4

© 2013-2024 Tencent

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