使用FFTW與矩形域的泊松方程

Question

我想通過使用具有矩形域的FFTW庫來求解泊松方程（-4 < = x < = 4和-2 < = y < = 2）。如果域是平方的，我有正確的結果，如果域是矩形的，那麼它是錯誤的。請給我一些建議。非常感謝。 這是我的代碼。

#include <stdio.h> 
#include <math.h> 
#include <cmath> 
#include <fftw3.h> 
#include <iostream> 
#include <vector> 

using namespace std; 
int main() { 

    int N1=64; 
    int N2=32; 

    double pi = 3.141592653589793; 
    double L1 = 8.0; 
    double dx = L1/(double)(N1-1); 
    double L2= 4.0; 
    double dy=L2/(double)(N2-1); 

    std::vector<double> in1(N1*N2,0.0); 
    std::vector<double> in2(N1*N2,0.0); 
    std::vector<double> out1(N1*N2,0.0); 
    std::vector<double> out2(N1*N2,0.0); 
    std::vector<double> X(N1,0.0); 
    std::vector<double> Y(N2,0.0); 


    fftw_plan p, q; 
    int i,j; 
    p = fftw_plan_r2r_2d(N1,N2, in1.data(), out1.data(), FFTW_REDFT00, FFTW_REDFT00, FFTW_EXHAUSTIVE); 
    q = fftw_plan_r2r_2d(N1,N2, in2.data(), out2.data(), FFTW_REDFT00, FFTW_REDFT00, FFTW_EXHAUSTIVE); 

    int l=-1; 
    for(i = 0;i <N1;i++){ 
     X[i] =-4.0+(double)i*dx ;   
     for(j = 0;j<N2;j++){ 
      l=l+1; 
      Y[j] =-2.0+ (double)j*dy ; 
      in1[l]= cos(pi*X[i]) + cos(pi*Y[j]) ; // row major ordering 
     } 
    } 

    fftw_execute(p); 

    l=-1; 
    for (i = 0; i < N1; i++){ // f = g/(kx² + ky²) 
     for(j = 0; j < N2; j++){ 
       l=l+1; 
      double fact=0; 
      in2[l]=0; 

      if(2*i<N1){ 
       fact=((double)i*i); 
      }else{ 
       fact=((double)(N1-i)*(N1-i)); 
      } 
      if(2*j<N2){ 
       fact+=((double)j*j); 
      }else{ 
       fact+=((double)(N2-j)*(N2-j)); 
      } 
      if(fact!=0){ 
       in2[l] = out1[l]/fact; 
      }else{ 
       in2[l] = 0.0; 
      } 
     } 
    } 

    fftw_execute(q); 
    l=-1; 
    double erl1 = 0.; 
    for (i = 0; i < N1; i++) { 
     for(j = 0; j < N2; j++){ 
      l=l+1; 

      erl1 += fabs(in1[l]- 8.*out2[l]/((double)(N1-1))/((double)(N2-1))); 
      printf("%3d %10.5f %10.5f\n", l, in1[l], 8.*out2[l]/((double)(N1-1))/((double)(N2-1))); 

     } 
    } 

    cout<<"error=" <<erl1 <<endl ; 
    fftw_destroy_plan(p); fftw_destroy_plan(q); fftw_cleanup(); 



    return 0; 
} 

Answer 1

在傅立葉空間中，頻率k_x和k_y必須取決於域的大小。由於該領域是矩形的，這就成爲一個重要問題。

嘗試：

double invL1s=1.0/(L1*L1); 
double invL2s=1.0/(L2*L2); 
... 
     if(2*i<N1){ 
      fact=((double)i*i)*invL1s; 
     }else{ 
      fact=((double)(N1-i)*(N1-i))*invL1s; 
     } 
     if(2*j<N2){ 
      fact+=((double)j*j)*invL2s; 
     }else{ 
      fact+=((double)(N2-j)*(N2-j))*invL2s; 
     } 

輸出應更接近所預期的（除了一個比例因子）。