使用具有矩形域的 FFTW 的泊松方程

Question

我正在尝试通过使用具有矩形域（-4

#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;
        }

输出应该更接近预期的输出（比例因子除外）。