积分无法计算

hcj

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

查看代码

打印?

001	PROGRAM MAIN

002	EXTERNAL D

003	DOUBLE PRECISION D,e,D1

004	DO e=0.1,1.0,0.01

005

D1=D(e)

006	WRITE(*,*)D1

007

END DO

008

009

END

010

011	FUNCTION D(e)

012	EXTERNAL F5

013	DOUBLE PRECISION D,F5,G5,e

014	CALL FLAGS2(F5,G5,e)

015	C        WRITE(*,*)G5

016

D=G5

017

018

END

019

020	FUNCTION F5(b,e)

021	EXTERNAL IS,IC

022	DOUBLE PRECISION F5,IS,IC,b,e

023	F5=(1/3.1415926)*(COS(b*(e-IC(b))))*EXP(-1*b*IS(b))

024

RETURN

025

026

END

027

028	FUNCTION IS(b)

029	EXTERNAL F4

030	DOUBLE PRECISION IS,F4,G4,b

031	CALL FLAGS2(F4,G4,b)

032	C        WRITE(*,*)G4

033

IS=G4

034

035

END

036

037	FUNCTION F4(V1,b)

038	EXTERNAL N

039	DOUBLE PRECISION F4,N,V1,b

040	F4=(SIN(b*V1))*N(V1)

041

RETURN

042

043

END

044

045	FUNCTION IC(b)

046	EXTERNAL F3

047	DOUBLE PRECISION IC,F3,G3,b

048	CALL FLAGS2(F3,G3,b)

049	C        WRITE(*,*)G3

050

IC=G3

051

052

END

053

054	FUNCTION F3(V1,b)

055	EXTERNAL N

056	DOUBLE PRECISION F3,b,V1

057	F3=(COS(b*V1))*N(V1)

058

RETURN

059

060

END

061

062	FUNCTION N(V1)

063	EXTERNAL F1,F2

064	DOUBLE PRECISION N,V1,F1,F2,V

065	N=F1(V)-F2(V1)

066

RETURN

067

068

END

069

070	FUNCTION F2(V1)

071	EXTERNAL F

072	DOUBLE PRECISION F2,V1,F,A,B,G1

073

A=0.0

074

B=V1

075

EPS=0.005

076	CALL FLAGS1(A,B,F,EPS,G1)

077	C        WRITE(*,*)G1

078

F2=G1

079

080

END

081

082	FUNCTION F1(V)

083	EXTERNAL F

084	DOUBLE PRECISION F1,F,G,V

085	CALL FLAGS(F,G)

086	C        WRITE(*,*)G

087

F1=G

088

089

END

090

091	FUNCTION F(V)

092	DOUBLE PRECISION F,V

093	F=(3/(2*3.1415926*V))*LOG((COSH(SQRT(1/(2*V)))**2-

094	*        COS(SQRT(1/(2*V)))**2)/(2*V))

095

RETURN

096

097

END

098

099	SUBROUTINE FLAGS(F,G)

100	DIMENSION T(5),C(5)

101	DOUBLE PRECISION F,G,T,C,X

102	DATA C/0.6790941054,1.638487956,2.769426772,

103	*       4.31594400,7.104896230/

104	DATA T/0.26355990,1.41340290,3.59642600,

105	*       7.08580990,12.64080000/

106

107

108

G=0.0D0

109	DO 10 I=1,5

110

X=T(I)

111	G=G+F(X)*C(I)

112	10    CONTINUE

113

114

END

115

116	SUBROUTINE FLAGS1(A,B,F,EPS,G)

117	DIMENSION T(5),C(5)

118	DOUBLE PRECISION A,B,F,G,T,C,S,P,H,AA,BB,W,X,Q

119	DATA T/-0.061798459,-0.5384693101,0.0,

120	*       0.5384693101,0.9061798459/

121	DATA C/0.2369268851,0.4786286705,0.5688888889,

122	*       0.4786286705,0.2369268851/

123

M=1

124	S=(B-A)*0.001

125

P=0.0

126	10    H=(B-A)/M

127

G=0.0

128	DO 30 I=1,M

129	AA=A+(I-1)*H

130

BB=A+I*H

131

W=0.0

132	DO 20 J=1,5

133	X=((BB-AA)*T(J)+(BB+AA))/2.0

134	W=W+F(X)*C(J)

135	20         CONTINUE

136

G=G+W

137	30    CONTINUE

138

G=G*H/2.0

139	Q=ABS(G-P)/(1.0+ABS(G))

140	IF ((Q.GE.EPS).AND.(ABS(H).GT.ABS(S)))THEN

141

P=G

142

M=M+1

143

GOTO 10

144

END IF

145

RETURN

146

147

END

148

149	SUBROUTINE FLAGS2(F,G,b)

150	DIMENSION T(5),C(5)

151	DOUBLE PRECISION F,G,T,C,X,b

152	DATA C/0.6790941054,1.638487956,2.769426772,

153	*       4.31594400,7.104896230/

154	DATA T/0.26355990,1.41340290,3.59642600,

155	*       7.08580990,12.64080000/

156

157

158

G=0.0D0

159	DO 10 I=1,5

160

X=T(I)

161	G=G+F(X,b)*C(I)

162	10    CONTINUE

163

164

END

带着未知量积分，把未知量一直顺延下去，在最后的一个积分里积掉，剩余一个未知量。然后做一个do循环。
这个方法为什么无法计算。与已知数据无法拟合

楚香饭

1. 我建议你换一本教科书。
都什么年代了，还用这么古老的代码风格，真是痛心。
放弃固定格式吧，放弃全部大写吧，放弃DATA语句吧，放弃 DO 数字 CONTINUE 吧
老代码尚且可以忍受，自己书写代码，为什么还要用这种古老的风格呢？？

2. 一定要书写 Implicit None，否则你出错了自己都不知道。比如 F3 函数里，你忘了定义 N 是 double，编译器误以为是 integer，结果就出现不可预料的后果。

3. 你的问题结果非常大，double 已经不足以承受。
特别是 F5 函数
  F5=(1/3.1415926)*(COS(b*(e-IC(b))))*EXP(-1*b*IS(b))
当 b 到 7.8 以后，IS 返回就是 -190 多。EXP 已经 10的500次方了。

hcj

楚香饭 发表于 2018-11-9 09:53
1. 我建议你换一本教科书。
都什么年代了，还用这么古老的代码风格，真是痛心。
放弃固定格式吧，放弃全部 ...

谢谢。。。。。。。。。。。。