用FFTW3库评估FORTRAN中高斯函数的快速傅里叶变换
Moh*_*ssi 8 fortran physics fft fftw
我正在尝试编写FORTRAN代码来评估f(r)=exp(-(r^2))使用FFTW3库的高斯函数的快速傅里叶变换.众所周知,高斯函数的傅里叶变换是另一种高斯函数.
我考虑在球面坐标中评估高斯函数的傅里叶变换积分.
因此,可以将得到的积分简化为积分[r*exp(-(r^2))*sin(kr)]dr.
我编写了以下FORTRAN代码来评估离散SINE变换DST,它是使用PURELY实数输入数组的离散傅立叶变换DFT.DST通过执行C_FFTW_RODFT00现有的FFTW3,考虑到在位置空间中的离散值是R = I*增量(i = 1,2,...,1024),以及用于DST输入数组是函数r*exp(-(r^2))NOT高斯.积分中的正弦函数是由积分[r*exp(-(r^2))*sin(kr)]dr超过SPHERICAL坐标产生的,并且它不是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.也许我没有正确使用它,尤其是边界条件.
查看FFTW站点中的相关页面(Real-to-Real Transforms,转换种类,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 )
执行上页所示的"type-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更改结果如下所示,因此修改似乎有效(但请参阅下面的网格定义).
补充说明
- 我稍微修改了我的代码中的r和k的网格定义(行(2-a)和(2-b)),这被发现可以提高准确性.但我仍不确定上述定义是否与FFTW使用的定义相符,所以请阅读手册了解详情......
在
fftw3.f03头文件中给出了界面fftw_plan_r2r_1d
Run Code Online (Sandbox Code Playgroud)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))通过exp(-i*(k1 x + k2 y + k3 z)),其中k =(k1,k2,k3)).所以,虽然有点反直觉,但结果变成了正高斯.- 我猜r网格可以选择更粗糙(例如
delta高达0.4),只要覆盖变换函数的频域(这里exp(-r^2)),它就能提供几乎相同的精度.