Я пытаюсь численно решить уравнение Свифта-Хоэнберга http://en.wikipedia.org/wiki/Swift%E2%80%93Hohenberg_equation , используя псевдоспектральную схему, где линейная термины неявно обрабатываются в пространстве Фурье, а нелинейность оценивается в реальном пространстве. Для интегрирования по времени используется простая схема Эйлера.
Моя проблема в том, что код Matlab, который я придумал, работает отлично, в то время как код C ++, который полагается на FFTW для преобразований Фурье, становится нестабильным и расходится через пару тысяч временных шагов. Я проследил это до того, как трактуется нелинейный термин (см. Комментарии в коде C ++). Если я использую только настоящую часть Phi, возникает нестабильность. Однако у Phi должна быть лишь незначительная мнимая часть из-за ошибок численного округления, и Matlab делает нечто подобное, сохраняя Phi чисто реальным. Код Matlab также отлично работает под Octave. Начальное условие может быть чем-то вроде
R = 0,02 * (rand (256,256) -0,5);
в Matlab (небольшие колебания амплитуды).

TL; DR;

Почему эти фрагменты кода делают разные вещи? В частности, как я могу заставить код C ++ работать так же, как версия Matlab?

Редактировать 1:

Для полноты я добавил код, используя функции R2C / C2R, предоставляемые FFTW. Подробнее см. http://fftw.org/fftw3_doc/Multi_002dDimensional-DFTs-of-Real-Data.html (надеюсь, я правильно понял структуру данных). Этот код всегда показывает нестабильность примерно через 3100 шагов по времени.Если я уменьшу dt до, например, 0.01, это происходит в 10 раз позже.

Код C ++ с использованием сложных ДПФ

#include 
#include 
#include 
#include 

int main() {

const int N=256, nSteps=10000;
const double k=2.0*M_PI/N, dt=0.1, eps=0.25;

double *Buf=(double*)fftw_malloc(N*N*sizeof(double));
double *D0=(double*)fftw_malloc(N*N*sizeof(double));

// complex arrays
fftw_complex *Phi=(fftw_complex*)fftw_malloc(N*N*sizeof(fftw_complex));
fftw_complex *PhiF=(fftw_complex*)fftw_malloc(N*N*sizeof(fftw_complex));
fftw_complex *NPhiF=(fftw_complex*)fftw_malloc(N*N*sizeof(fftw_complex));

// plans for Fourier transforms
fftw_plan phiPlan=fftw_plan_dft_2d(N, N, Phi, PhiF, FFTW_FORWARD, FFTW_ESTIMATE);
fftw_plan nPhiPlan=fftw_plan_dft_2d(N, N, NPhiF, NPhiF, FFTW_FORWARD, FFTW_ESTIMATE);
fftw_plan phiInvPlan=fftw_plan_dft_2d(N, N, Phi, Phi, FFTW_BACKWARD, FFTW_ESTIMATE);

std::ifstream fin("R.dat", std::ios::in | std::ios::binary); // read initial condition
fin.read(reinterpret_cast(Buf), N*N*sizeof(double));
fin.close();
for(int i=0; i(Buf), N*N*sizeof(double));
fout.close();

for(int i=0; i(Buf), N*N*sizeof(double));
fout.close();


fftw_free(D0);
fftw_free(Buf);

fftw_free(Phi);
fftw_free(PhiF);
fftw_free(NPhiF);

fftw_destroy_plan(phiPlan);
fftw_destroy_plan(phiInvPlan);
fftw_destroy_plan(nPhiPlan);

return EXIT_SUCCESS;
}

Код C ++ с использованием R2C / C2R

#include 
#include 
#include 
#include 

int main() {

const int N=256, nSteps=3100;
const int w=N/2+1;
const double k=2.0*M_PI/N, dt=0.1, eps=0.25;

double *Buf=(double*)fftw_malloc(N*N*sizeof(double));
double *D0=(double*)fftw_malloc(N*w*sizeof(double));

fftw_complex *Phi=(fftw_complex*)fftw_malloc(N*w*sizeof(fftw_complex));
fftw_complex *PhiF=(fftw_complex*)fftw_malloc(N*w*sizeof(fftw_complex));
fftw_complex *NPhi=(fftw_complex*)fftw_malloc(N*w*sizeof(fftw_complex));

fftw_plan phiPlan=fftw_plan_dft_r2c_2d(N, N, (double*)PhiF, PhiF, FFTW_ESTIMATE);
fftw_plan nPhiPlan=fftw_plan_dft_r2c_2d(N, N, (double*)NPhi, NPhi, FFTW_ESTIMATE);
fftw_plan phiInvPlan=fftw_plan_dft_c2r_2d(N, N, Phi, (double*)Phi, FFTW_ESTIMATE);

std::ifstream fin("R.dat", std::ios::in | std::ios::binary);
fin.read(reinterpret_cast(Buf), N*N*sizeof(double));
fin.close();
for(int j=0; j(Buf), N*N*sizeof(double));
fout.close();

fftw_destroy_plan(phiPlan);
fftw_destroy_plan(phiInvPlan);
fftw_destroy_plan(nPhiPlan);

fftw_free(D0);
fftw_free(Buf);
fftw_free(Phi);
fftw_free(PhiF);
fftw_free(NPhi);
}

Код Matlab

function Phi=SwiHoEuler(Phi, nSteps)
epsi=0.25;
dt=0.1;

[nR nC]=size(Phi);
if mod(nR, 2)==0
    kR=[0:nR/2-1 -nR/2:-1]*2*pi/nR;
else
    kR=[0:nR/2 -floor(nR/2):-1]*2*pi/nR;
end
Ky=repmat(kR.', 1, nC);

if mod(nC, 2)==0
    kC=[0:nC/2-1 -nC/2:-1]*2*pi/nC;
else
    kC=[0:nC/2 -floor(nC/2):-1]*2*pi/nC;
end
Kx=repmat(kC, nR, 1); % frequencies
K2=Kx.^2+Ky.^2; % used for Laplacian in Fourier space
D0=1.0./(1.0-dt*(epsi-1.0+2.0*K2-K2.*K2)); % linear factors combined

PhiF=fft2(Phi);

for n=0:nSteps
    NPhiF=fft2(Phi.^3); % nonlinear term, evaluated in real space
    if mod(n, 100)==0
        fprintf('n = %i\n', n);
    end
    PhiF=(PhiF - dt*NPhiF).*D0; % update

    Phi=ifft2(PhiF); % inverse transform
end
return