Solving the 1D heat equation using FFTW in C++

Question

I recently just picked up programming again for school and I am running across an issue with my code that involves the Fast Fourier Transform package FFTW.

In my code, I start with an initial function (in this case u(x,t=0) = sin(x) + sin(3*x) and will use RK4 to attempt to solve U_t of the heat equation. For anyone who has experience with the FFTW package, I am sending in an array of size N and transforming it, taking two derivatives, and then taking an inverse transform to solve for the U_xx of the heat equation.

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

When I plot the graph of y[i] after a few iterations of t, the values explode quite rapidly and alternate between positive and negative.

If anyone has any suggestions on syntax problems that may be causing these weird results, or if find errors in my mathematical methods, I would greatly appreciate any information. Thank you.

Answer 1

This is my first answer here so please avoid any silly mistakes. This question is asked almost 7 years ago but I found the problem in this code. Hope this will help someone.

1. Stability criteria is not considered resulting, the numerical solution isn't converging as mentioned by #duffymo. So I have added the length of the domain for normalizing the grid size and reduced the number of points in space as per the stability criteria.

2. Frequencies in Fourier-space is not defined properly. I added this part in the code and now this code is working properly.

I am attaching the modified code with this comment.

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

Answer 2

Do you know how to solve this PDE without using FFT? I would advise that you do so. It helps to know something about the answer before you start.

You said it was the heat/diffusion equation, and you gave an initial temperature distribution, but you said nothing about what the boundary conditions are. Since you said nothing, is it correct to assume that both end points are insulated (U_x = 0 at x = 0, x = L)? Or are they fixed in value at the ends? It will help if you know what those are - no solution without proper boundary conditions.