Given the definition of DST-I, DCT-I, their relation to the DFT ([1],[2]) as well as the differentiation property of the DFT;
[1]: "A DST-I is exactly equivalent to a DFT of a real sequence that is odd around the zero-th and middle points, scaled by 1/2. For example, a DST-I of N=3 real numbers (a,b,c) is exactly equivalent to a DFT of eight real numbers (0,a,b,c,0,-c,-b,-a) (odd symmetry), scaled by 1/2. "
[2]: "The DCT-I is exactly equivalent (up to an overall scale factor of 2), to a DFT of 2N-2 real numbers with even symmetry. For example, a DCT-I of N=5 real numbers abcde is exactly equivalent to a DFT of eight real numbers abcdedcb (even symmetry), divided by two."
How do we go about calculating the Divergence of a velocity field on a square domain (with DST-I boundary condition)?
. . .
My unsuccessful tackle with this problem thus far:
I generate a 2D Function, similar to [0,a,b,c,0,-c,-b,-a] with Matlab:
N=8;
x5 = linspace(3, -3,N+1); x5 = x5(1:N);
[X1,X2] = ndgrid(x5, x5);
Z= sin(X1) .* sin(X2);
The resulting, somewhat hand-modified output look like:
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
Initialize the above matrix as:
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
The DFT of this(B2D_7) matrix is obtained as:
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]
We note here that this matrix has 0 imaginary component and odd-symmetry in y axis.
We now execute divergence_dft(B2D_6_DFT, B2D_6_DFT, DFT_N), with the same DFT matrix(seen above) for FDX and FDY:
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();
}}
}
When we examine the output of divergence_dft, we notice that the symmetry in Y-axis has been lost. It complies with neither DST [0abc0-c-ba] nor DCT [abcdedcb] pattern and we cannot apply either real transforms to invert the matrix:
[0.000000, 0.000000 i] [0.000000, 0.000000 i] [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, 0.000000 i]
[0.000000, 0.000000 i] [0.000000, 61.595959 i] [-0.000000, -12.566371 i] [0.000000, 8.412147 i] [0.000000, 0.000000 i]
[0.000000, 0.000000 i] [-0.000000, -31.415932 i] [0.000000, 8.412146 i] [-0.000000, -11.319254 i] [0.000000, 0.000000 i]
[0.000000, 0.000000 i] [0.000000, 0.000000 i] [0.000000, 0.000000 i] [0.000000, 0.000000 i] [0.000000, 0.000000 i]
[0.000000, 0.000000 i] [0.000000, -15.707966 i] [0.000000, 1.682429 i] [0.000000, 0.000000 i] [0.000000, 0.000000 i]
[0.000000, 0.000000 i] [0.000000, 20.531986 i] [0.000000, 0.000000 i] [0.000000, -1.682429 i] [0.000000, 0.000000 i]
[0.000000, 0.000000 i] [0.000000, 0.000000 i] [0.000000, -20.531986 i] [0.000000, 15.707962 i] [0.000000, 0.000000 i]