我在Fortran中编写了以下代码，使用r2r（实数到真实）类型的FFTW sin变换来求解泊松方程。在这段代码中，首先我完成了正常的FFTW，即数学函数27*sin(3x)*sin(3y)*sin(3z)的r2c（实数到复数）类型，然后将其除以27（3 * 3 + 3 * 3 + 3 * 3）以计算输入函数的二阶导数。输入函数幅度的三维图表显示其正确的幅度.amplitude along z-axis against x-,y- co-ordinate。
然后，c2r（复数到实数）类型的逆FFTW重新生成输入函数，但幅度现在减小到1，如3-D曲线2所示。这说明我的泊松方程求解器代码在正常FFTW的情况下工作正常。
Program Derivative
! To run this code: gcc dst_3d.c -c -std=c99 && gfortran derivative.f95 dst_3d.o -lfftw3 && ./a.out
implicit none
include "fftw3.f"
integer ( kind = 4 ), parameter :: Nx = 64
integer ( kind = 4 ), parameter :: Ny = 64
integer ( kind = 4 ), parameter :: Nz = 64
real ( kind = 8 ), parameter :: pi=3.14159265358979323846d0
integer ( kind = 4 ) i,j,k
real ( kind = 8 ) Lx,Ly,Lz,dx,dy,dz,kx,ky,kz
real ( kind = 8 ) x(Nx),y(Ny),z(Nz)
real ( kind = 8 ) in_dst(Nx,Ny,Nz),out_dst(Nx,Ny,Nz) ! DST
real ( kind = 8 ) in_k_dst(Nx,Ny,Nz),out_k_dst(Nx,Ny,Nz) ! DST
real ( kind = 8 ) in_dft(Nx,Ny,Nz),out_dft(Nx,Ny,Nz) ! DFT
complex ( kind = 8 ) in_k_dft(Nx/2+1,Ny,Nz),out_k_dft(Nx/2+1,Ny,Nz) ! DFT
integer ( kind = 8 ) plan_forward,plan_backward ! DFT
! System Size.
Lx = 2.0d0*pi; Ly = 2.0d0*pi; Lz = 2.0d0*pi
! Grid Resolution.
dx = Lx/dfloat(Nx); dy = Ly/dfloat(Ny); dz = Lz/dfloat(Nz)
! =================================== INPUT ===========================================
! Initial Profile Details.
do i = 1, Nx
x(i) = dfloat(i-1)*dx
do j = 1, Ny
y(j) = dfloat(j-1)*dy
do k = 1, Nz
z(k) = dfloat(k-1)*dz
in_dst(i,j,k) = 27.0d0*sin(3.0d0*x(i))*sin(3.0d0*y(j))*sin(3.0d0*z(k))
in_dft(i,j,k) = in_dst(i,j,k)
write(10,*) x(i), y(j), z(k), in_dft(i,j,k)
enddo
enddo
enddo
! =================================== 3D DFT ===========================================
call dfftw_plan_dft_r2c_3d_ (plan_forward, Nx, Ny, Nz, in_dft, in_k_dft, FFTW_ESTIMATE)
call dfftw_execute_ (plan_forward)
call dfftw_destroy_plan_ (plan_forward)
do i = 1, Nx/2+1
do j = 1, Ny/2
do k = 1, Nz/2
kx = 2.0d0*pi*dfloat(i-1)/Lx
ky = 2.0d0*pi*dfloat(j-1)/Ly
kz = 2.0d0*pi*dfloat(k-1)/Lz
out_k_dft(i,j,k) = in_k_dft(i,j,k)/(kx*kx+ky*ky+kz*kz)
enddo
do k = Nz/2+1, Nz
kx = 2.0d0*pi*dfloat(i-1)/Lx
ky = 2.0d0*pi*dfloat(j-1)/Ly
kz = 2.0d0*pi*dfloat((k-1)-Nz)/Lz
out_k_dft(i,j,k) = in_k_dft(i,j,k)/(kx*kx+ky*ky+kz*kz)
enddo
enddo
do j = Ny/2+1, Ny
do k = 1, Nz/2
kx = 2.0d0*pi*dfloat(i-1)/Lx
ky = 2.0d0*pi*dfloat((j-1)-Ny)/Ly
kz = 2.0d0*pi*dfloat(k-1)/Lz
out_k_dft(i,j,k) = in_k_dft(i,j,k)/(kx*kx+ky*ky+kz*kz)
enddo
do k = Nz/2+1, Nz
kx = 2.0d0*pi*dfloat(i-1)/Lx
ky = 2.0d0*pi*dfloat((j-1)-Ny)/Ly
kz = 2.0d0*pi*dfloat((k-1)-Nz)/Lz
out_k_dft(i,j,k) = in_k_dft(i,j,k)/(kx*kx+ky*ky+kz*kz)
enddo
enddo
enddo
out_k_dft(1,1,1) = in_k_dft(1,1,1)
call dfftw_plan_dft_c2r_3d_ (plan_backward, Nx, Ny, Nz, out_k_dft, out_dft, FFTW_ESTIMATE)
call dfftw_execute_ (plan_backward)
call dfftw_destroy_plan_ (plan_backward)
do i = 1, Nx
do j = 1, Ny
do k = 1, Nz
out_dft(i,j,k) = out_dft(i,j,k)/dfloat(Nx*Ny*Nz)
write(20,*) x(i), y(j), z(k), out_dft(i,j,k)
enddo
enddo
enddo
! =================================== 3D DST ===========================================
call Forward_FFT(Nx, Ny, Nz, in_dst, in_k_dst)
do k = 1, Nz
do j = 1, Ny/2
do i = 1, Nx/2
kz = 2.0d0*pi*dfloat((k-1))/Lz
ky = 2.0d0*pi*dfloat((j-1))/Ly
kx = 2.0d0*pi*dfloat((i-1))/Lx
out_k_dst(i,j,k) = in_k_dst(i,j,k)/(kx*kx+ky*ky+kz*kz)
enddo
do i = Nx/2+1, Nx
kz = 2.0d0*pi*dfloat((k-1))/Lz
ky = 2.0d0*pi*dfloat((j-1))/Ly
kx = 2.0d0*pi*dfloat(Nx-(i-1))/Lx
out_k_dst(i,j,k) = in_k_dst(i,j,k)/(kx*kx+ky*ky+kz*kz)
enddo
enddo
do j = Ny/2+1, Ny
do i = 1, Nx/2
kz = 2.0d0*pi*dfloat((k-1))/Lz
ky = 2.0d0*pi*dfloat(Ny-(j-1))/Ly
kx = 2.0d0*pi*dfloat((i-1))/Lx
out_k_dst(i,j,k) = in_k_dst(i,j,k)/(kx*kx+ky*ky+kz*kz)
enddo
do i = Nx/2+1, Nx
kz = 2.0d0*pi*dfloat((k-1))/Lz
ky = 2.0d0*pi*dfloat(Ny-(j-1))/Ly
kx = 2.0d0*pi*dfloat(Nx-(i-1))/Lx
out_k_dst(i,j,k) = in_k_dst(i,j,k)/(kx*kx+ky*ky+kz*kz)
enddo
enddo
enddo
out_k_dst(1,1,1) = in_k_dst(1,1,1)
call Backward_FFT(Nx, Ny, Nz, out_k_dst, out_dst)
do i = 1, Nx
do j = 1, Ny
do k = 1, Nz
out_dst(i,j,k) = out_dst(i,j,k)/dfloat(8*Nx*Ny*Nz)
write(30,*) x(i), y(j), z(k), out_dst(i,j,k)
enddo
enddo
enddo
end program Derivative
! ================================== FFTW SUBROUTINES
====================================================
! ================================================================= !
! Wrapper Subroutine to call forward fftw c functions for 3d arrays !
! ================================================================= !
subroutine Forward_FFT(Nx, Ny, Nz, in, out)
implicit none
integer ( kind = 4 ) Nx,Ny,Nz,i,j,k
real ( kind = 8 ) in(Nx, Ny, Nz),out(Nx, Ny, Nz),dum(Nx*Ny*Nz)
do i = 1, Nx
do j = 1, Ny
do k = 1, Nz
dum(((i-1)*Ny+(j-1))*Nz+k) = in(i,j,k)
enddo
enddo
enddo
call dst3f(Nx, Ny, Nz, dum)
do i = 1, Nx
do j = 1, Ny
do k = 1, Nz
out(i,j,k) = dum(((i-1)*Ny+(j-1))*Nz+k)
enddo
enddo
enddo
end subroutine
! ================================================================== !
! Wrapper Subroutine to call backward fftw c functions for 3d arrays !
! ================================================================== !
subroutine Backward_FFT(Nx, Ny, Nz, in, out)
implicit none
integer ( kind = 4 ) Nx,Ny,Nz,i,j,k
real ( kind = 8 ) in(Nx, Ny, Nz),out(Nx, Ny, Nz),dum(Nx*Ny*Nz)
do i = 1, Nx
do j = 1, Ny
do k = 1, Nz
dum(((i-1)*Ny+(j-1))*Nz+k) = in(i,j,k)
enddo
enddo
enddo
call dst3b(Nx, Ny, Nz, dum)
do i = 1, Nx
do j = 1, Ny
do k = 1, Nz
out(i,j,k) = dum(((i-1)*Ny+(j-1))*Nz+k)
enddo
enddo
enddo
end subroutine
! ==================================================================
此代码使用下面的C-wrapper来计算前向3D FFTW正弦变换和后向3D FFTW正弦变换，
#include <fftw3.h>
int dst3f_(int *n0, int *n1, int *n2, double *in3cs)
{
double *out3cs;
out3cs = (double*) fftw_malloc(sizeof(double) * (*n0) * (*n1) * (*n2));
fftw_plan p3cs;
p3cs = fftw_plan_r2r_3d(*n0, *n1, *n2, in3cs, out3cs, FFTW_RODFT10, FFTW_RODFT10, FFTW_RODFT10, FFTW_ESTIMATE);
fftw_execute(p3cs);
fftw_destroy_plan(p3cs);
for(int i0=0;i0<*n0;i0++) {
for(int i1=0;i1<*n1;i1++) {
for(int i2=0;i2<*n2;i2++) {
in3cs[(i0*(*n1)+i1)*(*n2)+i2] = out3cs[(i0*(*n1)+i1)*(*n2)+i2];
}
}
}
return 0;
}
int dst3b_(int *n0, int *n1, int *n2, double *in3cs)
{
double *out3cs;
out3cs = (double*) fftw_malloc(sizeof(double) * (*n0) * (*n1) * (*n2));
fftw_plan p3cs;
p3cs = fftw_plan_r2r_3d(*n0, *n1, *n2, in3cs, out3cs, FFTW_RODFT01, FFTW_RODFT01, FFTW_RODFT01, FFTW_ESTIMATE);
fftw_execute(p3cs);
fftw_destroy_plan(p3cs);
for(int i0=0;i0<*n0;i0++) {
for(int i1=0;i1<*n1;i1++) {
for(int i2=0;i2<*n2;i2++) {
in3cs[(i0*(*n1)+i1)*(*n2)+i2] = out3cs[(i0*(*n1)+i1)*(*n2)+i2];
}
}
}
return 0;
}
然后我尝试使用正弦变换FFTW即r2r（真实到真实）类型来解决相同的泊松方程。当我在这里显示3的三维输出图时，幅度现在减小到小于1.我无法找出代码中的错误在哪里，因为在正弦变换的情况下这个幅度减小了。
使用实数到实数变换来求解泊松方程是非常有吸引力的，因为它允许使用各种边界条件。然而，感测点并不明显对应于考虑周期性边界条件的感测点。
对于周期性边界条件，感测点位于规则网格上，例如0,1 ..，n-1。如果单位晶胞的尺寸是Lx，那么这些点之间的间距是Lx / n。故事的结局。
现在，让我们考虑边界条件，使得DST III用于正向变换，即标志RODFT01。如there所示，它在documentation of FFTW发出信号：
FFTW_RODFT01（DST-III）：j = -1附近的奇数，甚至j = n-1附近。
传感点仍为0,1 ..，n-1。如果长度是Lx，则间距仍为dx = Lx /（n）。但是DST III的输入函数和基函数在j = -1附近是奇数，甚至在j = n-1附近。这解释了幅度的差异：
in_dst(i,j,k) = 27.0d0*sin(3.0d0*x(i))*sin(3.0d0*y(j))*sin(3.0d0*z(k))
实际上，这个输入在i = -1附近并不是奇数，甚至在i = n-1附近。 i = 0和i = n都是奇数。因此，以下操作可能会有所帮助：
FFTW_RODFT00 (DST-I): odd around j=-1 and odd around j=n.。最后，可能需要对输出进行评分，因为FFTW transforms are not normalized.对于DST-I，FFTW_RODFT00，它是N=2(n+1)。弗拉基米尔的建议无疑是一个很好的建议。实际上，虽然测试单个频率是理解和实现算法的理想选择，但最终测试必须涵盖多个频率以确保程序可靠。