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

Question

Я пытаюсь решить уравнение Яда с граничным условием Дирихле для четырех сторон вычислительной области. Как известно, я должен использовать FFTW_RODFT00 для выполнения условия. Однако результат неверен. Не могли бы вы мне помочь?

#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=100; 

double pi = 3.141592653589793; 
double L1 = 2.0; 
double dx = L1/(double)(N1-1); 
double L2= 2.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=-1; 
for(i = 0;i <N1;i++){ 
    X[i] =-1.0+(double)i*dx ;   
    for(j = 0;j<N2;j++){ 
     l=l+1; 
     Y[j] =-1.0+ (double)j*dy ; 
     in1[l]= sin(pi*X[i]) + sin(pi*Y[j]) ; // row major ordering 
    } 
} 

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; 
     in2[l]=0; 

     if(2*i<N1){ 
      fact=((double)i*i)*invL1s;; 
     }else{ 
      fact=((double)(N1-i)*(N1-i))*invL1s; 
     } 
     if(2*j<N2){ 
      fact+=((double)j*j)*invL2s; 
     }else{ 
      fact+=((double)(N2-j)*(N2-j))*invL2s; 
     } 
     if(fact!=0){ 
      in2[l] = out1[l]/fact; 
     }else{ 
      in2[l] = 0.0; 
     } 
    } 
} 

fftw_execute(q); 
l=-1; 
double erl1 = 0.; 
for (i = 0; i < N1; i++) { 
    for(j = 0; j < N2; j++){ 
     l=l+1; 

     erl1 +=1.0/pi/pi*fabs(in1[l]- 0.25*out2[l]/((double)(N1-1))/((double)(N2-1))); 
     printf("%3d %10.5f %10.5f\n", l, in1[l], 0.25*out2[l]/((double)(N1-1))/((double)(N2-1))); 

    } 
} 

cout<<"error=" <<erl1 <<endl ; 
fftw_destroy_plan(p); fftw_destroy_plan(q); fftw_cleanup(); 

return 0; 

}

Answer 1

Я признаю трюк я предоставил вам в Poisson equation using FFTW with rectanguar domain и код я предоставил в своем ответе на Confusion testing fftw3 - poisson equation 2d test, который был адаптирован из кода Аскер @Charles_P! Пожалуйста, подумайте о добавлении ссылок на этот URL в будущих вопросах!

Ответ на вопрос Confusion testing fftw3 - poisson equation 2d test был посвящен случаю периодических граничных условий. Итак, вот несколько модификаций для решения случая граничных условий Дирихле.

fftw_plan_r2r_2d(N1,N2, in1.data(), out1.data(), FFTW_RODFT00, FFTW_RODFT00,FFTW_EXHAUSTIVE) соответствует типу I дискретного преобразования синус as defined by the FFTW library:

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

и частоты f_k типа 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.

• Термин источник получают непосредственно из раствора

Наконец, вот код решения уравнения 2D Пуассона с использованием типа I DST библиотеки 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.

EDIT: Как вычислить реакцию на каждую частоту?

Идея БПФ основе заключается в представлении решения в виде линейной комбинации функций:

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

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

  

• В случае граничных условий Неймана, тип-я ДКП используется. Базовые функции:

  

их производные равны нулю в x=0 и x=L1.

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

Следовательно, componant k от DST раствора соответствует componant k от DST из источник, масштабированный в