计算科学 - 如何在 Fortran 中高效准确地计算径向函数的傅里叶变换 - 吾爱随笔录

如何在 Fortran 中高效准确地计算径向函数的傅里叶变换

计算科学正则傅立叶分析傅里叶变换

2021-11-28 15:57:30

正如我的问题所述，我想计算傅里叶变换 $F(q)$ 径向函数 $f(r)$ （定义在 $[0,\infty)$ 并且像指数一样衰减 $\exp(-Ar+b)$ 大体上 $r$ ) 在 Fortran 中尽可能准确。函数值来自数据文件（例如，我可以很容易地通过三次插值对其进行插值，并推断出整个行为 $r$ 已知）。

我在 3D 中使用傅里叶变换的“物理”定义，它给出了（因为 $f$ 是径向的）：

\int d r f (r) \exp (i q r) = 4 π \int_{0}^{\infty} r^{2} \frac{\sin (q r)}{q r} f (r) d r

$\int d\mathbf{r}\,f(\mathbf{r})\exp(i\mathbf{q}\mathbf{r}) = 4\pi\int\limits_0^\infty r^2\frac{\sin(qr)}{qr}f(r)\,dr$

我首先尝试为一些选定的值计算这个积分 $q$ 通过使用 Gauss-Legendre 求积，通过 NAG 例程 D01BCF（D01BCF 链接）生成大约 60 或 100 个横坐标和权重。在高斯勒让德求积的情况下，问题是选择区间 $[0,B]$ 在其上集成。而函数 $f$ 损失 4 到 5 个数量级 $r=10$ 到 $r=20$ （例），选择 $B$ 作为对计算结果的强烈影响......当我将结果与“近乎精确”的计算进行比较时（使用 MATLAB 进行但计算时间很长），我发现实际上这仅适用于的小值 $q$ （大约 5，当我必须处理大到 150 的值时）。Gauss-Laguerre 求积并没有给出更好的结果，可能是因为被积函数的振荡部分。

然后我尝试为某些给定的值计算这个傅里叶变换 $q$ 使用例程 D01ASF （D01ASF 链接）。它是一个“一维求积、自适应、半无限区间、权重函数 $\cos(ωx)$ 或 $\sin(ωx)$ ”，正是我所需要的。如果我输入 10E-5 的绝对误差容限，对于高达 80 或 100的，结果非常令人信服。 $q$ 问题是：我需要更大的 $q$ ，并且傅里叶变换 $F(q)$ 处以 ~10E-6 的幅度振荡。将容差降低到 10E-5 已经需要一些时间，甚至会使整个事情从子程序输出一些错误消息，所以我不知道 10E-6 是否可行。 $q$

因此，我目前想知道尝试用 FFT 计算这个傅里叶变换是否不是一个好主意？我面临的问题是我不知道如何用 FFT 计算径向波函数（而且我什至不知道如何正确使用 FFT，因为变换的定义甚至不一样（指数符号和论点）并且我以前从未使用过它）。

你有想法吗？

编辑：这来自https://stackoverflow.com/q/39127340/2320757，有人建议我把它放在 stackoverflow 上。

编辑：我尝试了 FFT（使用 NAG 库中的例程C06FAF）。值，它工作得很好。我面临的问题是总是有一些恒定的标准化因素需要考虑。我不明白为什么。这个归一化因子随着网格中使用它具有幂律的：归一化因子近似（见图，其中黑线是“精确”傅里叶变换，绿色、品红色、蓝色、红色和黄色值计算的 FFT ）。 $q$ $N$ $F = N^(-0.5) \exp(9.9)$ $N$

2个回答

的以下代码，这似乎以某种方式工作。那么，您能否将您的代码与它进行比较，看看是否存在一些规范化差异？（请注意，下面的fft()没有进行归一化，即不存在，而 NAG 库中的 FFT 例程可能会乘以某个归一化因子。） $f(r) = exp(-r)$ $1/sqrt(N)$

function test( N, rmax )

    dr = rmax / N
    r = linspace( 0.0, rmax - dr, N )       # r[1:N] = 0, dr, 2dr, ..., (N-1)dr

    rexp = r .* exp( - r )                  # rexp[i] = r[i] * exp(-r[i]), i=1:N
    integral = - imag( fft( rexp ) ) * dr   # \int_0^rmax dr sin(qr) r exp(-r)

    file = open( "test.dat", "w" )

    for m = 1 : N

        if m < div( N, 2 )
            q = (2 * pi / rmax) * (m - 1)
        else
            q = -(2 * pi / rmax) * (N - m + 1)
        end

        H = integral[ m ]                   # obtained from FFT
        H_exact = 2 * q / ( q^2 + 1 )^2     # analytical solution

        if 0 <= q <= 200.0
            @printf( file, "%20.10e %20.10e %20.10e\n",
                     q, H, H_exact )
        end
    end

    close( file )
end

test( 10^6, 50.0 )

我一直在使用这些类型的转换。我选择的工具不是完整的 FFT，而是离散傅里叶正弦变换 (DST)。并非所有的傅立叶分析软件包都有，但看起来 NAG 的C06RAF可以解决问题。NAG 提供的示例代码应该是一个很好的起点，但这是我遵循的大纲：

确保我的和数组满足奈奎斯特条件。（*见脚注） $q$ $r$ $\Delta q \ge \frac{\pi}{r_{max}}$
上调用 DST 例程。 $r f(r)$
缩放结果。您至少需要一个因子并且可以凭经验确定其余的因子。是因为您的积分从运行到，而不是到。 $\frac{\Delta r}{2q}$ $1/2$ $0$ $\infty$ $-\infty$ $\infty$

倍的差异。这是因为您的实际域是，但 DST 将转换视为因此您可以免费获得两倍的分辨率. $2$ $r$ $[0, r_{max}]$ $[-r_{max}, r_{max}]$ $q$

其它你可能感兴趣的问题

上一篇pde数值解和非线性方程三步下一篇如何改进非对称矩阵的这种双移位 QR 算法？