1D Heat Equation using DFT produces incorrect results (FFTW)

Question

I am trying to solve the 1D heat equation using a complex to complex IDFT. The problem is that the output after a single timestep does not seem to be correct. I have included a simple example below to illustrate the problem.

I initialize the temperature state as follows:

The initial modes in the frequency domain are:
k[ 0] = 12.5 + 0i
k[ 1] = 12.5 + 0i
k[ 2] = 12.5 + 0i
k[ 3] = 12.5 + 0i
k[ 4] = 12.5 + 0i
k[-3] = 12.5 + 0i
k[-2] = 12.5 + 0i
k[-1] = 12.5 + 0i

I then advance the state of the frequency domain to t=0.02 using the standard 1D heat equation:

double alpha = 0.2; // Thermal conductivity constant
double timestep = 0.02;

for (int i = 0; i < N; i++) {
    int k = (i <= N / 2) ? i : i - N;

    F[i][REAL] *= exp(-alpha * k * k * timestep); // Decay the real part
    F[i][IMAG] *= exp(-alpha * k * k * timestep); // Decay the imaginary part
}

The frequency modes at t=0.02 become:
k[ 0] = 12.5 + 0i
k[ 1] = 12.45 + 0i
k[ 2] = 12.3 + 0i
k[ 3] = 12.06 + 0i
k[ 4] = 11.73 + 0i
k[-3] = 12.06 + 0i
k[-2] = 12.3 + 0i
k[-1] = 12.45 + 0i

After performing the IDFT to obtain the temperature domain state at t=0.02 I get:

The spatial and frequency domain both seem to be correctly periodic. However, the heat (values in spatial domain) does not seem to dissipate according to a Gaussian curve. Even more surprisingly, some temperatures dip below their initial value (they become negative!).

The conservation of energy does seem to hold correctly: adding all the temperatures together still nets 100.

This is my full heat equation code:

double alpha = 0.2;     // Thermal conductivity constant
double timestep = 0.02; // Physical heat equation timestep
int N = 8;              // Number of data points

fftw_complex* T = (fftw_complex*)fftw_alloc_complex(N); // Temperature domain
fftw_complex* F = (fftw_complex*)fftw_alloc_complex(N); // Frequency domain

fftw_plan plan = fftw_plan_dft_1d(N, F, T, FFTW_BACKWARD, FFTW_MEASURE); // IDFT from frequency to temperature domain

// Initialize all frequency modes such that there is a peak of 100 at x=0 in the temperature domain
// All other other points in the temperature domain are 0
for (int i = 0; i < N; i++) {
    F[i][REAL] = 100.0 / N;
    F[i][IMAG] = 0.0;
}

// Perform the IDFT to obtain the initial state in the temperature domain
fftw_execute(plan);
printTime1d(T, N);
printFrequencies1d(F, N);

// Perform a single timestep of the heat equation to obtain the frequency domain state at t=0.02
for (int i = 0; i < N; i++) {
    int k = (i <= N / 2) ? i : i - N;

    F[i][REAL] *= exp(-alpha * k * k * timestep); // Decay the real part
    F[i][IMAG] *= exp(-alpha * k * k * timestep); // Decay the imaginary part
}

// Perform the IDFT to obtain the temperature domain state at t=0.02
fftw_execute(plan);
printTime1d(T, N);
printFrequencies1d(F, N);

The definition of printTime(...) and printFrequencies(...) is:

void printTime1d(fftw_complex* data, int N) {
    int rounding_factor = pow(10, 2);

    for (int i = 0; i < N; i++) {
        std::cout << std::setw(8) << round(data[i][REAL] * rounding_factor) / rounding_factor;
    }

    std::cout << std::endl;
}

void printFrequencies1d(fftw_complex* data, int N) {
    int rounding_factor = pow(10, 2);

    for (int i = 0; i < N; i++) {
        int k = (i <= N / 2) ? i : i - N;

        double R = round(data[i][REAL] * rounding_factor) / rounding_factor;
        double I = round(data[i][IMAG] * rounding_factor) / rounding_factor;

        std::cout << "k[" << std::setw(2) << k << "]: " << std::setw(2) << R << ((I < 0) ? " - " : " + ") << std::setw(1) << abs(I) << "i" << std::endl;
    }

    std::cout << std::endl;
}

Perhaps good to note is that I have also conducted this experiment using a complex to real IDFT (with fftw's fftw_plan_dft_c2r_1d()) and it gave the exact same results.

Answer 1

Your problem is that you don't resolve the needed frequencies, getting instead the following Fourier image of the function after multiplication by the decay coefficients:

The above result is too far from what you should get instead—a Gaussian—at least something like this (using 80 points instead of 8):

Notice how the amplitudes in the first chart above don't even have any chance to come even close to zero, instead bumping into the Nyquist frequency. It's then obvious that you'll get artifacts resembling the Gibbs phenomenon: it's the usual behavior of Fourier partial sums.

The inverse Fourier transform of the 80-point data version is then as follows:

This result still does have negative components (since we use a finite number of harmonics), but these are much smaller in amplitude than what you got with only 8 harmonics.

Note that this does mean that, if you increase the value of time you're interested in, you could reduce the number of harmonics taken into account. This might be unexpected at first, but it's simply because the upper harmonics decay much faster than the lower ones, and they never ever increase back.