257 lines
6.8 KiB
Fortran
257 lines
6.8 KiB
Fortran
module numint
|
|
|
|
use const_and_precisions, only : wp_, zero, one
|
|
implicit none
|
|
|
|
contains
|
|
|
|
subroutine simpson (n,h,fi,s)
|
|
! subroutine for integration over f(x) with the simpson rule. fi:
|
|
! integrand f(x); h: interval; s: integral. copyright (c) tao pang 1997.
|
|
implicit none
|
|
integer, intent(in) :: n
|
|
real(wp_), intent(in) :: h
|
|
real(wp_), dimension(n), intent(in) :: fi
|
|
real(wp_), intent(out) :: s
|
|
integer :: i
|
|
real(wp_) :: s0,s1,s2
|
|
|
|
s = zero
|
|
s0 = zero
|
|
s1 = zero
|
|
s2 = zero
|
|
do i = 2, n-1, 2
|
|
s1 = s1+fi(i-1)
|
|
s0 = s0+fi(i)
|
|
s2 = s2+fi(i+1)
|
|
end do
|
|
s = h*(s1+4.0_wp_*s0+s2)/3.0_wp_
|
|
! if n is even, add the last slice separately
|
|
if (mod(n,2).eq.0) s = s+h*(5.0_wp_*fi(n)+8.0_wp_*fi(n-1)-fi(n-2))/12.0_wp_
|
|
end subroutine simpson
|
|
|
|
subroutine trapezoid(n,xi,fi,s)
|
|
! subroutine for integration with the trapezoidal rule.
|
|
! fi: integrand f(x); xi: abscissa x;
|
|
! s: integral Int_{xi(1)}^{xi(n)} f(x)dx
|
|
implicit none
|
|
integer, intent(in) :: n
|
|
real(wp_), dimension(n), intent(in) :: xi,fi
|
|
real(wp_), intent(out) :: s
|
|
integer :: i
|
|
|
|
s = zero
|
|
do i = 1, n-1
|
|
s = s+(xi(i+1)-xi(i))*(fi(i+1)-fi(i))
|
|
end do
|
|
s = 0.5_wp_*s
|
|
end subroutine trapezoid
|
|
|
|
subroutine quanc8(fun,a,b,abserr,relerr,result,errest,nofun,flag)
|
|
implicit none
|
|
real(wp_), intent(in) :: a, b, abserr, relerr
|
|
real(wp_), intent(out) :: result, errest, flag
|
|
integer, intent(out) :: nofun
|
|
!
|
|
! estimate the integral of fun(x) from a to b
|
|
! to a user provided tolerance.
|
|
! an automatic adaptive routine based on
|
|
! the 8-panel newton-cotes rule.
|
|
!
|
|
! input ..
|
|
!
|
|
! fun the name of the integrand function subprogram fun(x).
|
|
! a the lower limit of integration.
|
|
! b the upper limit of integration.(b may be less than a.)
|
|
! relerr a relative error tolerance. (should be non-negative)
|
|
! abserr an absolute error tolerance. (should be non-negative)
|
|
!
|
|
! output ..
|
|
!
|
|
! result an approximation to the integral hopefully satisfying the
|
|
! least stringent of the two error tolerances.
|
|
! errest an estimate of the magnitude of the actual error.
|
|
! nofun the number of function values used in calculation of result.
|
|
! flag a reliability indicator. if flag is zero, then result
|
|
! probably satisfies the error tolerance. if flag is
|
|
! xxx.yyy , then xxx = the number of intervals which have
|
|
! not converged and 0.yyy = the fraction of the interval
|
|
! left to do when the limit on nofun was approached.
|
|
!
|
|
real(wp_) :: w0,w1,w2,w3,w4,area,x0,f0,stone,step,cor11,temp
|
|
real(wp_) :: qprev,qnow,qdiff,qleft,esterr,tolerr
|
|
real(wp_), dimension(31) :: qright
|
|
real(wp_), dimension(16) :: f,x
|
|
real(wp_), dimension(8,30) :: fsave,xsave
|
|
integer :: levmin,levmax,levout,nomax,nofin,lev,nim,i,j
|
|
|
|
interface
|
|
function fun(x)
|
|
use const_and_precisions, only : wp_
|
|
implicit none
|
|
real(wp_), intent(in) :: x
|
|
real(wp_) :: fun
|
|
end function fun
|
|
end interface
|
|
!
|
|
! *** stage 1 *** general initialization
|
|
! set constants.
|
|
!
|
|
levmin = 1
|
|
levmax = 30
|
|
levout = 6
|
|
nomax = 5000
|
|
nofin = nomax - 8*(levmax-levout+2**(levout+1))
|
|
!
|
|
! trouble when nofun reaches nofin
|
|
!
|
|
w0 = 3956.0_wp_ / 14175.0_wp_
|
|
w1 = 23552.0_wp_ / 14175.0_wp_
|
|
w2 = -3712.0_wp_ / 14175.0_wp_
|
|
w3 = 41984.0_wp_ / 14175.0_wp_
|
|
w4 = -18160.0_wp_ / 14175.0_wp_
|
|
!
|
|
! initialize running sums to zero.
|
|
!
|
|
flag = zero
|
|
result = zero
|
|
cor11 = zero
|
|
errest = zero
|
|
area = zero
|
|
nofun = 0
|
|
if (a .eq. b) return
|
|
!
|
|
! *** stage 2 *** initialization for first interval
|
|
!
|
|
lev = 0
|
|
nim = 1
|
|
x0 = a
|
|
x(16) = b
|
|
qprev = zero
|
|
f0 = fun(x0)
|
|
stone = (b - a) / 16.0_wp_
|
|
x(8) = (x0 + x(16)) / 2.0_wp_
|
|
x(4) = (x0 + x(8)) / 2.0_wp_
|
|
x(12) = (x(8) + x(16)) / 2.0_wp_
|
|
x(2) = (x0 + x(4)) / 2.0_wp_
|
|
x(6) = (x(4) + x(8)) / 2.0_wp_
|
|
x(10) = (x(8) + x(12)) / 2.0_wp_
|
|
x(14) = (x(12) + x(16)) / 2.0_wp_
|
|
do j = 2, 16, 2
|
|
f(j) = fun(x(j))
|
|
end do
|
|
nofun = 9
|
|
!
|
|
! *** stage 3 *** central calculation
|
|
! requires qprev,x0,x2,x4,...,x16,f0,f2,f4,...,f16.
|
|
! calculates x1,x3,...x15, f1,f3,...f15,qleft,qright,qnow,qdiff,area.
|
|
!
|
|
do
|
|
do
|
|
x(1) = (x0 + x(2)) / 2.0_wp_
|
|
f(1) = fun(x(1))
|
|
do j = 3, 15, 2
|
|
x(j) = (x(j-1) + x(j+1)) / 2.0_wp_
|
|
f(j) = fun(x(j))
|
|
end do
|
|
nofun = nofun + 8
|
|
step = (x(16) - x0) / 16.0_wp_
|
|
qleft = (w0*(f0 + f(8)) + w1*(f(1)+f(7)) + w2*(f(2)+f(6)) &
|
|
+ w3*(f(3)+f(5)) + w4*f(4)) * step
|
|
qright(lev+1)=(w0*(f(8)+f(16))+w1*(f(9)+f(15))+w2*(f(10)+f(14)) &
|
|
+ w3*(f(11)+f(13)) + w4*f(12)) * step
|
|
qnow = qleft + qright(lev+1)
|
|
qdiff = qnow - qprev
|
|
area = area + qdiff
|
|
!
|
|
! *** stage 4 *** interval convergence test
|
|
!
|
|
esterr = abs(qdiff) / 1023.0_wp_
|
|
tolerr = max(abserr,relerr*abs(area)) * (step/stone)
|
|
if (lev .ge. levmin) then
|
|
!
|
|
! *** stage 6 *** trouble section
|
|
! number of function values is about to exceed limit.
|
|
!
|
|
if (lev .ge. levmax) then
|
|
!
|
|
! current level is levmax.
|
|
!
|
|
flag = flag + one
|
|
exit
|
|
end if
|
|
if (nofun .gt. nofin) then
|
|
nofin = 2*nofin
|
|
levmax = levout
|
|
flag = flag + (b - x0) / (b - a)
|
|
exit
|
|
end if
|
|
if (esterr .le. tolerr) exit
|
|
end if
|
|
!
|
|
! *** stage 5 *** no convergence
|
|
! locate next interval.
|
|
!
|
|
nim = 2*nim
|
|
lev = lev+1
|
|
!
|
|
! store right hand elements for future use.
|
|
!
|
|
do i = 1, 8
|
|
fsave(i,lev) = f(i+8)
|
|
xsave(i,lev) = x(i+8)
|
|
end do
|
|
!
|
|
! assemble left hand elements for immediate use.
|
|
!
|
|
qprev = qleft
|
|
do i = 1, 8
|
|
j = -i
|
|
f(2*j+18) = f(j+9)
|
|
x(2*j+18) = x(j+9)
|
|
end do
|
|
end do
|
|
!
|
|
! *** stage 7 *** interval converged
|
|
! add contributions into running sums.
|
|
!
|
|
result = result + qnow
|
|
errest = errest + esterr
|
|
cor11 = cor11 + qdiff / 1023.0_wp_
|
|
!
|
|
! locate next interval.
|
|
!
|
|
do
|
|
if (nim .eq. 2*(nim/2)) exit
|
|
nim = nim/2
|
|
lev = lev-1
|
|
end do
|
|
nim = nim + 1
|
|
if (lev .le. 0) exit
|
|
!
|
|
! assemble elements required for the next interval.
|
|
!
|
|
qprev = qright(lev)
|
|
x0 = x(16)
|
|
f0 = f(16)
|
|
do i = 1, 8
|
|
f(2*i) = fsave(i,lev)
|
|
x(2*i) = xsave(i,lev)
|
|
end do
|
|
end do
|
|
!
|
|
! *** stage 8 *** finalize and return
|
|
!
|
|
result = result + cor11
|
|
!
|
|
! make sure errest not less than roundoff level.
|
|
!
|
|
if (errest .eq. zero) return
|
|
do
|
|
temp = abs(result) + errest
|
|
if (temp .ne. abs(result)) return
|
|
errest = 2.0_wp_*errest
|
|
end do
|
|
end subroutine quanc8
|
|
|
|
end module numint |