Fortran 中的 Filon 正交

Question

我正在尝试使用 Filon 方法计算调和函数。我的代码看起来像这样：

program function1
implicit none
 real(10) :: g, a, x, b, k, s
 ! frequency intervals for function and her transform
 real(10) :: step, dstep
 ! number of inervals
 integer, parameter :: nmax=10
 integer :: j
 ! function and result of her transformation
 real(kind=10), dimension(0:nmax) :: f, wynik
 step=0.01
 dstep=0.1
 !constant value for function
 a=1
 b=-1
 s=2
 open(1, file='wykresik.dat', action='write', status='replace')
 ! loop over all frequencies in the function
do j=0,1000
  k=float(j)*step
  f=-1*cos(k*s)
enddo
  call filonc(step,dstep,nmax,f,wynik)
write(1,*) k, wynik
! filon subroutine call
end program function1
include 'filon.f95'

包含的子程序是用于高振荡函数的 filons 正交函数，它看起来像这样：

SUBROUTINE FILONC ( DT, DOM, NMAX, D, CHAT )
    INTEGER :: NMAX
        REAL(kind=10) :: DT, DOM
        REAL(kind=10), dimension(0:NMAX) ::D, CHAT
    REAL(kind=10) :: TMAX, OMEGA, THETA, SINTH, COSTH, CE, CO
        REAL(kind=10) :: SINSQ, COSSQ, THSQ, THCUB, ALPHA, BETA, GAMMA
        INTEGER    TAU, NU

      IF ( MOD ( NMAX, 2 ) .NE. 0 ) then

      stop ' NMAX SHOULD BE EVEN '

      ENDIF

        TMAX = float(NMAX) * DT

       DO NU = 0, NMAX

           OMEGA = float(NU) * DOM
           THETA = OMEGA * DT



           SINTH = SIN ( THETA )
           COSTH = COS ( THETA )
           SINSQ = SINTH * SINTH
           COSSQ = COSTH * COSTH
           THSQ  = THETA * THETA
           THCUB = THSQ * THETA

           IF ( THETA .EQ. 0.0 ) THEN

              ALPHA = 0.0
              BETA  = 2.0 / 3.0
              GAMMA = 4.0 / 3.0

            ELSE

              ALPHA = ( 1.0 / THCUB )* ( THSQ + THETA * SINTH * COSTH - 2.0 * SINSQ )

              BETA  = ( 2.0 / THCUB )* ( THETA * ( 1.0 + COSSQ ) -2.0 * SINTH * COSTH )

              GAMMA = ( 4.0 / THCUB ) * ( SINTH - THETA * COSTH )

           ENDIF


       CE = 0.0

           DO  TAU = 0, NMAX, 2

              CE = CE + D(TAU) * COS ( THETA * float( TAU ) )

       ENDDO



           CE = CE - 0.5 * ( D(0) + D(NMAX) * COS ( OMEGA * TMAX ) )


           CO = 0.0

           DO  TAU = 1, NMAX - 1, 2

              CO = CO + D(TAU) * COS ( THETA * REAL ( TAU ) )

       ENDDO



           CHAT(NU) = 2.0 * ( ALPHA * D(NMAX) * SIN ( OMEGA * TMAX )+ BETA * CE + GAMMA * CO ) * DT

           ENDDO
ENDSUBROUTINE FILONC

我想知道，为什么在该代码中我会得到 11 个 1 k 的结果，我不知道为什么会这样工作。

Answer 1

这一行

write(1,*) k, wynik

写出标量 k 和整个 rank-1 数组 wynik 的值。在这两行

 integer, parameter :: nmax=10
 ...
 real(kind=10), dimension(0:nmax) :: f, wynik

wynik 被声明为有 11 个元素。

因此，自然地，程序会为 k 写出一个值，为 wynik

写出 11 个值

这看起来太明显了，我怀疑我误解了你真正要问的内容 - 在这种情况下，请尽量更清楚地了解你想知道的内容。

Fortran 中的 Filon 正交

Filon's quadrature in fortran

fortran

fft

fortran90