用FFTW在c++中求解一维热方程
Solving the 1D heat equation using FFTW in C++
我最近刚刚为学校重新开始编程,我遇到了一个涉及快速傅里叶变换包FFTW的代码问题。
在我的代码中,我从一个初始函数开始(在这种情况下,u(x,t=0) = sin(x) + sin(3*x)),并将使用RK4尝试解决热方程的U_t。对于任何使用过FFTW包的人来说,我将发送一个大小为N的数组并对其进行变换,对其进行两次导数,然后进行反变换以求解热量方程的U_xx。
#include <iostream>
#include <cmath>
#include <fstream>
#include <fftw3.h>
using namespace std;
const double Pi= 3.141592653589;
const int N = 2048;
void Derivative(double& num1, double& num2, int J){
//takes two derivatives
num1 = (-1)*J*J*num1;
num2 = (-1)*J*J*num2;
}
double Function(double X){ // Initial function
double value = sin(X) + sin(3*X);
return value;
}
void FFT(double *in, double *Result){
double *input, *outI;
input = new double [N];
outI = new double [N];
for(int i=0; i<N; i++){
input[i] = in[i];
}
//Declaring transformed matricies;
fftw_complex *out;
out= (fftw_complex*) fftw_malloc(sizeof(fftw_complex)*((N/2)+1));
//Plans for execution
fftw_plan p;
p= fftw_plan_dft_r2c_1d(N,input, out, FFTW_ESTIMATE);
fftw_execute(p);
//Derivative
for(int i=0; i<(N/2+1);i++){
Derivative(out[i][0],out[i][1],i);
}
fftw_plan pI;
//Execution of Inverse FFT
pI = fftw_plan_dft_c2r_1d(N,out,outI,FFTW_PRESERVE_INPUT);
fftw_execute(pI);
for(int i=0;i<N; i++){//Dividing by the size of the array
Result[i] = (1.0/N)*outI[i];
}
fftw_free(outI); fftw_free(out);
fftw_free(input);
fftw_destroy_plan(pI); fftw_destroy_plan(p);
}
int main()
{
//Creating and delcaring function variables
double *initial,*w1,*w2,*w3,*w4,*y,*holder1,*holder2, *holder3;
double dx= 2*Pi/N, x=0, t=0;
double h = pow(dx,2), tmax= 100*h;
//creating arrays for the RK4 procedure
initial = new double [N];
w1 = new double [N]; w2 = new double [N]; y= new double [N];
w3 = new double [N]; w4 = new double [N]; holder1 = new double [N];
holder2 = new double [N]; holder3 = new double [N];
double *resultArray=new double[N];
int j =0;
//Output files
ofstream sendtofileINITIAL("HeatdataINITIAL.dat");
ofstream sendtofile5("Heatdata5.dat");
ofstream sendtofile25("Heatdata25.dat");
ofstream sendtofile85("Heatdata85.dat");
ofstream sendtofileFINAL("HeatdataFINAL.dat");
//Initial data
for(int i=0; i<N; i++){
y[i] = Function(x);
sendtofileINITIAL << i*2*Pi/N <<" "<< y[i] << endl;
x+=dx;
}
//RK4
// y[i+1] = y[i] + (h/2)*(w1 + 2*w2 + 2*w3 + w4)
// where the w1,w2,w3,w4 is the standard RK4 calculations
// w1= f(y[i]), w2= f(y[i]+(h/2)*w1[i])
// w3= f(y[i] + (h/2)*w2[i]) and w4 = f(y[i]+h*w3[i])
while(t<=tmax){
FFT(y,resultArray);
for(int i=0; i<N;i++){ //Calculating w1
w1[i] = resultArray[i];
}
for(int i=0; i<N; i++){
holder1[i] = y[i] + (h/2)*w1[i];}
FFT(holder1,resultArray);
for(int i=0; i<N; i++){ //Calculating w2
w2[i] = resultArray[i];
}
for(int i=0; i<N; i++){
holder2[i] = y[i] + (h/2)*w2[i];
}
FFT(holder2,resultArray);
for(int i=0; i<N; i++){ //Calculating w3
w3[i] = resultArray[i];
}
for(int i=0; i<N; i++){
holder3[i]= y[i] + h*w3[i];
}
FFT(holder3,resultArray);
for(int i=0; i<N; i++){ //Calculating w4
w4[i] = resultArray[i];
}
for(int i=0; i<N; i++){
y[i] += (h/6.0)*(w1[i] + 2*w2[i] + 2*w3[i] +w4[i]);
//Outputing data at certain time periods
if(j==5)
sendtofile5 << i*2*Pi/N <<" "<< y[i] << endl;
if(j==25)
sendtofile25 << i*2*Pi/N <<" "<< y[i] << endl;
if(j==85)
sendtofile85 << i*2*Pi/N <<" "<< y[i] << endl;
if(j==100)
sendtofileFINAL << i*2*Pi/N <<" "<< y[i] << endl;
}
j++;// counter
t+=h;// time counter
}
sendtofileINITIAL.close();
sendtofile5.close();
sendtofile25.close();
sendtofile85.close();
sendtofileFINAL.close();
return 0;
}
当我在t的几次迭代后绘制y[I]的图时,值迅速爆炸,并在正负之间交替。
如果有人对可能导致这些奇怪结果的语法问题有任何建议,或者如果发现我的数学方法错误,我将非常感谢任何信息。谢谢你。
你知道如何在不使用FFT的情况下解决这个PDE吗?我建议你这样做。在你开始之前对答案有所了解是有帮助的。
你说这是热/扩散方程,你给出了一个初始温度分布,但你没有说边界条件是什么。既然你什么也没说,那么假设两个端点都是绝缘的(在x = 0, x = L时U_x = 0)是否正确?或者它们最终的价值是固定的?如果你知道它们是什么就会有帮助——没有适当的边界条件就没有解。
这是我在这里的第一个答案,所以请避免任何愚蠢的错误。这个问题几乎是7年前问的,但我在这个代码中发现了问题。
-
稳定性标准不被认为是结果,数值解不收敛,如#duffymo所提到的。因此,我添加了域的长度,用于规范化网格大小,并根据稳定性标准减少了空间中的点数。
-
傅里叶空间中的频率没有被正确定义。我在代码中添加了这一部分,现在这段代码工作正常。
我将修改后的代码与此注释附加在一起。
#include <iostream>
#include <cmath>
#include <fstream>
#include <fftw3.h>
using namespace std;
const double Pi = 3.141592653589;
const int N = 64;
const double L = 5.0;
void Derivative(double& num1, double& num2, double J) {
//takes two derivatives
num1 = -1 * pow(J, 2)*num1;
num2 = -1 * pow(J, 2)*num2;
}
double Function(double X) { // Initial function
double value = 1.0 / cosh(3 * (X - M_PI));//sin(X) + sin(3*X);
return value;
}
void FFT(double *in, double *Result) {
double *input, *outI, *kx;
input = new double[N];
outI = new double[N];
kx = new double[N];
for (int i = 0; i<N; i++) {
input[i] = in[i];
}
//Declaring transformed matricies;
fftw_complex *out;
out = (fftw_complex*)fftw_malloc(sizeof(fftw_complex)*(N));
//Plans for execution
fftw_plan p;
p = fftw_plan_dft_r2c_1d(N, input, out, FFTW_ESTIMATE);
fftw_execute(p);
//----This part calculate frequencies--------//
//Frequencies
for (int i = 0; i< N / 2; ++i) {
kx[i] = i / L;
}
kx[N / 2] = 0.00;
for (int i = 0; i < ((N / 2) - 1); ++i) {
kx[i + 1 + N / 2] = -kx[N / 2 - i - 1];
}
//Derivative
for (int i = 0; i<N; i++) {
//the frequencies is passed to this function to calculate
//the derivatives
Derivative(out[i][0], out[i][1], kx[i]);
}
fftw_plan pI;
//Execution of Inverse FFT
pI = fftw_plan_dft_c2r_1d(N, out, outI, FFTW_PRESERVE_INPUT);
fftw_execute(pI);
for (int i = 0; i<N; i++) {//Dividing by the size of the array
Result[i] = (1.0 / N)*outI[i];
}
fftw_free(outI); fftw_free(out);
fftw_free(input);
fftw_destroy_plan(pI); fftw_destroy_plan(p);
}
int main()
{
//Creating and delcaring function variables
double *initial, *w1, *w2, *w3, *w4, *y, *holder1, *holder2, *holder3;
double dx = 2 * Pi / N, x = 0, t = 0, dt = 0.01;
double h = pow(dx, 2), tmax = 100 * h;
//creating arrays for the RK4 procedure
initial = new double[N];
w1 = new double[N]; w2 = new double[N]; y = new double[N];
w3 = new double[N]; w4 = new double[N]; holder1 = new double[N];
holder2 = new double[N]; holder3 = new double[N];
double *resultArray = new double[N];
double *resultnext = new double[N];
int j = 1;
// The stability parameter are declared but not used
float alpha = 0.1, CFL = alpha * dt / (10 * pow(dx, 2));
//Output files
ofstream sendtofileINITIAL("HeatdataINITIAL.dat");
ofstream sendtofile5("Heatdata5.dat");
ofstream sendtofile25("Heatdata25.dat");
ofstream sendtofile85("Heatdata85.dat");
ofstream sendtofileFINAL("HeatdataFINAL.dat");
//Initial data
for (int i = 0; i < N; i++) {
y[i] = Function(x);
sendtofileINITIAL << i * 2 * Pi / N << " " << y[i] << endl;
x += dx;
}
//RK4
// y[i+1] = y[i] + (h/2)*(w1 + 2*w2 + 2*w3 + w4)
// where the w1,w2,w3,w4 is the standard RK4 calculations
// w1= f(y[i]), w2= f(y[i]+(h/2)*w1[i])
// w3= f(y[i] + (h/2)*w2[i]) and w4 = f(y[i]+h*w3[i])
while (t <= tmax) {
FFT(y, resultArray);
for (int i = 0; i < N; i++) { //Calculating w1
w1[i] = resultArray[i];
}
for (int i = 0; i < N; i++) {
holder1[i] = y[i] + (h / 2)*w1[i];
}
FFT(holder1, resultArray);
for (int i = 0; i < N; i++) { //Calculating w2
w2[i] = resultArray[i];
}
for (int i = 0; i < N; i++) {
holder2[i] = y[i] + (h / 2)*w2[i];
}
FFT(holder2, resultArray);
for (int i = 0; i < N; i++) { //Calculating w3
w3[i] = resultArray[i];
}
for (int i = 0; i < N; i++) {
holder3[i] = y[i] + h * w3[i];
}
FFT(holder3, resultArray);
for (int i = 0; i < N; i++) { //Calculating w4
w4[i] = resultArray[i];
}
for (int i = 0; i < N; i++) {
y[i] += (h / 6.0)*(w1[i] + 2 * w2[i] + 2 * w3[i] + w4[i]);
//Outputing data at certain time periods
if (j == 5)
sendtofile5 << i * 2 * Pi / N << " " << y[i] << endl;
if (j == 25)
sendtofile25 << i * 2 * Pi / N << " " << y[i] << endl;
if (j == 85)
sendtofile85 << i * 2 * Pi / N << " " << y[i] << endl;
if (j == 100)
sendtofileFINAL << i * 2 * Pi / N << " " << y[i] << endl;
}
j++;// counter
t += h;// time counter
}
sendtofileINITIAL.close();
sendtofile5.close();
sendtofile25.close();
sendtofile85.close();
sendtofileFINAL.close();
return 0;
}
相关文章:
- 如何将三维尺寸不固定的三维阵列展平为一维阵列
- STL算法函数在多个一维容器上的使用
- 将二维数组的所有元素插入到一维数组中
- C++语法差异:二维和一维数组(指针算术)
- 将一维数组写入 CSV C++中的不同列?
- C++:将矩阵存储在一维数组中
- 当表示为对象的一维向量时,有效地旋转 NxM 矩阵 (C++)
- 如何在一维数组中的每个元素中都有多个int值
- 以C++填充一维数组
- 一维阵列的运动检测(神经网络或其他选项?
- 编写所需的代码以创建动态一维整数数组
- 用于在一维数组上嵌套循环操作的正确 openmp 指令
- 使用两个不同大小的一维阵列制作 2D 阵列
- 如何在 <threads> c++ 中使用和一维数组进行矩阵乘法?
- C++按内存地址将多维数组更改为一维数组
- 在 c++ 中返回一维数组时出错
- 新的一个一维阵列,非常大,例如60000*60000
- 如何使用嵌套初始化构造函数中的一维向量初始化矩阵
- 随时间变化的一维薛定谔方程C++
- 用FFTW在c++中求解一维热方程