利用傅立叶正弦变换进行矢量场发散计算

编程问答 2022-06-12

问题描述

鉴于DST-I，DCT-I的定义，它们与DFT的关系（[1]，[2]）以及DFT的区别property；

[1]：“ DST-I完全等于实数序列的DFT，该实数序列在第零和中间点周围是奇数，缩放比例为1/2。例如，DST-I为N = 3个实数（a，b，c）完全等于八个实数（0，a，b，c，0，-c，-b，-a）的DFT（奇对称），​​按1/2缩放。

[2]：“ DCT-I完全等效（最高为2的比例因子），具有对称对称的2N-2个实数的DFT。例如，DCT-I的N = 5实数abcde完全等于八个实数abcdedcb（偶数对称）除以二的DFT。“

如何计算平方域（具有DST-I边界条件）上的速度场发散？

。 。

到目前为止，我未能成功解决此问题：

我使用Matlab生成2D函数，类似于[0，a，b，c，0，-c，-b，-a]：

N=8;
x5 =  linspace(3,-3,N+1); x5 = x5(1:N);
[X1,X2] = ndgrid(x5,x5); 
Z= sin(X1) .* sin(X2);

经过手工修改的输出结果如下：

0   0   0   0   0   0   0   0
0   a   b   c   0   -c  -b  -a
0   b   e   f   0   -f  -e  -b
0   c   f   g   0   -g  -f  -c
0   0   0   0   0   0   0   0
0   -c  -f  -g  0   g   f   c
0   -b  -e  -f  0   f   e   b
0   -a  -b  -c  0   c   b   a

将以上矩阵初始化为：

0   0   0   0   0   0   0   0
0   1   2   3   0   -3  -2  -1
0   2   4   5   0   -5  -4  -2
0   3   5   6   0   -6  -5  -3
0   0   0   0   0   0   0   0
0   -3  -5  -6  0   6   5   3
0   -2  -4  -5  0   5   4   2
0   -1  -2  -3  0   3   2   1

此（B2D_7）矩阵的DFT为：

fftwf_plan gplan_fw= fftwf_plan_dft_r2c_2d(DFT_N,DFT_N,B2D_7,(fftwf_complex*)B2D_6_DFT,FFTW_ESTIMATE);
fftwf_execute(gplan_fw); 


[0.000000]  [0.000000]      [0.000000]  [0.000000]  [0.000000]
[0.000000]  [-81.597977]    [26.142136] [-9.999999] [0.000000]
[0.000000]  [26.142136]     [-4.000000] [2.142136]  [0.000000]
[0.000000]  [-10.000002]    [2.142135]  [-2.402020] [0.000000]
[0.000000]  [0.000000]      [0.000000]  [0.000000]  [0.000000]
[0.000000]  [10.000002]     [-2.142135] [2.402020]  [0.000000]
[0.000000]  [-26.142136]    [4.000000]  [-2.142136] [0.000000]
[0.000000]  [81.597977]     [-26.142136][9.999999]  [0.000000]

我们在这里注意到，这个矩阵在y轴上具有0个虚构分量和奇对称性。

我们现在使用与FDX和FDY相同的DFT矩阵（如上所示）执行divergence_dft（B2D_6_DFT，B2D_6_DFT，DFT_N）：

static void divergence_dft  (float* const FDX,float* const FDY,cuint& N)
{
cfloat fM_PI= (float) M_PI;

cmplx *CMFDX= (cmplx*)  &FDX[0];
cmplx *CMFDY= (cmplx*)  &FDY[0];
const cmplx iota(0,1);

fComplx * const CFDX= (fComplx*) &FDX[0];
fComplx * const CFDY= (fComplx*) &FDY[0];

cfloat rcpL=1.f/N;
cfloat fN=(float) N;
cuint hNZ1 = N/2+1;
cuint hN = N/2;
cfloat fhN= (float) N/2;

uint i,j;

float fcX,fcY;
float nfcX=0.f,nfcY=0.f;

#define fPST 0.0f

for(j=0,fcY=fPST; j<N; ++j,fcY++ ){
    const uint idJ=(j*hNZ1);
    if (fcY<= fhN) nfcY=fcY; else nfcY=fcY-fN ;
                                                 

    for(i=0,fcX=fPST; i<hNZ1; ++i,++fcX ){     
    const uint idc = i+idJ;


        cfloat kwtX=2.f*fM_PI*rcpL* fcX;
        cfloat kwtY=2.f*fM_PI*rcpL* nfcY;

        cmplx deltaFDX= CMFDX[idc] * iota* kwtX;
        cmplx deltaFDY= CMFDY[idc] * iota* kwtY;
        cmplx dotFD= (deltaFDX+deltaFDY);

        CFDX[idc].r = dotFD.real(); 
        CFDX[idc].i = dotFD.imag();
    

    }}


}

当我们检查divergence_dft的输出时，我们注意到Y轴上的对称性已丢失。它既不符合DST [0abc0-c-ba]模式，也不符合DCT [abcdedcb]模式，因此我们无法应用任何实数变换来反转矩阵：

[0.000000,0.000000 i]  [0.000000,0.000000 i]      [0.000000,0.000000 i]
[0.000000,0.000000 i]  [-0.000000,-128.173813 i]  [0.000000,61.595959 i] [-0.000000,-31.415924 i]   [0.000000,61.595959 i]     [-0.000000,-12.566371 i]   [0.000000,8.412147 i]  [0.000000,-31.415932 i]   [0.000000,8.412146 i]  [-0.000000,-11.319254 i]   [0.000000,-15.707966 i]    [0.000000,1.682429 i]  [0.000000,20.531986 i]     [0.000000,-1.682429 i] [0.000000,-20.531986 i]    [0.000000,15.707962 i] [0.000000,0.000000 i]

解决方法

暂无找到可以解决该程序问题的有效方法，小编努力寻找整理中！

如果你已经找到好的解决方法，欢迎将解决方案带上本链接一起发送给小编。

小编邮箱:dio#foxmail.com (将#修改为@）

cdifferential-equationsfftfftfftwmatlabmatlab