我试图通过使用具有矩形域(-4 ＜＝ x ＜＝ 4和-2 ＜＝ y ＜＝ 2)的FFTW库来求解泊松方程.如果域是正方形，我有正确的结果，如果域是矩形，那是错误的.请给我一些建议.非常感谢.这是我的代码.
I am trying to solve Poisson equation by using FFTW library with rectangular domain(-4<=x<=4 and -2<=y<=2). I have correct result if domain is square and it is wrong if domain is rectangular. Please give me some suggestion. Thank you so much. Here is my code.
#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; } 推荐答案在傅立叶空间中，频率k_x和k_y必须取决于域的大小.由于域是矩形的，因此变得很重要.
In the Fourier space, the frequencies k_x and k_y must depend on the size of the domain. As the domain is rectangular, it becomes a matter of importance.
尝试:
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; }输出应接近预期的输出(比例因子除外).
The output should be closer to the expected one (except for a scale factor).