fftw3 для пуассона с граничным условием Дирихле для всех сторон расчетной области

Я узнаю трюк, который я вам предоставил Уравнение Пуассона с использованием БПФ с прямоугольной областью и код, который я предоставил в своем ответе Испытание путаницы , который был адаптирован из кода аскера @Charles_P! Пожалуйста, рассмотрите возможность добавления ссылок на эти URL в будущих вопросах!

Ответ на Испытание путаницы был посвящен случаю периодических граничных условий. Итак, вот несколько модификаций для решения случая граничных условий Дирихле.

fftw_plan_r2r_2d(N1,N2, in1.data(), out1.data(), FFTW_RODFT00, FFTW_RODFT00,FFTW_EXHAUSTIVE) соответствует дискретному синусоидальному преобразованию типа I как определено библиотекой FFTW:

Это значение хорошо описано в https://en.wikipedia.org/wiki/Discrete_sine_transform . Если размер массива FFTW N1=4 и его значения [a, b, c, d], полный массив, включая границы, равен [0, a, b, c, d, 0]. Отсюда пространственный шаг:

И частоты f_k типа I DST являются:

Инверсией DST типа I является DST типа I, за исключением масштабного коэффициента (см. http://www.fftw.org/doc/1d-Real_002dodd-DFTs-_0028DSTs_0029.html#g_t1d-Real_002dodd-DFTs-_0028DSTs_0029), Вот 4.(N1+1).(N2+1),

Наконец, контрольный пример должен быть адаптирован к случаю граничных условий Дирихле. Действительно, на коробке размера L1,L2 функция не соблюдает граничные условия Дирихле. Действительно, даже если исходный член один и тот же, решение, удовлетворяющее периодическим граничным условиям, может отличаться от решения, удовлетворяющего граничным условиям Дирихелта. Вместо этого можно проверить два исходных термина:

• Исходный термин соответствует одной частоте ДСТ.

• Исходный термин напрямую выводится из решения

Наконец, вот код, решающий двумерное уравнение Пуассона с использованием DST типа I библиотеки FFTW:

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

using namespace std;
int main() {

int N1=100;
int N2=200;

double pi = 3.141592653589793;
double L1  = 1.0;
double dx = L1/(double)(N1+1);//+ instead of -1
double L2= 5.0;
double dy=L2/(double)(N2+1);

double invL1s=1.0/(L1*L1);
double invL2s=1.0/(L2*L2);

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_RODFT00, FFTW_RODFT00, FFTW_EXHAUSTIVE);
q = fftw_plan_r2r_2d(N1,N2, in2.data(), out2.data(), FFTW_RODFT00, FFTW_RODFT00, FFTW_EXHAUSTIVE);

int l=0;

for(i = 0;i <N1;i++){
for(j = 0;j<N2;j++){
X[i] =dx+(double)i*dx ;
Y[j] =dy+ (double)j*dy ;
//in1[l]= sin(pi*X[i])*sin(pi*Y[j]) ; // row major ordering
in1[l]=2*Y[j]*(L2-Y[j])+2*X[i]*(L1-X[i]);
l=l+1;
}
}

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;

fact=pi*pi*((double)(i+1)*(i+1))*invL1s;

fact+=pi*pi*((double)(j+1)*(j+1))*invL2s;

in2[l] = out1[l]/fact;

}
}

fftw_execute(q);
l=-1;
double erl1 = 0.;
for ( i = 0; i < N1; i++) {
for( j = 0; j < N2; j++){
l=l+1;
X[i] =dx+(double)i*dx ;
Y[j] =dy+ (double)j*dy ;
//double res=0.5/pi/pi*in1[l];
double res=X[i]*(L1-X[i])*Y[j]*(L2-Y[j]);
erl1 +=pow(fabs(res-  0.25*out2[l]/((double)(N1+1))/((double)(N2+1))),2);
printf("%3d %10.5g %10.5g\n", l, res,  0.25*out2[l]/((double)(N1+1))/((double)(N2+1)));

}
}
erl1=erl1/((double)N1*N2);
cout<<"error=" <<sqrt(erl1) <<endl ;
fftw_destroy_plan(p); fftw_destroy_plan(q); fftw_cleanup();

return 0;
}

Составлено g++ main.cpp -o main -lfftw3 -Wall,

РЕДАКТИРОВАТЬ : Как рассчитать отклик на каждую частоту?

Идея на основе FFT состоит в том, чтобы представить решение в виде линейной комбинации функций:

• В случае периодических граничных условий используется БПФ. Основные функции:

• В случае граничных условий Дирихле используется DST типа I. Базовые функции, равные 0 при x=0 а также x=L1, являются:

  

• В случае граничных условий Неймана используется DCT типа I. Основные функции:

  

Их производные нулевые в x=0 а также x=L1,

Давайте вычислим вторую производную компонента f_k типа I DST:

Следовательно, компонент k ДСТ раствора соответствует компоненту k DST исходного слагаемого, масштабированный с коэффициентом