混淆测试 fftw3 - 泊松方程 2d 测试答案

【问题标题】：Confusion testing fftw3 - poisson equation 2d test混淆测试 fftw3 - 泊松方程 2d 测试
【发布时间】：2014-06-02 15:23:28
【问题描述】：

我无法解释/理解以下现象：为了测试 fftw3，我使用了 2d poisson 测试用例：

laplacian(f(x,y)) = - g(x,y) 具有周期性边界条件。

对方程应用傅里叶变换后，我们得到：F(kx,ky) = G(kx,ky) /(kx² + ky²) (1)

如果我取 g(x,y) = sin (x) + sin(y) , (x,y) \in [0,2 \pi] 我立即得到 f(x,y) = g(x ,y)

这是我想用 fft 获得的：

我通过正向傅里叶变换从 g 计算 G

由此我可以用 (1) 计算 f 的傅里叶变换。

最后，我使用反向傅里叶变换计算 f（不要忘记按 1/(nx*ny) 进行归一化）。

在实践中，结果很糟糕？

（例如，N = 256 的幅度是 N = 512 获得的幅度的两倍）

更糟糕的是，如果我尝试 g(x,y) = sin(x)*sin(y) ，曲线甚至没有相同形式的解。

（请注意，我必须更改方程式；在这种情况下，我将拉普拉斯算子除以二：（1）变为 F(kx,ky) = 2*G(kx,ky)/(kx²+ky²)

代码如下：

/*
* fftw test -- double precision
*/
#include <iostream>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <fftw3.h>
using namespace std;

int main()
{
 int N = 128;
 int i, j ;
 double pi = 3.14159265359;
 double *X, *Y  ; 
 X = (double*) malloc(N*sizeof(double));
   Y = (double*) malloc(N*sizeof(double));
   fftw_complex  *out1, *in2, *out2, *in1;
   fftw_plan     p1, p2;
   double L  = 2.*pi;
   double dx = L/(N - 1);

   in1 = (fftw_complex*) fftw_malloc(sizeof(fftw_complex)*(N*N) );
   out2 = (fftw_complex*) fftw_malloc(sizeof(fftw_complex)*(N*N) );
   out1 = (fftw_complex*) fftw_malloc(sizeof(fftw_complex)*(N*N) );
   in2 = (fftw_complex*) fftw_malloc(sizeof(fftw_complex)*(N*N) );

   p1 = fftw_plan_dft_2d(N, N, in1, out1, FFTW_FORWARD,FFTW_MEASURE ); 
   p2 = fftw_plan_dft_2d(N, N, in2, out2, FFTW_BACKWARD,FFTW_MEASURE);

   for(i = 0; i < N; i++){
       X[i] = -pi + i*dx ;
       for(j = 0; j < N; j++){
            Y[j] = -pi + j*dx ;
        in1[i*N + j][0] = sin(X[i]) + sin(Y[j]) ; // row major ordering
        //in1[i*N + j][0] = sin(X[i]) * sin(Y[j]) ; // 2nd test case
        in1[i*N + j][1] = 0 ; 
       }
   }

     fftw_execute(p1); // FFT forward 

     for ( i = 0; i < N; i++){   // f = g / ( kx² + ky² )  
       for( j = 0; j < N; j++){
         in2[i*N + j][0] = out1[i*N + j][0]/ (i*i+j*j+1e-16); 
         in2[i*N + j][1] = out1[i*N + j][1]/ (i*i+j*j+1e-16); 
         //in2[i*N + j][0] = 2*out1[i*N + j][0]/ (i*i+j*j+1e-16); // 2nd test case
         //in2[i*N + j][1] = 2*out1[i*N + j][1]/ (i*i+j*j+1e-16); 
       }
     }

     fftw_execute(p2); //FFT backward

     // checking the results computed

     double erl1 = 0.;
     for ( i = 0; i < N; i++) {
       for( j = 0; j < N; j++){
         erl1 += fabs( in1[i*N + j][0] -  out2[i*N + j][0]/N/N )*dx*dx; 
         cout<< i <<" "<< j<<" "<< sin(X[i])+sin(Y[j])<<" "<<  out2[i*N+j][0]/N/N <<" "<< endl; // > output
        }
      }
      cout<< erl1 << endl ;  // L1 error

      fftw_destroy_plan(p1);
      fftw_destroy_plan(p2);
      fftw_free(out1);
      fftw_free(out2);
      fftw_free(in1);
      fftw_free(in2);

      return 0;
    }

我在我的代码中找不到任何（更多）错误（我上周安装了 fftw3 库），我也没有看到数学问题，但我认为这不是 fft 的错。因此我的困境。我完全没有想法，也完全没有谷歌。

任何帮助解决这个难题将不胜感激。

注意：

编译：g++ test.cpp -lfftw3 -lm

执行：./a.out > 输出

我使用 gnuplot 来绘制曲线：（在 gnuplot 中） splot “输出” u 1:2:4（用于计算的解决方案）

【问题讨论】：

