fortran - 使用 FFTW3 库评估 FORTRAN 中高斯函数的快速傅里叶变换

Question

我正在尝试编写 FORTRAN 代码来评估高斯函数的快速傅里叶变换 f(r)=exp(-(r^2))使用 FFTW3图书馆。众所周知，高斯函数的傅里叶变换是另一个高斯函数。

我考虑评估球坐标中高斯函数的傅里叶变换积分。

因此，得到的积分可以简化为 [r*exp(-(r^2))*sin(kr)]dr 的积分。 .

我编写了以下 FORTRAN 代码来评估离散正弦变换 DST，它是使用纯实输入数组的离散傅里叶变换 DFT。夏令时由 C_FFTW_RODFT00 执行存在于 FFTW3 , 考虑到位置空间中的离散值是 r=i*delta (i=1,2,...,1024)，而 DST 的输入数组是函数 r*exp(-(r^2))不是高斯。 [r*exp(-(r^2))*sin(kr)]dr 积分中的正弦函数由球坐标上的积分产生，它不是 exp(ik.r) 的虚部通常在进行解析傅里叶变换时出现。

然而，结果不是动量空间中的高斯函数。

Module FFTW3
 use, intrinsic :: iso_c_binding
include 'fftw3.f03'
end module  

program sine_FFT_transform
use FFTW3
implicit none
integer, parameter :: dp=selected_real_kind(8)

real(kind=dp), parameter :: pi=acos(-1.0_dp)
integer, parameter :: n=1024 
real(kind=dp) :: delta, k
real(kind=dp) :: numerical_F_transform
integer :: i
type(C_PTR) ::  my_plan
real(C_DOUBLE), dimension(1024) :: y
real(C_DOUBLE), dimension(1024) :: yy, yk
integer(C_FFTW_R2R_KIND) :: C_FFTW_RODFT00

my_plan= fftw_plan_r2r_1d(1024,y,yy,FFTW_FORWARD, FFTW_ESTIMATE)

delta=0.0125_dp
do i=1, n        !inserting the input one-dimension position function
y(i)= 2*(delta)*(i-1)*exp(-((i-1)*delta)**2) 
! I multiplied by 2 due to the definition of C_FFTW_RODFT00 in FFTW3
end do

call fftw_execute_r2r(my_plan, y,yy)   
do i=2, n
k = (i-1)*pi/n/delta 
yk(i) = 4*pi*delta*yy(i)/2  !I divide by 2 due to the definition of 
                            !C_FFTW_RODFT00
numerical_F_transform=yk(i)/k
write(11,*) i,k,numerical_F_transform
end do
call fftw_destroy_plan(my_plan)

end program 

执行前面的代码给出了以下不适用于高斯函数的图。

谁能帮我理解问题是什么？我猜这个问题主要是由于 FFTW3 .也许我没有正确使用它，尤其是在边界条件方面。

Answer 1

查看 FFTW 站点中的相关页面(Real-to-Real Transforms、transform kinds、Real-odd DFT (DST))和 Fortran 的头文件，似乎 FFTW 期望 FFTW_RODFT00等等而不是FFTW_FORWARD用于指定类型
真实到真实的转换。例如，

! my_plan= fftw_plan_r2r_1d( n, y, yy, FFTW_FORWARD, FFTW_ESTIMATE )
my_plan= fftw_plan_r2r_1d( n, y, yy, FFTW_RODFT00, FFTW_ESTIMATE )

执行上页所示的“I 型”离散正弦变换 (DST-I)。这种修改似乎解决了这个问题(即，使傅里叶变换成为具有正值的高斯)。

以下是 OP 代码的略微修改版本，用于实验上述修改:

! ... only the modified part is shown...
real(dp) :: delta, k, r, fftw, num, ana
integer :: i, j, n
type(C_PTR) ::  my_plan
real(C_DOUBLE), allocatable :: y(:), yy(:)

delta = 0.0125_dp ; n = 1024   ! rmax = 12.8
! delta = 0.1_dp    ; n = 128    ! rmax = 12.8
! delta = 0.2_dp    ; n = 64    ! rmax = 12.8
! delta = 0.4_dp    ; n = 32    ! rmax = 12.8

allocate( y( n ), yy( n ) )

! my_plan= fftw_plan_r2r_1d( n, y, yy, FFTW_FORWARD, FFTW_ESTIMATE )
my_plan= fftw_plan_r2r_1d( n, y, yy, FFTW_RODFT00, FFTW_ESTIMATE )

! Loop over r-grid
do i = 1, n
    r = i * delta              ! (2-a)
    y( i )= r * exp( -r**2 )
end do

call fftw_execute_r2r( my_plan, y, yy )

! Loop over k-grid
do i = 1, n

    ! Result of FFTW
    k = i * pi / ((n + 1) * delta)    ! (2-b)
    fftw = 4 * pi * delta * yy( i ) / k / 2   ! the last 2 due to RODFT00

    ! Numerical result via quadrature
    num = 0
    do j = 1, n
        r = j * delta
        num = num + r * exp( -r**2 ) * sin( k * r )
    enddo
    num = num * 4 * pi * delta / k

    ! Analytical result
    ana = sqrt( pi )**3 * exp( -k**2 / 4 )

    ! Output
    write(10,*) k, fftw
    write(20,*) k, num
    write(30,*) k, ana
end do

编译(使用 gfortran-8.2 + FFTW3.3.8 + OSX10.11):

$ gfortran -fcheck=all -Wall sine.f90 -I/usr/local/Cellar/fftw/3.3.8/include -L/usr/local/Cellar/fftw/3.3.8/lib -lfftw3

如果我们使用 FFTW_FORWARD如在原始代码中，我们得到



它有一个负瓣(其中 fort.10、fort.20 和 fort.30 对应于 FFTW、正交和分析结果)。修改代码以使用 FFTW_RODFT00将结果更改如下，因此修改似乎正在工作(但请参阅下面的网格定义)。



补充笔记

我在我的代码(第 (2-a) 和 (2-b) 行)中稍微修改了 r 和 k 的网格定义，发现这可以提高准确性。但是我还是不确定上面的定义是否与FFTW使用的定义相符，所以请阅读手册了解详情...

fftw3.f03头文件给出了fftw_plan_r2r_1d的接口(interface)

type(C_PTR) function fftw_plan_r2r_1d(n,in,out,kind,flags) bind(C, name='fftw_plan_r2r_1d')
  import
  integer(C_INT), value :: n
  real(C_DOUBLE), dimension(*), intent(out) :: in
  real(C_DOUBLE), dimension(*), intent(out) :: out
  integer(C_FFTW_R2R_KIND), value :: kind
  integer(C_INT), value :: flags
end function fftw_plan_r2r_1d

(因为没有Tex支持，这部分很丑。。。)4 pi r^2 * exp(-r^2) * sin(kr)/(kr)的积分对于 r = 0 -> 无限是 pi^(3/2) * exp(-k^2 / 4) (从 Wolfram Alpha 获得或注意到这实际上是 exp(-(x^2 + y^2 + z^2)) 的 3-D 傅里叶变换，由 exp(-i*(k1 x + k2 y + k3 z)) 其中 k =(k1,k2,k3))。因此，虽然有点违反直觉，但结果变成了正高斯。

我猜 r-grid 可以选择得更粗一些(例如 delta 高达 0.4)，只要它覆盖变换函数的频域(这里是 exp(-r^2))，它就可以提供几乎相同的精度。