!! DVEC_EVEN returns N values, evenly spaced between ALO and AHI.
! Parameters:
!
! Input, integer N, the number of values.
!
! Input, real ( kind = 8 ) ALO, AHI, the low and high values.
!
! Output, real ( kind = 8 ) A(N), N evenly spaced values.
! Normally, A(1) = ALO and A(N) = AHI.
! However, if N = 1, then A(1) = 0.5*(ALO+AHI).
subroutine dvec_even ( n, alo, ahi, a )
implicit none
integer*4 n
real*8 a(n), ahi, alo
integer*4 i
if ( n .eq. 1 ) then
a(1) = 0.5d0 * ( alo + ahi )
else
do i = 1, n
! a(i) = ( real ( n - i, kind = 8 ) * alo &
! + real ( i - 1, kind = 8 ) * ahi ) &
! / real ( n - 1, kind = 8 )
a(i) = ( dble(n - i) * alo + dble(i - 1) * ahi ) / dble(n - 1)
enddo
endif
return
end subroutine dvec_even
subroutine filon_cos ( ftab, a, b, ntab, t, result )
!xtab
!***********************************************************************
!
!! FILON_COS uses Filon's method on integrals with a cosine factor.
!
! Discussion:
!
! The integral to be approximated has the form:
!
! Integral ( A <= X <= B ) F(X) * COS(T*X) dX
!
! where T is user specified.
!
! The function is interpolated over each subinterval by
! a parabolic arc.
!
! Reference:
!
! Abramowitz and Stegun,
! Handbook of Mathematical Functions,
! pages 890-891.
!
! S M Chase and L D Fosdick,
! Algorithm 353, Filon Quadrature,
! Communications of the Association for Computing Machinery,
! Volume 12, 1969, pages 457-458.
!
! Philip Davis and Philip Rabinowitz,
! Methods of Numerical Integration,
! Blaisdell Publishing, 1967.
!
! Parameters:
!
! Input, real ( kind = 8 ) FTAB(NTAB), contains the value of the function
! at A, A+H, A+2*H, ... , B-H, B, where H = (B-A)/(NTAB-1).
!
! Input, real ( kind = 8 ) A, B, the limits of integration.
!
! Input, integer NTAB, the number of data points.
! NTAB must be odd, and greater than 1.
!
! Input, real ( kind = 8 ) T, the multiplier of the X argument of the cosine.
!
! Output, real ( kind = 8 ) RESULT, the approximate value of the integral.
!
implicit none
external dvec_even
integer*4 ntab
real*8 ftab(ntab), a, b, t, result
real*8 xtab(ntab)
integer*4 i
real*8 alpha, beta, gamma
real*8 c2n, c2nm1, cost, sint, theta, h
if ( a .eq. b ) then
result = 0.0d0
return
endif
if ( ntab .le. 1 ) then
result = 0.0d0
return
endif
! NTAB must be odd
if ( mod ( ntab, 2 ) .ne. 1 ) then
result = 0.0d0
return
end if
! set up a vector of the NTAB X values
call dvec_even ( ntab, a, b, xtab )
h = ( b - a ) / dble( ntab - 1)
theta = t * h
sint = dsin ( theta )
cost = dcos ( theta )
! small theta
if ( 6.0d0 * dabs ( theta ) .le. 1.0d0 ) then
alpha = 2.0d0 * theta**3 / 45.0d0
& - 2.0d0 * theta**5 / 315.0d0
& + 2.0d0 * theta**7 / 4725.0d0
beta = 2.0d0 / 3.0d0
& + 2.0d0 * theta**2 / 15.0d0
& - 4.0d0 * theta**4 / 105.0d0
& + 2.0d0 * theta**6 / 567.0d0
& - 4.0d0 * theta**8 / 22275.0d0
gamma = 4.0d0 / 3.0d0
& - 2.0d0 * theta**2 / 15.0d0
& + theta**4 / 210.0d0
& - theta**6 / 11340.0d0
else
alpha = ( theta**2 + theta * sint * cost - 2.0d0 * sint**2 )
& / theta**3
beta = ( 2.0d0 * theta + 2.0d0 * theta * cost**2
& - 4.0d0 * sint * cost ) / theta**3
gamma = 4.0d0 * ( sint - theta * cost ) / theta**3
endif
!!! c2n = sum ( ftab(1:ntab:2) * dcos ( t * xtab(1:ntab:2) ) )
!!! & - 0.5d0 * ( ftab(ntab) * dcos ( t * xtab(ntab) )
!!! & + ftab(1) * dcos ( t * xtab(1) ) )
!??? c2n = -0.5d0 * ( ftab(ntab) * dcos(t * xtab(ntab))
!??? & + ftab(1) * dcos(t * xtab(1)) )
c2n = 0.d0
do i = 1, ntab, 2
c2n = c2n + ftab(i) * dcos(t * xtab(i))
enddo
c2n = c2n - 0.5d0 * ( ftab(ntab) * dcos(t * xtab(ntab))
& + ftab(1) * dcos(t * xtab(1)) )
!!! c2nm1 = sum ( ftab(2:ntab-1:2) * dcos ( t * xtab(2:ntab-1:2) ) )
c2nm1 = 0.0d0
do i = 2, ntab-1, 2
c2nm1 = c2nm1 + ftab(i) * dcos (t * xtab(i))
enddo
result = h * (
& alpha * ( ftab(ntab) * dsin ( t * xtab(ntab) )
& - ftab(1) * dsin ( t * xtab(1) ) )
& + beta * c2n
& + gamma * c2nm1 )
return
end
! main
real*8 PI
parameter (PI = 3.14159265358979324d0)
integer*4 menu
integer*4 i, j, k, n
real*8 exact, result
real*8 a, b, t
real*8 h
real*8 rint_cs, rint_sn
integer*4 ntab
real*8 ftab(5000), xtab(5000)
!
a = 0.0d0
b = 2.0d0 * PI
!odd
ntab = 11!10 for old
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
if ( mod ( ntab, 2 ) .ne. 1 ) then
! add (?)
! ntab = ntab + 1
! or remove the last [f, x] point
ntab = ntab - 1
endif
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
call dvec_even ( ntab, a, b, xtab )
h = ( b - a ) / ( ntab - 1)
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'filon demo'
write ( *, '(a)') ' FILON_COS estimates the integral of.'
write ( *, '(a)' ) ' F(X) * COS ( T * X )'
write ( *, '(a)' ) ' '
write ( *, '(a)' ) ' Integrate F(X)*COS(T*X):'
write ( *, '(a)' ) ' with F(X)=1, X, X**2.'
write ( *, '(a)' ) ' '
write ( *, '(a,g24.16)' ) ' A = ', a
write ( *, '(a,g24.16)' ) ' B = ', b
write ( *, '(a,i6)' ) ' NTAB = ', ntab
write ( *, '(a,g24.16)' ) ' H = ', h
write ( *, '(a)' ) ' '
write ( *, '(a)' )
& ' T Approximate Exact'
write ( *, '(a)' ) ' '
do k = 1, 3 ! 3 different T
if (k .eq. 1) then
t = 1.0d0
else if (k .eq. 2) then
t = 2.0d0
else if (k .eq. 3) then
t = 10.0d0
endif
do i = 1, 3 ! 3 functions: 1, x, x^2
!-------------------
do j = 1, ntab
if (i .eq. 1)then
ftab(j) = 1.0d0
else
ftab(j) = xtab(j)**(i - 1)
endif
enddo!j
!-------------------
call filon_cos(ftab, a, b, ntab, t, result)
! menu = i
! call filon_old(menu, a, b, t, n, result, rint_sn)
if ( i .eq. 1 ) then
exact = ( dsin ( t * b ) - dsin ( t * a ) ) / t
else if ( i .eq. 2 ) then
exact = ( ( dcos ( t * b ) + t * b * dsin ( t * b ) ) -
& ( dcos ( t * a ) + t * a * dsin ( t * a ) ) )/t**2
else if ( i .eq. 3 ) then
exact = ( ( 2.0d0 * t * b * dcos ( t * b )
& + ( t * t * b**2 - 2.0d0 ) * dsin ( t * b ) )
& - ( 2.0d0 * t * a * dcos ( t * a )
& + ( t * t * a**2 - 2.0d0 ) * dsin ( t * a ) ) )
& / t**3
endif
write ( *, '(2x,3g24.16)' ) t, result, exact
end do!i
write ( *, '(a)' ) ' '
enddo!k
end
|
