How can I use fast Fourier transform (FFT) to solve a PDE (heat equation)?

I'm trying to solve a one-dimensional heat equation (PDE) with the Fourier transform numerically, in the way it was done here. The equation:

is subject to the initial condition: where U(x,t) is temperature, x is space, a is heat conductivity, and t is time.

I want to solve this equation using fast Fourier transform (FFT).

The key observation here is concerning the derivatives:

where k=2 pi/L[-N/2,N/2] is a spatial frequency or wave number. So, u(k,t) is a vector of Fourier coefficients and k square is a vector of frequency, so that gives n decoupled ODEs, one for each these kj.

The respective code is:

n = 1000;(*Number of discretization points*)
L = 100; (*Length of domain**)
T = 20; (*Time Integration*)
a = 1; (*thermal diffusivity constant*)
k = (2 Pi)/L Table[i, {i, -n/2, n/2 - 1}];(*define discrete wave number*)
kt = Fourier[k];(*Re-order FFT wave number*)
ic1 = Table[If[400 < i < 600, 1, 0], {i, Length[k]}]; (*initial condition*)
ict = Fourier[ic1];(*Fourier Transform of initial condition*)
ic = Table[Subscript[u, i] == ict[[i]], {i,Length[k]}];(*vector of initial condition in Fourier transform domain*)
vars = Table[Subscript[u, i][t], {i, Length[k]}]; (*vector of variables*)
eqns = Table[{Subscript[u, i]'[t] == -a (kt[[i]])^2 Subscript[u, i][t]}, {i,Length[k]}];(*model of ODEs system*)
eqn = Join[eqns, ic];
sol = NDSolve[eqn,vars, {t, 0, T}];(*Simulate in Fourier frequency domain*)

This error message appears:

General::munfl: (-1.4792510^-23+4.7229210^-303> I)+(1.5440710^-23-4.7229110^-303 I) is too small to represent as a> normalized machine number; precision may be lost. General::stop:

Further output of General::munfl will be suppressed during this > calculation.

How can I fix this?

Next, I intend to plot the solution by applying the inverse Fourier transform (IFFT).

s = Table[  Table[Evaluate[Subscript[u, j][t] /. sol], {t, 0, T}], {j,Length[k]}];(*select the solution of temperature in frequency domain*)
fin = Table[{i, j, sin[[j, i, 1]]}, {i, Length[sin[]]}, {j,Length[sin]}];
Table[ListPlot3D[fin[[i]]], {i, Length[fin]}];

The plot of the solution I'm looking for should be something like: • Why is the underflow an error? (It just means a result was so small it was converted to zero, which sometimes is perfectly fine.) What line of code causes the warning message, if it is important? Nov 4 '21 at 16:27
• @MichaelE2 I believe this message is an error because the simulation takes a long time.
– SAC
Nov 4 '21 at 18:03
• Wow what a beautiful educational video! Nov 4 '21 at 19:12

There are three mistakes here, mainly related to FFT.

First, you misunderstood the meaning of fftshift in the MATLAB code. It's not Fourier, but a shift. This has been discussed detailedly in the following post:

What's the correct way to shift zero frequency to the center of a Fourier Transform?

So we need to modify

kt = Fourier[k];

to

kt = fftshift@k;

The next mistake lies in inverse FFT. You've typed

sin = InverseFourier[s];

which leads to a 2D (to be precise, 3D, because you haven't yet stripped out the redundant {} in this step) inverse FFT, but we actually only need a bunch of 1D inverse FFT in $$k$$ direction! So the line should be modified to

sin = InverseFourier /@ Transpose@s;

Finally, visualization. ListPlot3D isn't the correct tool for the expected plot. You need the new-in-12.3 ListLinePlot3D:

ListLinePlot3D@Transpose@fin I would not use a numerical method if the problem can easily be solved analytically:

sol[u_, t_] =
u[x, t] /. DSolve[{D[u[x, t], t] == D[u[x, t], {x, 2}],
u[x, 0] == If[Abs[x] < 10, 1, 0]},
u, {x, -10, 10}, {t, 0, 20}][] The result can be plotted:

Plot3D[sol[x, t], {x, -30, 30}, {t, 0, 20}] • As I understand it, the PDE problem should be solved by FFT for didactic purposes. Nov 4 '21 at 17:45
• @user64494 Still, this answer gives an interesting alternative perspective, which on its own is plenty didactic on my book. So, +1 from me. Nov 5 '21 at 9:40