在 C 中用 FFTW 求解泊松方程时失去径向对称性

Question

我对使用 FFTW 求解泊松方程有疑问。我要求解的方程式如下：

问题是我得到的解决方案不再是径向对称的，特别是沿轴失去了对称性，正如您从输出图中看到的那样：

代码如下：

    double dx = 25. / ( N*1. ), dy = 25. / ( N*1. ), dkx = 2*pi/(N*dx), dky = 2*pi/(N*dy);//define steps in real and Fourier space
    double Lx = N*dx, Ly = N*dy, Lkx = N*dkx, Lky = N*dky;//define box sides in real and Fourier space
    int i, j, Nh = N/2 + 1;
    double *inr, *in2r, *in2d;
    fftw_complex *outr;
    fftw_complex *outd;
    fftw_plan rplan_forward;
    fftw_plan dplan_backward;
    fftw_plan rplan_backward;

    for ( i = 0; i < N; i++ ) 
    { 
        x[i] = dx*i - Lx/2.;
        kx[i] = dkx*i - Lkx/2.;
        y[i] = dy*i - Ly/2.;
        ky[i] = dky*i - Lky/2.;
    }

    outd = (fftw_complex*)fftw_malloc ( sizeof ( fftw_complex ) * N * Nh );

    inr = ( double * ) malloc ( sizeof ( double ) * N * N );


    for ( i = 0; i < N; i++ ) //Initialize ine to rho
    {
        for ( j = 0; j < N; j++ )
        {
            inr[i*N+j] = rho[i*N+j] * pow( -1, i + j );                            
        }
    }

    outr = (fftw_complex*)fftw_malloc ( sizeof ( fftw_complex ) * N * Nh );

    rplan_forward = fftw_plan_dft_r2c_2d ( N, N, inr, outr, FFTW_ESTIMATE );

    fftw_execute ( rplan_forward ); //fft rho

    for ( i = 0; i < N; i++ ) //Evaluate d in Fourier Space according to Poisson eq
    {
        for ( j = 0; j < Nh; j++ ) 
        {
            outd[i*Nh+j][0] = g * outr[i*Nh+j][0] / ( kx[i]*kx[i] + ky[j]*ky[j] ); 
            outd[i*Nh+j][1] = g * outr[i*Nh+j][1] / ( kx[i]*kx[i] + ky[j]*ky[j] ); 
        }                      
    }

    outd[Nh*N/2+N/2][0] = 0.; //Avoid singularity at k = 0
    outd[Nh*N/2+N/2][1] = 0.;

    in2d = ( double * ) malloc ( sizeof ( double ) * N * N );

    dplan_backward = fftw_plan_dft_c2r_2d ( N, N, outd, in2d, FFTW_ESTIMATE );

    fftw_execute ( dplan_backward ); //Antitrasform d

    in2r = ( double * ) malloc ( sizeof ( double ) * N * N );

    rplan_backward = fftw_plan_dft_c2r_2d ( N, N, outr, in2r, FFTW_ESTIMATE );

    fftw_execute ( rplan_backward ); //Antitrasform rho

    for( i = 0; i < N; i++ ) 
    {
        for ( j = 0; j < N; j++ ) 
        {
            printf( " %le %le %le \n", x[i], y[j], in2d[i*N+j] * pow( -1, i + j ) / ( N*N*1. ) - in2d[0] / ( N*N*1. ) ); //print d in such a way that it is zero at the boundaries
        }

    }

    fftw_destroy_plan ( rplan_forward );
    fftw_destroy_plan ( dplan_backward );
    fftw_destroy_plan ( rplan_backward );
    fftw_free ( outd );
    fftw_free ( inr );
    fftw_free ( outr ); 
    fftw_free ( in2d ); 
    fftw_free ( in2r); 
    fftw_cleanup();

请注意，我想要一个在边界处为零的解。

你能告诉我问题出在哪里吗？我怀疑这只是一个数值问题，但我不确定。

这是边界的放大图（特别是 -12.5 < y < -11.5 区域）。

Answer 1

我认为这里有一个问题：

outd[Nh*N/2+N/2][0] = 0.; //Avoid singularity at k = 0

k=0 在 (0,0)，而不是 (N/2,N/2)。这就是 DFT 的定义方式。 k 从 0 运行到 N-1。 ~~您正在空间域中引入相移。我不知道这会如何影响您以后的计算。~~ 您正试图通过将输入乘以 pow(-1,i) 来使频谱居中。这个技巧仅甚至 N!

如果正确计算 kx 和 ky 会更好：

for ( i = 0; i < N; i++ ) {
   kx[i] = ( i > N/2 ? i-N : i ) * 2*pi/N;
}

这适用于任何 N。另请注意，我上面定义的 k 对应的是径向频率 (ω)，而不是 DFT 的整数 k。它也与你的不同，我认为（如果是这样，你的不正确！）。

另一个问题是DFT不是FT。 DFT 是在有界空间域上计算的，假定是周期性的。这个域是一个正方形。很难在正方形网格上做旋转对称的事情。与轴相比，沿对角线的奈奎斯特频率（最大可能频率）是 sqrt(2) 倍。并且图案的周期也是沿对角线的 sqrt(2) 倍。因此，您会沿轴获得更多混叠，并且您的图案更有可能沿轴重叠。

如果以上原因导致您看到的效果，请尝试增加空间域的大小，and/or 采样密度。

Answer 2

径向对称性的损失与域的形状有关：它是正方形，而不是圆形。此外，由于应用了离散傅里叶变换，因此研究的问题是使用周期性边界条件的方形域上的泊松方程。因此，边界附近的温度可能不会降低到零。

这就是为什么不使用零频率（即源项的平均值）的原因。否则，就无法求解稳态热方程。实际上，由于周期性边界条件，从一侧流出的热通量会从另一侧进入。因为当域完全绝缘时，温度会上升并且永远不会达到稳定的解决方案。

要应用零狄利克雷边界解，必须使用 正弦变换 之一，而不是傅里叶变换。例如参见：

http://www.ds.ewi.tudelft.nl/vakken/in4049TU/sheets/lecture10.pdf

这是一个答案（我的！），其中包含一段代码，通过 FFTW 的正弦变换求解二维泊松方程，零狄利克雷边界条件。它是 c++，但将其转换为 c 很简单。

要恢复接近径向对称的东西，源的支撑必须远离域的边缘。如果您知道域、边界条件和源项遵循径向对称性，最好的解决方案可能是切换到一维径向方程。

在 C 中用 FFTW 求解泊松方程时失去径向对称性

Losing radial symmetry in solving Poisson equation with FFTW in C

c

fft

fftw