主题:[讨论]求龙贝格数值积分法程序
fjbcuit
[专家分:0] 发布于 2011-06-12 16:27:00
在下不才,不会编写程序,求编一个用龙贝格数值积分程序,计算任意函数,任意区段的积分,请高手拔刀相助啊
在下感激不尽!
回复列表 (共4个回复)
沙发
bshine1225 [专家分:720] 发布于 2011-06-13 09:46:00
SUBROUTINE FROMB(A,B,F,EPS,T)
DIMENSION Y(10)
DOUBLE PRECISION A,B,F,T,Y,H,P,S,Q
H=B-A
Y(1)=H*(F(A)+F(B))/2.0
M=1
N=1
10 P=0.0
DO 20 I=0,N-1
20 P=P+F(A+(I+0.5)*H)
P=(Y(1)+H*P)/2.0
S=1.0
DO 30 K=1,M
S=4*S
Q=(S*P-Y(K))/(S-1)
Y(K)=P
P=Q
30 CONTINUE
IF ((ABS(Q-Y(M)).GE.EPS).AND.(M.LE.9)) THEN
M=M+1
Y(M)=Q
N=N+N
H=H/2.0
GOTO 10
END IF
T=Q
RETURN
END
A.B一一均为双实型变量,输入参数。积分的下限与上限。
F一~双精度实型函数:子程序名,输入参数。用于计算被积函数值f(x)。在主程序中
必须用外部语句及类型说明语句对相应的实参进行说明。该子程序由用户自编,其语句形
式为
DOUBLE PRECISI0N FUNCTION F(X)
其中:X 为双精度实型变最,自变量值;前数名F 返回双精度其型函数值。
EPS一一实型变量,输入参数。控制精度要求。
T一一双精度实型变盘,输出参数。返回积分值。
板凳
fjbcuit [专家分:0] 发布于 2011-06-13 17:43:00
subroutine rombrg(f,a,b,n,er)
dimension r(n,n)
h=b-a
r(1,1)=0.5*h*(f(b)+f(a))
open(10,file='out1.dat')
write(10,5) r(1,1)
l=1
do 4 i=2,n
h=0.5*h
l=l+l
sum=0
do 2 k=1,l-1,2
sum=sum+f(a+real(k)*h)
2 continue
r(1,1)=0.5*r(l-1,1)+h*sum
m=1
n=1
do 3 j=2,l
m=4*m
r(l,j)=r(l,j-1)+(r(l,j-1)-r(i-1,j))/real(m-l)
3 continue
write(10,5)(r(i,j),j=1,i)
if(abs(r(i,j)-r(i,j-1)).le.er) return
4 continue
close(10)
5 format(5x,5e14.7)
return
end
function f(x)
if(x==0) then
f=1
else
f=sin(x)/x
endif
return
end
external f
dimension r(5,5)
data er,a,b,n/0.00001,0,1,5/
call rombrg(f,a,b,n,r,er)
stop
end
这个是龙贝格求积分程序,好像有错误啊,我看不出来,帮帮忙看看啊,大哥!
错误如下:
C:\MSDEV\Projects\long\Text1.f90 : error FOR2233: wrong number of arguments to procedure ROMBRG invoked from main: 6 found, 5 expected
3 楼
fortran2008 [专家分:750] 发布于 2011-06-13 20:46:00
subroutine rombrg(f,a,b,n,er)
dimension r(n,n)
h=b-a
r(1,1)=0.5*h*(f(b)+f(a))
open(10,file='out1.dat')
write(10,5) r(1,1)
l=1
do 4 i=2,n
h=0.5*h
l=l+l
sum=0
do 2 k=1,l-1,2
sum=sum+f(a+real(k)*h)
2 continue
r(1,1)=0.5*r(l-1,1)+h*sum
m=1
n=1
do 3 j=2,l
m=4*m
r(l,j)=r(l,j-1)+(r(l,j-1)-r(i-1,j))/real(m-l)
3 continue
write(10,5)(r(i,j),j=1,i)
if(abs(r(i,j)-r(i,j-1)).le.er) return
4 continue
close(10)
5 format(5x,5e14.7)
return
end
function f(x)
if(x==0) then
f=1
else
f=sin(x)/x
endif
return
end
external f
!dimension r(5,5) !这可删除,是无用的变量
data er,a,b,n/0.00001,0,1,5/
call rombrg(f,a,b,n,er)
stop
end
倒数第三行的call rombrg(f,a,b,n,r,er) 改为 call rombrg(f,a,b,n,er)
4 楼
bshine1225 [专家分:720] 发布于 2011-06-14 10:06:00
楼上正解,楼主再试试
我来回复