急,高手帮我看看这个fortran求积分的程序哪错了,为什么运行不出来结果.
急,高手帮我看看这个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