Question

急，高手帮我看看这个fortran求积分的程序哪错了，为什么运行不出来结果。。。。。

program f915
implicit none
external f1,f2,f3
double precision a,b,simps
write(*,*)'a,b'
read(*,*)a,b

write(*,*)simps(a,b,1d-7,f1)
write(*,*)simps(a,b,1d-7,f2)
write(*,*)simps(a,b,1d-7,f3)
end program

function simps(a,b,eps,f)
implicit none
double precision f
double precision a,b,simps,eps
double precision h,s1,s2,d,rc,rp
double precision x
integer m,i
h=(b-a)/2.0
s1=(f(a)+f(b)+4*f(a+h))*h/3.0
rc=f(a+h)
rp=f(a)+f(b)
m=2
do
rp=rp+2*rc
x=a-h/2.0
rc=0.0
do i=1,m
x=x+h
rc=rc+f(x)
enddo
s2=(rp+4.0*rc)*h/6.0
d=s2-s1
if(abs(s2)>=1.0)d=d/s2
if(abs(d)==0) s1=s2
h=h/2.0
m=m+m
enddo
simps=s2
end function simps

function f1(x)
implicit none
double precision x,f1
f1=1.0/(1.0+x)
end function f1

function f2(x)
implicit none
double precision x,y,f2
y=1.0/3.0
f2=0.5*((1.0+x*x)**y-(1.0-x*x)**y)
end function f2

function f3(x)
implicit none
double precision x,f3
f3=4.0/(1.0+x*x)
end function f3

Answer 1

这是想写一维变步长Simpson积分吧
1、你给代码里有一个do后面什么都没有写啊，这肯定运行不出来的。
2、if(abs(s2)>=1.0)d=d/s2
if(abs(d)==0) s1=s2 这两句不太明白想表达什么。
我把simps(a,b,eps,f)修改了下，应该可以运行出来了，LZ看下
××××××××××××××××××××××××××××××××××××××××××××××××××××××××××××××
function simps(a,b,eps,f)
implicit none
double precision f
double precision a,b,simps,eps
double precision h,s1,s2,d,rc,rp
double precision x
integer m,i
h=(b-a)/2.0
s1=(f(a)+f(b)+4*f(a+h))*h/3.0
rc=f(a+h)
rp=f(a)+f(b)
m=2
d=eps+1.0
do while(d.ge.eps)
rp=rp+2*rc
x=a-h/2.0
rc=0.0
do i=1,m
x=x+h
rc=rc+f(x)
enddo
s2=(rp+4.0*rc)*h/6.0
d=abs(s2-s1)
s1=s2
h=h/2.0
m=m+m
enddo
simps=s2
end function simps

Answer 2

do
rp=rp+2*rc
x=a-h/2.0
rc=0.0
do i=1,m
x=x+h
rc=rc+f(x)
enddo
s2=(rp+4.0*rc)*h/6.0
d=s2-s1
if(abs(s2)>=1.0)d=d/s2
if(abs(d)==0) s1=s2
h=h/2.0
m=m+m
enddo
这一段是死循环，楼主自己先改了吧。