DE積分 オイラーの定数

Fortranを使って振動型の積分
\int_0^\infty f(x)\sim\omg x dx
hが刻み幅
nは打ち切る項数
omgは周期
Mはhとomgから決まる数
詳しくは数値解析、森正武、共立出版p214


module euler_lib
implicit none
real(8) pi

contains

subroutine calf(f,x)
real(8),intent(in) :: x
real(8),intent(out) :: f
if(x == 0.d0) then
f=0.d0
else
f=log(x)
end if
end subroutine calf

subroutine de_int(sum,em,omg,n,calf)
integer,intent(in) :: n
real(8),intent(in) :: em,omg
real(8),intent(out) :: sum
real(8) h,t,e6,phi,dphi,x,f
integer i
h=pi/(em*omg)
t=0.d0
e6=1.d0
phi=1.d0/6.d0
dphi=1.d0/2.d0
x=em*phi
call calf(f,x)
sum=sum+f*sin(omg*x)*dphi
do i=1,n
t=i*h
e6=exp(-6.d0*sinh(t))
phi=t/(1.d0-e6)
dphi=(1.d0+(-6.d0*t*cosh(t)-1.d0)*e6)/(1.d0-e6)**2
x=em*phi
call calf(f,x)
sum=sum+f*sin(omg*x)*dphi
t=-t
e6=exp(-6.d0*sinh(t))
phi=t/(1.d0-e6)
dphi=(1.d0+(-6.d0*t*cosh(t)-1.d0)*e6)/(1.d0-e6)**2
x=em*phi
call calf(f,x)
sum=sum+f*sin(omg*x)*dphi
end do
sum=sum*em*h
end subroutine de_int

end module euler_lib

program euler
use euler_lib
implicit none
external calf
integer n
real(8) em,omg,sum
pi=acos(-1.d0)
sum=0.d0
em=50.d0
omg=1.d0
n=75
call de_int(sum,em,omg,n,calf)
write(*,*) sum
end program euler