数值不稳定性 FFTW <> Matlab
Numerical instability FFTW <> Matlab
我正在尝试使用伪谱方案数值解决Swift-Hohenberg方程http://en.wikipedia.org/wiki/Swift%E2%80%93Hohenberg_equation，其中线性项在傅里叶空间中隐式处理，而非线性在实空间中进行评估。时间积分采用简单的欧拉格式。我的问题是，我已经提出的Matlab代码工作完美，而c++代码，它依赖于FFTW的傅里叶变换，变得不稳定，并在几千个时间步后发散。我已经追踪到处理非线性项的方式(参见c++代码中的注释)。如果我只用实部，就会出现不稳定性。然而，由于数值舍入误差，Phi应该只有一个可以忽略不计的虚部，Matlab正在做类似的事情，保持Phi纯粹实数。Matlab代码在Octave下也运行良好。初始条件可以像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 condition
fin.read(reinterpret_cast<char*>(Buf), N*N*sizeof(double));
fin.close();
for(int i=0; i<N*N; i++) {
Phi[i][0]=Buf[i]; //initial condition
Phi[i][1]=0.0; //no imaginary part
}
fftw_execute(phiPlan); //PhiF contains FT of initial condition
for(int j=0; j<N; j++) {
for(int i=0; i<N; i++) {
double kx=(i-(i/(N-N/2)*N))*k;
double ky=(j-(j/(N-N/2)*N))*k;
double k2=kx*kx+ky*ky;
D0[j*N+i]=1.0/((1.0 - dt*(eps-1.0 + 2.0*k2 - k2*k2))); // array of prefactors
}
}
const double norm=1.0/(N*N);
for(int n=0; n<=nSteps; n++) {
if(n%100==0) {
std::cout<<"n = "<<n<<'n';
}
for(int j=0; j<N*N; j++) {
// nonlinear term Phi^3
//NPhiF[j][0]=Phi[j][0]*Phi[j][0]*Phi[j][0]; // unstable
//NPhiF[j][1]=0.0;
NPhiF[j][0]=Phi[j][0]*Phi[j][0]*Phi[j][0] - 3.0*Phi[j][0]*Phi[j][1]*Phi[j][1];
NPhiF[j][1]=-Phi[j][1]*Phi[j][1]*Phi[j][1] + 3.0*Phi[j][0]*Phi[j][0]*Phi[j][1];
}
fftw_execute(nPhiPlan); // NPhiF contains FT of Phi^3
for(int j=0; j<N*N; j++) {
PhiF[j][0]=(PhiF[j][0] - dt*NPhiF[j][0])*D0[j]; // update
PhiF[j][1]=(PhiF[j][1] - dt*NPhiF[j][1])*D0[j];
Phi[j][0]=PhiF[j][0]*norm; // FFTW does not normalize
Phi[j][1]=PhiF[j][1]*norm;
}
fftw_execute(phiInvPlan); // Phi contains the updated Phi in real space
}
for(int i=0; i<N*N; i++) {
Buf[i]=Phi[i][0]; // saving the real part of Phi
}
std::ofstream fout("Phi.dat", std::ios::trunc | std::ios::binary);
fout.write(reinterpret_cast<char*>(Buf), N*N*sizeof(double));
fout.close();
for(int i=0; i<N*N; i++) {
Buf[i]=Phi[i][1]; // saving the imag part of Phi
}
fout.open("PhiImag.dat", std::ios::trunc | std::ios::binary);
fout.write(reinterpret_cast<char*>(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;
}
#include <iostream>
#include <fstream>
#include <cmath>
#include <fftw3.h>
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<char*>(Buf), N*N*sizeof(double));
fin.close();
for(int j=0; j<N; j++) {
for(int i=0; i<N; i++) {
((double*)PhiF)[j*2*w+i]=Buf[j*N+i];
((double*)Phi)[j*2*w+i]=Buf[j*N+i];
}
}
fftw_execute(phiPlan); //PhiF contains FT of IC
for(int j=0; j<N; j++) {
for(int i=0; i<w; i++) {
double kx=(i-(i/(N-N/2)*N))*k;
double ky=(j-(j/(N-N/2)*N))*k;
double k2=kx*kx+ky*ky;
D0[j*w+i]=1.0/(1.0 - dt*(eps-1.0 + 2.0*k2 - k2*k2));
}
}
const double norm=1.0/(N*N);
//begin first Euler step
for(int n=0; n<=nSteps; n++) {
if(n%100==0) {
std::cout<<"n = "<<n<<'n';
}
for(int j=0; j<N; j++) {
for(int i=0; i<N; i++) {
((double*)NPhi)[j*2*w+i]=((double*)Phi)[j*2*w+i] *((double*)Phi)[j*2*w+i] * ((double*)Phi)[j*2*w+i];
}
}
fftw_execute(nPhiPlan); // NPhi contains FT of Phi^3
for(int j=0; j<N*w; j++) {
PhiF[j][0]=(PhiF[j][0] - dt*NPhi[j][0])*D0[j];
PhiF[j][1]=(PhiF[j][1] - dt*NPhi[j][1])*D0[j];
}
for(int j=0; j<N*w; j++) {
Phi[j][0]=PhiF[j][0]*norm;
Phi[j][1]=PhiF[j][1]*norm;
}
fftw_execute(phiInvPlan);
}
for(int j=0; j<N; j++) {
for(int i=0; i<N; i++) {
Buf[j*N+i]=((double*)Phi)[j*2*w+i];
}
}
std::ofstream fout("Phi.dat", std::ios::trunc | std::ios::binary);
fout.write(reinterpret_cast<char*>(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 = %in', n);
end
PhiF=(PhiF - dt*NPhiF).*D0; % update
Phi=ifft2(PhiF); % inverse transform
end
return
看这些行:
<>之前为…双kx = (i - (i/(N N/2) * N)) * k;肯塔基州双= (j - j/(N N/2) * (N)) * k;双k2 = kx * kx +肯塔基州*肯塔基州;…之前我必须承认我没有研究算法，但"I/(N-N/2)"由整数组成，我怀疑你的kx, ky和k2是按预期计算的。您可以尝试这样做，以避免潜在的整数舍入错误:
<>之前为…双kx =(双(i) -(双(i)/(0.5 *双(N * N))) * k;//在我们的例子中:(N-N/2)*N) = 0.5*N*N……EDIT下面是不正确的OP说对了
指向实部的指针存储在[0]中，虚部存储在[1]中(即NPhi[1][j]是你应该参考的——至少根据他们的网站)。这可能是一个问题。固定行(我认为)应该如下:
NPhiF[0][j]=Phi[0][j]*Phi[0][j]*Phi[0][j] - 3.0*Phi[0][j]*Phi[1][j]*Phi[1][j];
NPhiF[1][j]=-Phi[1][j]*Phi[1][j]*Phi[1][j] + 3.0*Phi[0][j]*Phi[0][j]*Phi[1][j];
这将解决一部分-你将不得不解决其余的。
相关文章：
- 在C++STL中是否有Polyval(Matlab函数)等价物？
- 有可能在Armadillo中复制MATLAB circshift方法吗
- 使用 MATLAB 编码器生成C++代码：编译错误"undefined reference to `rgb2gray_tbb_real64'"
- EASTL矢量<向量<int>>连续的
- MATLAB to C++: csvread() not supported by MATLAB Coder
- 加载由 MATLAB Coder 生成的带有函数的 DLL，该函数调用外部函数
- MATLAB:跟踪imufilter对象中的状态变化
- 如何将 c++ 中的客户端 TCP 的替身列表发送到 Matlab 中的 TCP 服务器?
- 将三角函数的正确值与matlab进行比较
- 连接 MATLAB 和 Visual Studios 的问题
- 如何将C++连接到 Matlab
- 如何发送 Mat H=findHomography 返回 Matlab
- 如何在C++中编写 MATLAB fix(X） 函数?
- 提供变量作为 MATLAB 系统命令的输入参数，以便C++可执行文件
- 有没有更好的方法可以使用特征/C++实现 matlab 的逻辑索引？
- 当我运行MEX文件时，MATLAB崩溃
- 从 MATLAB 到 C++：相当于带有选项 'remove' 的 bwmorph
- MATLAB的imfinfo in openCV
- 将C++中的多个参数传递给MatLab共享库函数
- 将 matlab 数组(MAT-文件）转换为C++数组
