回 帖 发 新 帖 刷新版面

主题:[讨论]求龙贝格数值积分法程序

在下不才,不会编写程序,求编一个用龙贝格数值积分程序,计算任意函数,任意区段的积分,请高手拔刀相助啊
在下感激不尽!

回复列表 (共4个回复)

沙发


    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一一双精度实型变盘,输出参数。返回积分值。

板凳


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 楼

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 楼

楼上正解,楼主再试试

我来回复

您尚未登录,请登录后再回复。点此登录或注册