FFTW.jl 用于 2D 阵列：扩散只发生在 1D答案

【问题标题】：FFTW.jl for 2D array: Diffusion only happening in 1DFFTW.jl 用于 2D 阵列：扩散只发生在 1D
【发布时间】：2021-05-13 02:48:22
【问题描述】：

根据我的阅读，当 A 是 2D 数组时，使用 FFTW.jl / AbstractFFTs.jl 的 fft(A) 应该在 2D 中执行 fft，而不是按列执行。知道为什么当（我认为）我向 u(t,x) 添加缩放的二阶空间导数时，我只看到列方向的扩散，就好像使用显式求解器一样？

谢谢！我对此很陌生。

code and heatmap screenshot

using Random
using FFTW
using Plots

gr()

N = (100,100)

# initialize with gaussian noise
u = randn(Float16, (N[1], N[2])).*0.4.+0.4
# include square of high concentration to observe diffusion clearly
u[40:50,40:50] .= 3

N = size(x)
L = 100
k1 = fftfreq(51)
k2 = fftfreq(51)
lap_mat = -(k1.^2 + k2.^2)

function lap_fft(x)
    lapF = rfft(x)
    lap = irfft(lap_mat.*lapF, 100)
    return lap
end

# ode stepper or Implicit-Explicit solver
for i in 1:100000
    u+=lap_fft(u)*0.0001
end
# plot state
heatmap(u)

【问题讨论】：

标签： julia fft fftw differentialequations.jl

【解决方案1】：

仅仅因为您正在执行真正的 FFT，并不意味着您可以对结果进行真正的逆 fft。 rfft 来自 R -> C。但是您可以执行以下操作：

function lap_fft(x)
    lapF = complex(zeros(100,100));        # only upper half filled
    lapF[1:51,1:100] = rfft(x) .* lap_mat; # R -> C
    return abs.(ifft(lapF));               # C -> R
end

真正的 FFT 到复频域（由于数据冗余，只填充了上半部分），在频域中乘以您的滤波器，将 FFT 反演到复数图像域并获得幅度abs.()，实部real.() 等。
但老实说，为什么要使用真正的 fft？

using Random
using FFTW
using Plots

gr()

N = (100,100)

# initialize with gaussian noise
u = randn(Float16, (N[1], N[2])).*0.4.+0.4;
# include square of high concentration to observe diffusion clearly
u[40:50,40:50] .= 3;

N = size(u);
L = 100;
k1 = fftfreq(100);
k2 = fftfreq(100);
tmp = -(k1.^2 + k2.^2);
lap_mat = sqrt.(tmp.*reshape(tmp,1,100));


function lap_fft(x)
    return abs.(ifftshift(ifft(fftshift(ifftshift(fft(fftshift(x))).*lap_mat))));
end

# ode stepper or Implicit-Explicit solver
for i in 1:100000
    u+=lap_fft(u)*0.001;
end
# plot state
heatmap(u)

【讨论】：