# Spectral method for 1D heat equation OrdinaryDiffEq.jl

Hello,

I would like to solve the heat equation \frac{\partial u(x,t)}{\partial t} = \alpha \frac{\partial^2 u}{\partial x^2} on the periodic domain [0, 1] with initial condition u(x, t = 0) = u_0(x) and periodic boundary conditions u(0, t) = u(1, t) with a spectral method in OrdinaryDiffEq.jl

I followed the note by @stevengj https://math.mit.edu/~stevenj/fft-deriv.pdf to compute the second order partial derivative.

Is it a correct one to solve the heat equation with FFT?
I tried to verify the result by convolution of the Green’s function for the heat equation with the initial condition. I get large discrepancies between the two approaches.

using FFTW
using LinearAlgebra
using OrdinaryDiffEq

L = 1.0
N = 128
α = 0.1
Δx = L/N
xgrid = collect(Δx*(0:1:N-1))
tf = 5.0

u0 = zeros(N)
for (j, xj) in enumerate(xgrid)
u0[j] = exp(-(xj-0.5)^2/(2*α^2))
end

u0hat = fft(u0)

params = Dict("N" => N,
"L" => L,
"α" => α)

function heat_equation!(duhat, uhat, p, t)

α = params["α"]
N = params["N"]
L = params["L"]

# Compute RHS in spectral space
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 = ODEProblem((duhat,uhat,p,t) ->heat_equation!(duhat,uhat,params, t),
u0hat, (0.0, tf))

sol = solve(prob, Tsit5())

# Validate solution
heat_kernel(x,t) = 1/sqrt(4*π*α*t)*exp(-(x)^2/(4*α*t))

discrete_kernel = zeros(ComplexF64, N)

for (j, xj) in enumerate(xgrid)
discrete_kernel[j] = heat_kernel(xj, 3.0)
end

fft!(discrete_kernel)

convolution_sol = ifft(discrete_kernel .* u0hat)

norm(real(ifft(sol(3.0))), real(convolution_sol))

1 Like