Я боролся с очень странной ошибкой почти месяц.Прошу вас, ребята, моя последняя надежда.Я написал программу на C, которая объединяет уравнение 2d Кан-Хиллиарда с использованием схемы неявного Эйлера (IE) в пространстве Фурье (или наоборот):
Где «шляпы» означают, что мы находимся в пространстве Фурье: h_q (t_n + 1) и h_q (t_n) - это FT h (x, y) в моменты времени t_n и t_ (n + 1)), N [h_q] - нелинейный оператор, применяемый к h_q в пространстве Фурье, а L_q - линейный оператор, опять же в пространстве Фурье.Я не хочу вдаваться в подробности используемого мной численного метода, поскольку уверен, что проблема не в этом (я пытался использовать другие схемы).
Мой код на самом деледостаточно просто.Вот начало, где я в основном объявляю переменные, выделяю память и создаю планы для подпрограмм FFTW.
# include <stdlib.h>
# include <stdio.h>
# include <time.h>
# include <math.h>
# include <fftw3.h>
# define pi M_PI
int main(){
// define lattice size and spacing
int Nx = 150; // n of points on x
int Ny = 150; // n of points on y
double dx = 0.5; // bin size on x and y
// define simulation time and time step
long int Nt = 1000; // n of time steps
double dt = 0.5; // time step size
// number of frames to plot (at denominator)
long int nframes = Nt/100;
// define the noise
double rn, drift = 0.05; // punctual drift of h(x)
srand(666); // seed the RNG
// other variables
int i, j, nt; // variables for space and time loops
// declare FFTW3 routine
fftw_plan FT_h_hft; // routine to perform fourier transform
fftw_plan FT_Nonl_Nonlft;
fftw_plan IFT_hft_h; // routine to perform inverse fourier transform
// declare and allocate memory for real variables
double *Linft = fftw_alloc_real(Nx*Ny);
double *Q2 = fftw_alloc_real(Nx*Ny);
double *qx = fftw_alloc_real(Nx);
double *qy = fftw_alloc_real(Ny);
// declare and allocate memory for complex variables
fftw_complex *dh = fftw_alloc_complex(Nx*Ny);
fftw_complex *dhft = fftw_alloc_complex(Nx*Ny);
fftw_complex *Nonl = fftw_alloc_complex(Nx*Ny);
fftw_complex *Nonlft = fftw_alloc_complex(Nx*Ny);
// create the FFTW plans
FT_h_hft = fftw_plan_dft_2d ( Nx, Ny, dh, dhft, FFTW_FORWARD, FFTW_ESTIMATE );
FT_Nonl_Nonlft = fftw_plan_dft_2d ( Nx, Ny, Nonl, Nonlft, FFTW_FORWARD, FFTW_ESTIMATE );
IFT_hft_h = fftw_plan_dft_2d ( Nx, Ny, dhft, dh, FFTW_BACKWARD, FFTW_ESTIMATE );
// open file to store the data
char acstr[160];
FILE *fp;
sprintf(acstr, "CH2d_IE_dt%.2f_dx%.3f_Nt%ld_Nx%d_Ny%d_#f%.ld.dat",dt,dx,Nt,Nx,Ny,Nt/nframes);
После этой преамбулы я инициализирую свою функцию h (x, y) равномерным случайным шумоми я также принимаю FT этого.Я установил мнимую часть h (x, y), которая в коде равна dh[i*Ny+j][1], равной 0, поскольку это реальная функция.Затем я вычисляю волновые векторы qx и qy, и с ними я вычисляю линейный оператор моего уравнения в пространстве Фурье, который в коде равен Linft.Я рассматриваю только четвертую производную от h как линейный член, так что FT линейного члена просто -q ^ 4 ... но опять же, я не хочу вдаваться в детали моего метода интегрирования.Вопрос не в этом.
// generate h(x,y) at initial time
for ( i = 0; i < Nx; i++ ) {
for ( j = 0; j < Ny; j++ ) {
rn = (double) rand()/RAND_MAX; // extract a random number between 0 and 1
dh[i*Ny+j][0] = drift-2.0*drift*rn; // shift of +-drift
dh[i*Ny+j][1] = 0.0;
}
}
// execute plan for the first time
fftw_execute (FT_h_hft);
// calculate wavenumbers
for (i = 0; i < Nx; i++) { qx[i] = 2.0*i*pi/(Nx*dx); }
for (i = 0; i < Ny; i++) { qy[i] = 2.0*i*pi/(Ny*dx); }
for (i = 1; i < Nx/2; i++) { qx[Nx-i] = -qx[i]; }
for (i = 1; i < Ny/2; i++) { qy[Ny-i] = -qy[i]; }
// calculate the FT of the linear operator
for ( i = 0; i < Nx; i++ ) {
for ( j = 0; j < Ny; j++ ) {
Q2[i*Ny+j] = qx[i]*qx[i] + qy[j]*qy[j];
Linft[i*Ny+j] = -Q2[i*Ny+j]*Q2[i*Ny+j];
}
}
Тогда, наконец, наступает петля времени.По сути, я делаю следующее:
Время от времени я сохраняю данные в файл и распечатываю некоторую информацию на терминале.В частности, я печатаю самое высокое значение FT нелинейного термина.Я также проверяю, расходится ли h (x, y) на бесконечность (это не должно происходить!),
Вычисляет h ^ 3 в прямом пространстве (то есть просто dh[i*Ny+j][0]*dh[i*Ny+j][0]*dh[i*Ny+j][0]),Опять же, мнимая часть установлена на 0,
Возьмите FT h ^ 3,
Получите полный нелинейный член в обратном пространстве(то есть N [h_q] в алгоритме IE, написанном выше) путем вычисления -q ^ 2 * (FT [h ^ 3] - FT [h]).В коде я имею в виду строки Nonlft[i*Ny+j][0] = -Q2[i*Ny+j]*(Nonlft[i*Ny+j][0] -dhft[i*Ny+j][0]) и одну ниже для мнимой части.Я делаю это потому, что:
- Ускорение во времени с использованием метода IE, преобразование обратно в прямое пространство изатем нормализуйте.
Вот код:
for(nt = 0; nt < Nt; nt++) {
if((nt % nframes)== 0) {
printf("%.0f %%\n",((double)nt/(double)Nt)*100);
printf("Nonlft %.15f \n",Nonlft[(Nx/2)*(Ny/2)][0]);
// write data to file
fp = fopen(acstr,"a");
for ( i = 0; i < Nx; i++ ) {
for ( j = 0; j < Ny; j++ ) {
fprintf(fp, "%4d %4d %.6f\n", i, j, dh[i*Ny+j][0]);
}
}
fclose(fp);
}
// check if h is going to infinity
if (isnan(dh[1][0])!=0) {
printf("crashed!\n");
return 0;
}
// calculate nonlinear term h^3 in direct space
for ( i = 0; i < Nx; i++ ) {
for ( j = 0; j < Ny; j++ ) {
Nonl[i*Ny+j][0] = dh[i*Ny+j][0]*dh[i*Ny+j][0]*dh[i*Ny+j][0];
Nonl[i*Ny+j][1] = 0.0;
}
}
// Fourier transform of nonlinear term
fftw_execute (FT_Nonl_Nonlft);
// second derivative in Fourier space is just multiplication by -q^2
for ( i = 0; i < Nx; i++ ) {
for ( j = 0; j < Ny; j++ ) {
Nonlft[i*Ny+j][0] = -Q2[i*Ny+j]*(Nonlft[i*Ny+j][0] -dhft[i*Ny+j][0]);
Nonlft[i*Ny+j][1] = -Q2[i*Ny+j]*(Nonlft[i*Ny+j][1] -dhft[i*Ny+j][1]);
}
}
// Implicit Euler scheme in Fourier space
for ( i = 0; i < Nx; i++ ) {
for ( j = 0; j < Ny; j++ ) {
dhft[i*Ny+j][0] = (dhft[i*Ny+j][0] + dt*Nonlft[i*Ny+j][0])/(1.0 - dt*Linft[i*Ny+j]);
dhft[i*Ny+j][1] = (dhft[i*Ny+j][1] + dt*Nonlft[i*Ny+j][1])/(1.0 - dt*Linft[i*Ny+j]);
}
}
// transform h back in direct space
fftw_execute (IFT_hft_h);
// normalize
for ( i = 0; i < Nx; i++ ) {
for ( j = 0; j < Ny; j++ ) {
dh[i*Ny+j][0] = dh[i*Ny+j][0] / (double) (Nx*Ny);
dh[i*Ny+j][1] = dh[i*Ny+j][1] / (double) (Nx*Ny);
}
}
}
Последняя часть кода: очистите память и уничтожьте планы FFTW.
// terminate the FFTW3 plan and free memory
fftw_destroy_plan (FT_h_hft);
fftw_destroy_plan (FT_Nonl_Nonlft);
fftw_destroy_plan (IFT_hft_h);
fftw_cleanup();
fftw_free(dh);
fftw_free(Nonl);
fftw_free(qx);
fftw_free(qy);
fftw_free(Q2);
fftw_free(Linft);
fftw_free(dhft);
fftw_free(Nonlft);
return 0;
}
ЕслиЯ запускаю этот код и получаю следующий вывод:
0 %
Nonlft 0.0000000000000000000
1 %
Nonlft -0.0000000000001353512
2 %
Nonlft -0.0000000000000115539
3 %
Nonlft 0.0000000001376379599
...
69 %
Nonlft -12.1987455309071730625
70 %
Nonlft -70.1631962517720353389
71 %
Nonlft -252.4941743351609204637
72 %
Nonlft 347.5067875825179726235
73 %
Nonlft 109.3351142318568633982
74 %
Nonlft 39933.1054502610786585137
crashed!
Код падает до достижения конца, и мы видим, что нелинейный термин расходится.
Теперь, вещь, которая недля меня не имеет смысла то, что если я изменю строки, в которых я вычисляю FT нелинейного термина следующим образом:
// calculate nonlinear term h^3 -h in direct space
for ( i = 0; i < Nx; i++ ) {
for ( j = 0; j < Ny; j++ ) {
Nonl[i*Ny+j][0] = dh[i*Ny+j][0]*dh[i*Ny+j][0]*dh[i*Ny+j][0] -dh[i*Ny+j][0];
Nonl[i*Ny+j][1] = 0.0;
}
}
// Fourier transform of nonlinear term
fftw_execute (FT_Nonl_Nonlft);
// second derivative in Fourier space is just multiplication by -q^2
for ( i = 0; i < Nx; i++ ) {
for ( j = 0; j < Ny; j++ ) {
Nonlft[i*Ny+j][0] = -Q2[i*Ny+j]* Nonlft[i*Ny+j][0];
Nonlft[i*Ny+j][1] = -Q2[i*Ny+j]* Nonlft[i*Ny+j][1];
}
}
Что означает, что я использую это определение:
вместо этого:
Тогда код является абсолютно стабильным идивергенции не происходит!Даже за миллиарды временных шагов! Почему это происходит, поскольку два способа вычисления Nonlft должны быть эквивалентны?
Большое спасибо всем, кто найдет время, чтобы прочитать все это и дать мне немногоhelp!
РЕДАКТИРОВАТЬ: Чтобы сделать вещи еще более странными, я должен отметить, что эта ошибка не происходит для той же системы в 1D.В 1D оба метода вычисления Nonlft стабильны.
РЕДАКТИРОВАТЬ: я добавляю короткую анимацию того, что происходит с функцией h (x, y) непосредственно перед сбоем.Также: я быстро переписал код в MATLAB, который использует функции быстрого преобразования Фурье на основе библиотеки FFTW, и ошибка НЕ происходит ... загадка углубляется.