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