我试图用伪谱方案数值求解Swift-Hohenberg方程http://en.wikipedia.org/wiki/Swift%E2%80%93Hohenberg_equation,其中线性项在傅立叶空间中隐式处理,而在现实空间中评估非线性.简单的Euler方案用于时间积分.
我的问题是我提出的Matlab代码完美无缺,而依赖FFTW进行傅里叶变换的C++代码变得不稳定,经过几千个步骤后就会发散.我已经跟踪了非线性项的处理方式(参见C++代码中的注释).如果我只使用Phi的真实部分,就会发生不稳定.然而,由于数值舍入误差,Phi应该只有一个可忽略的虚部,并且Matlab正在做类似的事情,保持Phi纯粹真实.在Octave下,Matlab代码也运行良好.初始条件可能类似于
R=0.02*(rand(256,256)-0.5);
Matlab(小振幅波动).

TLDR;

为什么这些代码片段做了不同的事情？具体来说,我如何使C++代码的工作方式与Matlab版本相同？

编辑1:

为了完整起见,我使用FFTW提供的R2C/C2R功能添加了代码.有关详细信息,请参阅http://fftw.org/fftw3_doc/Multi_002dDimensional-DFTs-of-Real-Data.html(我希望我的数据布局正确).此代码始终显示约3100个时间步后的不稳定性.如果我将dt减少到例如0.01,则会发生10次.

使用复杂DFT的C++代码

#include <iostream>
#include <fstream>
#include <cmath>
#include <fftw3.h>

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 …