从 FFTW 获得的真实数据的一阶、二阶、三阶和四阶偏导数不正确

Question

我计算了 sin(2.0*pi*j/lambda)（j 从 0 到 255 变化）的一阶、二阶和三阶导数，方法是使用 fftw、对其进行傅里叶变换、乘以适当的值并进行傅里叶逆变换。使用的计划不合适，带有虚拟输入和输出阵列。结果数组的实部是分析结果的一半。我通过选择适当的 kx/ky 值来处理混叠问题。只需取 fftw 和逆 fftw 就可以正确返回原始数组。

fftw 返回的值是所有导数解析值的一半。

如果我在 x 方向上做同样的事情（i 而不是 j），则获得的值仅对于 x 的一阶导数是正确的，y 导数不为零（应该为零，因为变化只是在 x) 中，拉普拉斯算子和双谐波（x 和 y 的四阶导数）是错误的。

在这两种情况下，导数的虚部都是非零且很大的。

我使用 r2c 和 c2r 尝试了不同的计划，但这也给出了错误的值。

任何建议都会有所帮助！

    void psiderivcalc(double *psi, double *dx, double *dy, double *del2, double *dxdel2, double *dydel2, double *del4)

{ int i, j, ir, ir1, k;

fftw_complex *fdx, *fdy, *fdxdel2, *fdydel2, *fdel2, *fdel4, *out, *bt, *in;
fftw_plan p, q; 

out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) *N);
in = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) *N);
bt = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) *N);
fdx = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N);
fdy = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N);
fdxdel2 = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) *N);
fdydel2 = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) *N);
fdel2 = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) *N);
fdel4 = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N);

p = fftw_plan_dft_2d(lx, ly, in, bt, FFTW_FORWARD,FFTW_MEASURE);
q = fftw_plan_dft_2d(lx, ly, bt, out, FFTW_BACKWARD,FFTW_MEASURE);

for(i = 0; i < lx; i++)
{ 
    for(j = 0; j < ly; j++)
    { 
        ir = j*lx + i; 
        in[ir] = psi[ir] + 0.0*I; 
    }
} 

fftw_execute(p); 

double kx, ky, pref; 
int y1, x1, y2, x2; 
kx = 2.0*pi/lx; 
ky = 2.0*pi/ly;
for(i=0;i<lx;i++)
{
    for(j = 0; j < ly; j++)
    {
        ir = j*lx + i;
        if(i<lx/2) 
        {
            x1 = i; 
            x2 = i;
        } 
        else if(i==lx/2.0)
        {
            x1 = 0; 
            x2 = i;
        } 
        else
        { 
            x1 = i - lx; 
            x2 = x1;
        }  
        if(j<ly/2) 
        {
            y1 = i; 
            y2 = i;
        } 
        else if(j==ly/2) 
        {
            y1 = 0; 
            y2 = j;
        } 
        else
        { 
            y1 = j - ly; 
            y2 = y1;
        }
        pref = -1.0*(pow((kx*x2),2.0) + pow((ky*y2),2.0));
        fdx[ir] = -1.0*kx*x1*I*bt[ir]; 
        fdy[ir] = -1.0*ky*y1*I*bt[ir]; 
        fdel2[ir] = pref*bt[ir]; 
        fdxdel2[ir] = -1.0*kx*I*x1*pref*bt[ir]; 
        fdydel2[ir] = -1.0*ky*I*y1*pref*bt[ir]; 
        fdel4[ir] = pref*pref*bt[ir];
    }
}

for(i = 0;i < lx; i++)
{ 
    for(j = 0;j < ly; j++)
    { 
        ir = j*lx + i; 
        bt[ir] = fdx[ir]; 
    }
}

fftw_execute(q); 

for(i = 0;i < lx; i++)
{ 
    for(j = 0;j < ly; j++)
    { 
        ir = j*lx + i; 
        dx[ir] = creal(out[ir]) / (double)N; 
    }
}

我通过每次在执行计划 q 之前更改 bt[ir] 和在执行之后更改 out[ir] 找到了与 dx[ir](1st导数 wrt x) 相同的其他导数。 del2 是拉普拉斯算子，dxdel2 是拉普拉斯算子 wrt x 的导数，del4 是双谐波。

psi[ir] 在传递给函数之前已正确初始化为 j 中的正弦函数。

    for(i = 0; i < lx; i++)
    { 
        for(j = 0; j < ly; j++)
        {
            ir = j*lx + i; 
            psi[ir] = sin(2.0*pi*j/lambda); 
        }
    }  

这是 del2(fftw)、del2a(analytical)、dx、dxa 和 ir(index) 的输出示例。 lambda = 8, lx = 16 and ly = 16, N = lx*ly, 所以 sin 函数在这个域中是周期性的。

    0.308425    -0.616850   0.000000    0.000000    175
    -0.218090   -0.436179   0.000000    0.000000    191
    -0.000000   -0.000000   0.000000    0.000000    207
    0.218090    0.436179    0.000000    0.000000    223
    0.308425    0.616850    0.000000    0.000000    239
    0.218090    0.436179    0.000000    0.000000    255

Answer 1

问题在于以下几行：

    if(j<ly/2) 
    {
        y1 = i; // should be j !!!
        y2 = i; // should be j !!!
    } 

另外，观察导数的符号：

    fdx[ir] = -1.0*kx*x1*I*bt[ir]; // should be a + !
    fdy[ir] = -1.0*ky*y1*I*bt[ir]; // should be a + !