DE積分 振動型積分 cos版

\int_0^\infty f(x)\cos\omg x dx
f(x)=-exp(x)

module exp_lib
implicit none
real(8) pi

contains

subroutine calf(f,x)
real(8),intent(in) :: x
real(8),intent(out) :: f
f=exp(-x)
end subroutine calf

subroutine de2_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+h*0.5d0
!!$ e6=1.d0
!!$ phi=1.d0/6.d0
!!$ dphi=1.d0/2.d0
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*cos(omg*x)*dphi
do i=1,n
t=i*h+h*0.5d0
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*cos(omg*x)*dphi
t=-t+h*0.5d0
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*cos(omg*x)*dphi
end do
sum=sum*em*h
end subroutine de2_int

end module exp_lib

program exp
use exp_lib
implicit none
external calf
integer n
real(8) em,omg,sum
pi=acos(-1.d0)
sum=0.d0
em=100.d0
omg=1.d0
n=75
call de2_int(sum,em,omg,n,calf)
write(*,*) sum
end program exp

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

DE積分

両端で特異点がある場合の積分
ちょっとした工夫で精度があがる。
jasosx.ils.uec.ac.jp/JSPF/JSPF_TEXT/jspf1990/jspf1990_05/jspf1990_05-397.pdf -

ただ、被積分関数がごちゃごちゃしている場合は、精度を落としてもいいから工夫無しで攻めるのも手かと思われる。

module lib_de
implicit none
contains
function f(x,opx,omx)
real(8) f,x,opx,omx
f=1.d0/sqrt(opx*omx) ! 1/\sqrt(1-x^2)
end function f
end module lib_de

program de
use lib_de
implicit none
real(8) h,t,x,dx,f1,f2,s,a,pih,e(3),err,opx,omx
integer i,n
s=0.d0
e(:)=1.d0
pih=asin(1.d0)
h=1.d0/8.d0
t=0.d0
x=0.d0
opx=1.d0
omx=1.d0
dx=pih
f1=f(x,opx,omx)*dx
s=s+f1*h
do i=1,100
t=i*h
a=pih*sinh(t)
x=tanh(a)
opx=exp(a)/cosh(a)
omx=exp(-a)/cosh(a)
dx=pih*cosh(t)/cosh(a)**2
f1=f(x,opx,omx)*dx
t=-i*h
a=pih*sinh(t)
x=tanh(a)
opx=exp(a)/cosh(a)
omx=exp(-a)/cosh(a)
dx=pih*cosh(t)/cosh(a)**2
f2=f(x,opx,omx)*dx
s=s+(f1+f2)*h
e(1)=e(2)
e(2)=e(3)
e(3)=max(abs(f1),abs(f2))*h
err=sum(e)/abs(s)
if(err < 1.d-17) then
n=2*(i-3)+1
exit
end if
end do
write(*,*) h,n,s
end program de