Быстрое косинус-преобразование Фурье для уравнения Пуассона со всеми нулевыми граничными условиями Неймана

Я пытаюсь решить уравнение Пуассона в прямоугольной области с помощью быстрого преобразования Фурье-косинуса с библиотекой FFTW3.

Граничное условие четырех сторон является нулевым граничным условием Неймана.

• Уравнение d^u/dx^2+d^2u/dy^2=-f
• с f=cos(x)+cos(y), домен [-pi pi,-pi pi]
• точное решение u=f,
• функция f удовлетворяют граничному условию Неймана.
• Управляем уравнением с помощью преобразования Фурье, имеем:
• U=F/(lamda(i)+lamda(j))
• Вот lamda(k)=2*(1-cos(i/(n0-1)))
• а также U инвертируется, чтобы стать маленьким u,

Тем не менее, я все еще получаю неправильный результат по сравнению с точным решением.

Не могли бы вы мне помочь? Огромное спасибо.

Вот мой код

#include <stdio.h>
#include <math.h>
#include <fftw3.h>
#include <iostream>

using namespace std;
int main() {

int n0=64;
int n1=64;

int N=n0*n1;

double pi = 3.14159265359;
double L  = 2.*pi;
double dx = L/(n0);

double *in1=new double [N];
double *in2=new double [N];
double *out1=new double [N];
double *out2=new double [N];
double *X=new double [N];
double *Y=new double [N];fftw_plan p, q;
int i,j;
p = fftw_plan_r2r_2d(n0,n1, in1, out1, FFTW_REDFT00, FFTW_REDFT00, FFTW_ESTIMATE);
q = fftw_plan_r2r_2d(n0,n1, in2, out2, FFTW_REDFT00, FFTW_REDFT00, FFTW_ESTIMATE);

for(i = 0;i <n0;i++){
X[i] =-pi+i*dx ;
for(j = 0;j<n1;j++){
Y[j] = -pi+j*dx ;
in1[i*n0 + j]= cos(X[i]) + cos(Y[j]) ; // row major ordering
}
}

fftw_execute(p);

double *lamda=new double [n0];

for (i=0;i<n0;i++){
lamda[i]=cos(pi*i/(n0-1));
}

for ( i= 0; i< n0; i++){   // f = g / ( lamda(ii)+lamda(jj) )
for( j = 0; j < n1; j++){

double fact=0;
in2[i*n0 + j]=0;

fact=2*(2-lamda[i]-lamda[j]);

if(fact!=0){
in2[i*n0 + j] = out1[i*n0 + j]/fact;
}else{
in2[i*n0 + j] =0;
}

}
}

fftw_execute(q);

double erl1 = 0.;
for ( i = 0; i < n0; i++) {
for( j = 0; j < n1; j++){

cout<< i <<" "<< j<<" "<< cos(X[i])+cos(Y[j])<<" "<<  dx*dx*out2[i*n0+j]/(2.*(double)(n0-1))/(2.0*(double)(n1-1)) <<" "<< endl;

}
}

fftw_destroy_plan(p); fftw_destroy_plan(q); fftw_cleanup();
delete [] in1;
delete [] in2;
delete [] out1;
delete [] out2;
delete [] X;
delete [] Y;
return 0;
}