Лапласовский расчет с fftw в c

Question

Я пытаюсь реализовать лапласиан в C путем вычисления преобразования Фурье (вперед и назад) экспоненциальной функции с пакетом fftw, однако после тестирования результатов с аналитической версией лапласиана я получаю совсем другие результаты. Ниже приведен код:

#include <stdlib.h>
#include <stddef.h>
#include <stdio.h>
#include <math.h>
#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif 

// See for reference:
// https://en.cppreference.com/w/c/numeric/complex
#include <complex.h>

// See for reference:
// http://www.fftw.org/fftw3_doc/
#include <fftw3.h>


/**
 * Test function
 * */
double function_xy(double x, double y)
{
    #define A -0.03
    #define B -0.01
    #define C -0.005 
    return exp(A*x*x + B*y*y + C*x*y);
}

/**
 * Analytical result
 * Use it for checking correctness!
 * */
double laplace_function_xy(double x, double y)
{
    double rdf2d_dx = function_xy(x,y)*(2.*A*x + C*y); // d/dx
    double rdf2d_dy = function_xy(x,y)*(2.*B*y + C*x); // d/dy

    double rlaplacef2d = rdf2d_dx*(2.*A*x + C*y)+function_xy(x,y)*(2.*A)  +  rdf2d_dy*(2.*B*y + C*x)+function_xy(x,y)*(2.*B); // laplace

    return rlaplacef2d;

    #undef A
    #undef B
    #undef C
}


/**
 * You can use this function to check diff between two arrays
 * */
void test_array_diff(int N, double *a, double *b)
{
    int ixyz=0;

    double d,d2;
    double maxd2 = 0.0;
    double sumd2 = 0.0;
    for(ixyz=0; ixyz<N; ixyz++)
    {
        d = a[ixyz]-b[ixyz];

        d2=d*d;
        sumd2+=d2;
        if(d2>maxd2) maxd2=d2;
    }

    printf("#    COMPARISON RESULTS:\n");
    printf("#           |max[a-b]| : %16.8g\n", sqrt(maxd2));
    printf("#         SUM[(a-b)^2] : %16.8g\n", sumd2);
    printf("# SQRT(SUM[(a-b)^2])/N : %16.8g\n", sqrt(sumd2)/N);
}


int main()
{

    // Settings
    int nx = 100; // number of points in x-direction
    int ny = 100; // number of points in y-direction

    int nthreads = 2;

    double Lx = 100.0; // width in x-direction
    double Ly = 100.0; // width in y-direction

    double x0 = -Lx/2;
    double y0 = -Ly/2;

    double dx = Lx/nx;
    double dy = Ly/ny;

    double *fxy; // function
    double *test_laplacefxy; // array with results computed according formula
    fxy = (double *) malloc(nx*ny*sizeof(double)); 
    if(fxy==NULL) { printf("Cannot allocate array fx!\n"); return 0;}

    test_laplacefxy = (double *) malloc(nx*ny*sizeof(double)); 
    if(test_laplacefxy==NULL) { printf("Cannot allocate array fx!\n"); return 0;}

    int ix, iy, ixy;

    ixy=0;
    for(ix=0; ix<nx; ix++) for(iy=0; iy<ny; iy++) 
    {
        // function
        fxy[ixy] = function_xy(x0 + dx*ix, y0 + dy*iy);

        // result for comparion
        test_laplacefxy[ixy] = laplace_function_xy(x0 + dx*ix, y0 + dy*iy);
        ixy++;
    }

    printf("Value of fxy, %2f\n", fxy[40]);


    double kx, ky;
    double *laplacefxy; // pointer to array with laplace

    double complex *fk; // its fourier transform
    fk = (double complex *) malloc(nx*ny*sizeof(double complex)); 
    if(fk==NULL) { printf("Cannot allocate array fk!\n"); return 0;}
    laplacefxy = (double *) malloc(nx*ny*sizeof(double)); 
    if(laplacefxy==NULL) { printf("Cannot allocate array fx!\n"); return 0;}

    // Create fftw plan

    fftw_plan plan_f = fftw_plan_dft_r2c_2d(nx, ny, fxy, fk, FFTW_ESTIMATE);

    // execute the plan

    fftw_execute(plan_f);

    ixy=0;
    for(ix=0; ix<nx; ix++) for(iy=0; iy<ny; iy++) 
    {
        if(ix<nx/2) kx=2.*M_PI/(dx*nx)*(ix   );
        else        kx=2.*M_PI/(dx*nx)*(ix-nx);

        if(iy<ny/2) ky=2.*M_PI/(dy*ny)*(iy   );
        else        ky=2.*M_PI/(dy*ny)*(iy-ny);

        fk[ixy] = (-(kx * kx) - (ky * ky)) * fk[ixy] / (nx * ny);
        ++ixy;
    }

    printf("Value of fk, %2f, fk imag %2f\n", creal(fk[4]), cimag(fk[4]));

     // Create fftw plan
     // execute the plan
    fftw_plan plan_f_backward = fftw_plan_dft_c2r_2d(nx, ny, fk, laplacefxy, FFTW_FORWARD);
    fftw_execute(plan_f_backward);

    // Uncomment to print results

    // for(ix=0; ix<nx; ix++) for(iy=0; iy<ny; iy++) 
    // {
    //     printf("%16.8f %16.8g %16.8g\n", kx, laplacefxy[ixy], laplace_function_xy(x0 + dx*ix, y0 + dy*iy));
    // }


    // Check correctness of computation
    test_array_diff(nx*ny, laplacefxy, test_laplacefxy);

    return 1;
}

Интересно, связано ли это с тем, что я использовал флаг FFTW_ESTIMATE. Я вижу, что в результатах трансармации есть примерно все нули. Буду признателен за любую помощь в этом.