我认为您的 N - 1 术语实际上应该是 N ？
如果我做了“for(i = 0; i
我删除了我的答案：您的公式是正确的...我在尝试重新计算此公式时出错...

标签： c++ fftw

【解决方案1】：

这里有几点需要修改：

您需要考虑所有小频率，包括负频率！索引i 对应于频率2PI i/N，也对应于频率2PI (i-N)/N。在傅立叶空间中，数组的结尾与开头一样重要！在我们的例子中，我们保持最小频率：数组的前半部分是2PI i/N，后半部分是 2PI(i-N)/N。
当然，正如保罗所说，N-1 应该是 Nin double dx = L/(N - 1); => double dx = L/(N ); N-1 不对应于一个连续的周期信号。做个测试用例就很难了……
缩放...我凭经验做到了

对于这两种情况，我获得的结果都更接近预期。这是代码：

    /*
 * fftw test -- double precision
 */
#include <iostream>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <fftw3.h>
using namespace std;

int main()
{
    int N = 128;
    int i, j ;
    double pi = 3.14159265359;
    double *X, *Y  ; 
    X = (double*) malloc(N*sizeof(double));
    Y = (double*) malloc(N*sizeof(double));
    fftw_complex  *out1, *in2, *out2, *in1;
    fftw_plan     p1, p2;
    double L  = 2.*pi;
    double dx = L/(N );


    in1 = (fftw_complex*) fftw_malloc(sizeof(fftw_complex)*(N*N) );
    out2 = (fftw_complex*) fftw_malloc(sizeof(fftw_complex)*(N*N) );
    out1 = (fftw_complex*) fftw_malloc(sizeof(fftw_complex)*(N*N) );
    in2 = (fftw_complex*) fftw_malloc(sizeof(fftw_complex)*(N*N) );

    p1 = fftw_plan_dft_2d(N, N, in1, out1, FFTW_FORWARD,FFTW_MEASURE ); 
    p2 = fftw_plan_dft_2d(N, N, in2, out2, FFTW_BACKWARD,FFTW_MEASURE);

    for(i = 0; i < N; i++){
        X[i] = -pi + i*dx ;
        for(j = 0; j < N; j++){
            Y[j] = -pi + j*dx ;
            in1[i*N + j][0] = sin(X[i]) + sin(Y[j]) ; // row major ordering
            //  in1[i*N + j][0] = sin(X[i]) * sin(Y[j]) ; // 2nd test case
            in1[i*N + j][1] = 0 ; 
        }
    }

    fftw_execute(p1); // FFT forward 

    for ( i = 0; i < N; i++){   // f = g / ( kx² + ky² )  
        for( j = 0; j < N; j++){
            double fact=0;
            in2[i*N + j][0]=0;
            in2[i*N + j][1]=0;
            if(2*i<N){
                fact=((double)i*i);
            }else{
                fact=((double)(N-i)*(N-i));
            }
            if(2*j<N){
                fact+=((double)j*j);
            }else{
                fact+=((double)(N-j)*(N-j));
            }
            if(fact!=0){
                in2[i*N + j][0] = out1[i*N + j][0]/fact;
                in2[i*N + j][1] = out1[i*N + j][1]/fact;
            }else{
                in2[i*N + j][0] = 0;
                in2[i*N + j][1] = 0;
            }
            //in2[i*N + j][0] = out1[i*N + j][0];
            //in2[i*N + j][1] = out1[i*N + j][1];
            //  in2[i*N + j][0] = out1[i*N + j][0]*(1.0/(i*i+1e-16)+1.0/(j*j+1e-16)+1.0/((N-i)*(N-i)+1e-16)+1.0/((N-j)*(N-j)+1e-16))*N*N; 
            //  in2[i*N + j][1] = out1[i*N + j][1]*(1.0/(i*i+1e-16)+1.0/(j*j+1e-16)+1.0/((N-i)*(N-i)+1e-16)+1.0/((N-j)*(N-j)+1e-16))*N*N; 
            //in2[i*N + j][0] = 2*out1[i*N + j][0]/ (i*i+j*j+1e-16); // 2nd test case
            //in2[i*N + j][1] = 2*out1[i*N + j][1]/ (i*i+j*j+1e-16); 
        }
    }

    fftw_execute(p2); //FFT backward

    // checking the results computed

    double erl1 = 0.;
    for ( i = 0; i < N; i++) {
        for( j = 0; j < N; j++){
            erl1 += fabs( in1[i*N + j][0] -  out2[i*N + j][0]/(N*N))*dx*dx; 
            cout<< i <<" "<< j<<" "<< sin(X[i])+sin(Y[j])<<" "<<  out2[i*N+j][0]/(N*N) <<" "<< endl; // > output
            //  cout<< i <<" "<< j<<" "<< sin(X[i])*sin(Y[j])<<" "<<  out2[i*N+j][0]/(N*N) <<" "<< endl; // > output
        }
    }
    cout<< erl1 << endl ;  // L1 error

    fftw_destroy_plan(p1);
    fftw_destroy_plan(p2);
    fftw_free(out1);
    fftw_free(out2);
    fftw_free(in1);
    fftw_free(in2);

    return 0;
}

这段代码远非完美，既不优化也不美观。但它几乎给出了预期的结果。

再见，

【讨论】：

您好，感谢您的回答。就我测试而言，它确实有效，但我不明白为什么会这样^^“。关于缩放，我认为你通过将“事实”乘以 NxN 来处理它。我对“N-1”有点困惑应该是N”，我的观点是，如果你有3分，你只需要2个步骤（但如果我改变你的代码它不起作用，所以我一定有问题）。我也觉得奇怪不考虑负面模式，但是您如何证明将其添加到“事实”中？再次非常感谢您的帮助。
我想我现在明白了为什么 N 而不是 N-1 ：在实际空间采样中，您不能重复第一点和最后一点，这就是我们停在 (N-1)/NL 的原因，那是对的吗？我仍然不知道为什么 ( 1/(ii) + 1/ ( (n-i)*(n-i) ) 有任何链接？( 到目前为止我所读到的关于该主题的所有内容都给了我 (kx²+ky² ) kx = 2 pi / Lx * i
你给了我们正确的公式，我为我的第一个答案道歉：我试图重新计算它，但我犯了一个错误。我更改了答案以提供更好的解决方案（正确的解决方案...）