hcj 发表于 2018-11-12 18:00:04

求大神帮助，使用fortran的版本很旧，是自学的

本帖最后由 hcj 于 2018-11-12 18:04 编辑

      PROGRAM MAIN
      EXTERNAL D
      DOUBLE PRECISION D,e,D1
      DO e=0.1,1.0,0.01
      D1=D(e)
      WRITE(*,*)D1
      END DO

      END

      FUNCTION D(e)
      EXTERNAL F5
      DOUBLE PRECISION D,F5,G5,e
      CALL FLAGS2(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 F1,F2
      DOUBLE PRECISION N,V1,F1,F2,V
      N=F1(V)-F2(V1)
      RETURN

      END

      FUNCTION F2(V1)
      EXTERNAL F
      DOUBLE PRECISION F2,V1,F,A,B,G1
      A=0.0
      B=V1
      EPS=0.005
      CALL FLAGS1(A,B,F,EPS,G1)
C      WRITE(*,*)G1
      F2=G1

      END

      FUNCTION F1(V)
      EXTERNAL F
      DOUBLE PRECISION F1,F,G,V
      CALL FLAGS(F,G)
C      WRITE(*,*)G
      F1=G

      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

liudy02 发表于 2018-11-12 19:51:19

新学FORTRAN，我觉得不要学老的那一套，把彭国伦的书和查普曼的书好好学学才是正经事

li913 发表于 2018-11-13 16:07:56

你的问题是啥?

hcj 发表于 2018-11-15 09:53:34

li913 发表于 2018-11-13 16:07
你的问题是啥?

解决这个积分

li913 发表于 2018-11-16 20:48:33

没法直接做无穷积分。你得想其他法儿，比如积分变换、傅氏变换等。

liudy02 发表于 2018-11-16 22:29:04

li913 发表于 2018-11-16 20:48
没法直接做无穷积分。你得想其他法儿，比如积分变换、傅氏变换等。

应该可以用有限区间积分+收敛性判定吧，不过这事没啥意思
应该有现成的软件或者函数库干这个事情吧
页: [1]
查看完整版本: 求大神帮助，使用fortran的版本很旧，是自学的