FFTW.jl 对于二维阵列：扩散仅发生在一维中

Question

根据我的阅读，当 A 是二维数组时使用 FFTW.jl / AbstractFFTs.jl 的 fft(A) 应该在二维中执行 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)

Answer 1

仅仅因为您正在执行真正的 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)

FFTW.jl 对于二维阵列：扩散仅发生在一维中

FFTW.jl for 2D array: Diffusion only happening in 1D

fft

fftw

julia

differentialequations.jl