数值不稳定性FFTW. Matlab的 [英] Numerical instability FFTW <> Matlab
问题描述
我正在尝试数值求解Swift-Hohenberg方程 http://en .wikipedia.org/wiki/Swift％E2％80％93Hohenberg_equation 使用伪光谱方案，其中线性项在傅里叶空间中被隐式处理，而非线性则在真实空间中进行评估.一个简单的Euler方案用于时间积分.
我的问题是，我提出的Matlab代码可以完美地工作，而依赖FFTW进行傅里叶变换的C ++代码在经过几千个时间步长后变得不稳定并发散.我已经找到了处理非线性项的方式(请参见C ++代码中的注释).如果仅使用Phi的实部，则会发生不稳定.但是，由于数值舍入误差，Phi的虚部只能忽略不计，而Matlab也在做类似的事情，使Phi保持纯真.
Matlab代码在Octave下也可以正常运行.初始条件可能类似于
R=0.02*(rand(256,256)-0.5);
在Matlab中(幅度波动较小).
I am trying to numerically solve the Swift-Hohenberg equation http://en.wikipedia.org/wiki/Swift%E2%80%93Hohenberg_equation using a pseudo-spectral scheme, where the linear terms are treated implicitly in Fourier space, while the nonlinearity is evaluated in real space. A simple Euler scheme is used for the time integration.
My problem is that the Matlab code I have come up with works perfectly, while the C++ code, which relies on FFTW for the Fourier transforms, becomes unstable and diverges after a couple thousand time steps. I have tracked it down to the way the nonlinear term is treated (see the comments in the C++ code). If I use only the real part of P the instability occurs. However, Phi should only have a negligible imaginary part due to numerical rounding errors, and Matlab is doing something similar, keeping Phi purely real.
The Matlab code also runs fine under Octave. The initial condition can be something like
R=0.02*(rand(256,256)-0.5);
in Matlab (small amplitude fluctuations).
为什么这些代码段会做不同的事情?具体来说，如何使C ++代码像Matlab版本一样工作?
Why do these pieces of code do different things? Specifically, how can I make the C++ code work the same way the Matlab version does?
为了完整起见，我已经使用FFTW提供的R2C/C2R功能添加了代码.参见 http://fftw.org/fftw3_doc/Multi_002dDimensional-DFTs-of-有关详细信息，请参见Real-Data.html (我希望我的数据布局正确).该代码始终显示经过约3100个时间步长后的不稳定性.如果我将dt减少到0.01，它会在10次后发生.
For completeness, I have added the code using the R2C/C2R functions provided by FFTW. See http://fftw.org/fftw3_doc/Multi_002dDimensional-DFTs-of-Real-Data.html for details (I hope I got the data layout right). This code always shows the instability after about 3100 time steps. If I reduce dt to e.g. 0.01, it occurs 10 times later.
#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;
}
使用R2C/C2R的C ++代码
#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 = %i\n', n);
end
PhiF=(PhiF - dt*NPhiF).*D0; % update
Phi=ifft2(PhiF); % inverse transform
end
return
推荐答案
看以下几行:
for ...
double kx=(i-(i/(N-N/2)*N))*k;
double ky=(j-(j/(N-N/2)*N))*k;
double k2=kx*kx+ky*ky;
...
我必须承认我没有研究算法，但是"i/(N-N/2)"由整数组成，我怀疑您的kx，ky和k2是按预期计算的.您可以尝试执行以下操作以避免潜在的整数舍入错误:
I have to admit that I did not look into the algos but "i/(N-N/2)" consists of integers and I suspect that your kx, ky, and k2 are calculated as expected. You may try something like this which avoids potential integer rounding errors:
for ...
double kx=( double(i) -( double(i)/(0.5*double(N*N)))*k; // where in our case: (N-N/2)*N) = 0.5*N*N
...
...
这篇关于数值不稳定性FFTW. Matlab的的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
