在 C++ 中使用 FFTW 求解一维热方程

Question

我最近刚刚在学校再次学习编程，我的代码遇到了一个涉及快速傅里叶变换包 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;
}

当我绘制 y[i 的图时 。 ] 在 t 的几次迭代之后，值会迅速爆炸并在正负之间交替

如果有人对可能导致这些奇怪结果的语法问题有任何建议，或者如果在我的数学方法中发现错误，我将不胜感激。谢谢。

Answer 1

你知道如何在不使用 FFT 的情况下求解这个偏微分方程吗？我建议你这样做。在开始之前了解一些答案会很有帮助。

您说这是热/扩散方程，并且给出了初始温度分布，但您没有提到边界条件是什么。既然您什么也没说，那么假设两个端点都是绝缘的（U_x = 0，x = 0，x = L）是否正确？或者它们的价值最终是固定的？如果您知道这些是什么，将会有所帮助 - 没有适当的边界条件就没有解决方案。

Answer 2

这是我在这里的第一个答案，所以请避免任何愚蠢的错误。这个问题大约 7 年前就被问到了，但我在这段代码中发现了问题。希望这会对某人有所帮助。

1. 稳定性标准不被认为是结果，数值解没有像 #duffymo 提到的那样收敛。因此，我添加了域的长度以标准化网格大小，并根据稳定性标准减少了空间中的点数。

2. 傅立叶空间中的频率未正确定义。我在代码中添加了这部分，现在该代码可以正常工作。

我将修改后的代码附加到此评论中。

#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;
}