Hello,

I am looking at the heat equation \partial_t u(x,t) = \alpha \partial^2_{x} u(x,t) on \mathbb{R} with initial condition u(x,0) = f(x)

The analytical solution is given by u(x,t) = G \star f, where G(x,t) = \frac{1}{\sqrt {4\pi \alpha t}} \exp(-\frac{1}{4 \alpha t}x^2) is the heat kernel and \star denotes the convolution.

I would like to compute u(x,t) by convolution using a Fourier transform. I also solved the heat equation in spectral domain for comparison, see note by Steven Johnson (https://math.mit.edu/~stevenj/fft-deriv.pdf).

The current result of the convolution with FFT does not diffuse the initial condition. I would greatly appreciate your help.

```
# Set heat problem
α = 1.0
tf = 5.0
L = 50.0
N = 256
Δx = L/N
xgrid = collect(Δx*(0:1:N-1));
# Initial function
f(x) = exp(-(x-L/2)^2)
# Heat kernel and its exact Fourier transform
G(x,t) = 1/sqrt(4*π*α*t)*exp(-(x)^2/(4*α*t))
Ĝ(k,t) = exp(-α*k^2*t)
# Discrete initial condition
fd = zeros(N)
for (i, xi) in enumerate(xgrid)
fd[i] = f(xi)
end
f̂d = fft(fd);
# Discrete heat kernel at t=1
G1d = zeros(N)
for (i, xi) in enumerate(xgrid)
G1d[i] = G(xi, 1.0)
end
Ĝ1d = fft(G1d);
kgrid = fftfreq(N, Δx);
Ĝ1_true = map(k-> Ĝ(k,1.0), kgrid);
# Compute solution by convolution
û1d = Ĝ1d .* f̂d;
û1_true = Ĝ1_true .* f̂d;
# Numerical solution of the heat equation in spectral domain
params = Dict("N" => N,
"L" => L,
"uhat" => zeros(ComplexF64, N),
"duhat" => zeros(ComplexF64, N),
"plan" => plan_fft(zeros(N)),
"α" => α)
function heat_equation_spectral!(duhat, uhat, p, t)
α = params["α"]
N = params["N"]
L = params["L"]
for k=0:N-1
if k <= N÷2
duhat[k+1] = α*(-(2*π*k/L)^2)*uhat[k+1]
elseif k > N÷2
duhat[k+1] = α*(-(2*π/L*(k-N))^2)*uhat[k+1]
end
end
end
prob_spectral = ODEProblem((duhat,uhat,p,t) ->heat_equation_spectral!(duhat,uhat,params, t),
f̂d, (0.0, tf));
sol_spectral = solve(prob_spectral, Tsit5());
# Comparison of the results
plot(xgrid,real(ifft(sol_spectral(0.0))), label = L"t = 0.0", linewidth= 3)
plot!(xgrid,real(ifft(sol_spectral(1.0))), label = L"t = 1.0", linewidth= 3)
plot!(xgrid, real(ifft(û1d)), label = "Convolution with discrete heat kernel", linewidth= 3)
plot!(xgrid, real(ifft(û1_true)), label = "Convolution with heat kernel", linewidth= 3)
```