利用MKL库FFT求导
在利用mkl提供的fft模块计算周期性函数的导数时进行测试1.利用傅里叶变换的微分性质:\mathscr{F}(f‘(x))=i\omega \mathscr{F}(f(x))
$$f'(x)=\mathscr{F}^{-1}\{\mathscr{F}\}=\mathscr{F}^{−1}\{i\omega F\}$$
2. 调用MKL库中的FFT模块,具体调用过程可以参考oneAPI官方手册,代码如下
module fft
use MKL_DFTI
implicit none
real(8), parameter :: pi = dacos(-1.0D0)
contains
subroutine fft1d(n, x_in, x_out)
implicit none
integer, intent(in) :: n
complex(8), intent(in) :: x_in(n)
complex(8), intent(out) :: x_out(n)
type(DFTI_DESCRIPTOR), pointer :: My_Desc_Handle
Integer :: Status
Status = DftiCreateDescriptor( My_Desc_Handle, DFTI_DOUBLE, DFTI_COMPLEX, 1, n )
Status = DftiSetValue( My_Desc_Handle, DFTI_PLACEMENT, DFTI_NOT_INPLACE)
Status = DftiCommitDescriptor( My_Desc_Handle )
Status = DftiComputeForward( My_Desc_Handle, x_in, x_out )
Status = DftiFreeDescriptor(My_Desc_Handle)
end subroutine
subroutine ifft1d(n, x_in, x_out)
implicit none
integer, intent(in) :: n
complex(8), intent(in) :: x_in(n)
complex(8), intent(out) :: x_out(n)
type(DFTI_DESCRIPTOR), pointer :: My_Desc_Handle
Integer :: Status
Status = DftiCreateDescriptor( My_Desc_Handle, DFTI_DOUBLE,DFTI_COMPLEX, 1, n )
Status = DftiSetValue( My_Desc_Handle, DFTI_PLACEMENT, DFTI_NOT_INPLACE)
Status = DftiSetValue( My_Desc_Handle, DFTI_BACKWARD_SCALE, 1.0D0/N)
Status = DftiCommitDescriptor( My_Desc_Handle )
Status = DftiComputeBackward( My_Desc_Handle, x_in, x_out )
Status = DftiFreeDescriptor(My_Desc_Handle)
end subroutine
subroutine fft_1st_deriv_1d(N, x_in, deriv_x, L)
implicit none
integer, intent(in) :: N
complex(8), intent(in) :: x_in(N)
complex(8), intent(out) :: deriv_x(N)
real(8), intent(in) :: L
complex(8) :: fft_x(N), dfft_x(N)
integer :: k
real(8) :: kx(N)
call fft1d(N, x_in, fft_x)
do k = 1, N
if(mod(N,2) .eq. 0) then
write(*,*) N
if( k .eq. N/2+1) then
kx(k) = 2*pi/L*0
dfft_x(k) = cmplx(0.0D0, 1.0D0,kind=8)*kx(k)*fft_x(k)
else if(k .le. N/2) then
kx(k) = 2*pi/L*(k-1)
dfft_x(k) = cmplx(0.0D0, 1.0D0,kind=8)*kx(k)*fft_x(k)
else
kx(k) = 2*pi/L*(k-1-N)
dfft_x(k) = cmplx(0.0D0, 1.0D0,kind=8)*kx(k)*fft_x(k)
end if
else
if(k .le. N/2) then
dfft_x(k) = 2*pi*dcmplx(0.0D0, 1.0D0)/L*(k-1)*fft_x(k)
else
dfft_x(k) = 2*pi*dcmplx(0.0D0, 1.0D0)/L*(k-N-1)*fft_x(k)
end if
end if
end do
call ifft1d(N, dfft_x, deriv_x)
open(unit = 16, file="dfft.txt")
end subroutine
end module fft
programtest
use fft
implicit none
integer, parameter :: n =128
complex(8) :: x_out(n), z(n), deriv_x(n), x_in(n)
! real(8), parameter :: pi = 4.0D0*datan(1.0D0)
real(8) :: l =2*pi
integer(8) :: i
real(8) :: dt
write(*,*) pi
do i = 1, n
! x_in(i) = exp(dcmplx(0.0D0,2.0D0*pi)*real(i)/n)
x_in(i) = cmplx(cos(2*pi*(i-1)/(n-1)),kind = 8)
end do
call fft1d(n,x_in,x_out)
call ifft1d(n,x_out,z)
call fft_1st_deriv_1d(n,x_in, deriv_x, l)
end programtest
在用此办法求导过程中,改变输入函数x_in(i)时发现:
x_in(i) = cmplx(sin(2*pi*(i-1)/(n-1)),kind = 8)
计算得到的导数两端会有振荡,如下图所示:
(https://pic1.zhimg.com/80/v2-b27a2f0d3c486fe9601981dd6dcc61eb_720w.png)
sinx的导数,红色为FFT计算得到的结果,黑色为cosx
但当输入函数改变为:
x_in(i) = cmplx(sin(2*pi*(i-1)/n),kind = 8)
振荡消失,结果如下图所示:
(https://pic1.zhimg.com/80/v2-2968068daa509021c62108b06b788d6e_720w.png)
输入函数为cosx没有遇到这一问题。结果均收敛。但是结果对比matlab的计算和python中numpy模块的fft精度上略有不足。所以这个是什么原因导致的?
cmplx函数有一个bug,反常识,就是不管参数是单精度还是双进度,cmplx(a,b) 的结果都是单精度,cmplx(a,b,8) 才是双精度。 li913 发表于 2022-11-24 22:31
cmplx函数有一个bug,反常识,就是不管参数是单精度还是双进度,cmplx(a,b) 的结果都是单精度,cmplx(a,b,8 ...
还真是这样,妈呀,吓得我赶快把以前的project全局搜索了一遍,发现了好多处cmplx,赶快都改成了dcmplx li913 发表于 2022-11-24 22:31
cmplx函数有一个bug,反常识,就是不管参数是单精度还是双进度,cmplx(a,b) 的结果都是单精度,cmplx(a,b,8 ...
貌似不是这个问题,我调用的时候都是用的cmplx(a,b,kind=8),而且求导后精度问题也很麻烦。
页:
[1]