# Solving 1D heat equation using FFT

I am trying to reproduce the solution of heat 1D equation solution from https://databookuw.com/.

u_t = \alpha^2u_{xx}

where u(t,\ x) is the temperature distribution in time and space.

\mathcal{F}(u(t,\ x)) = \hat u(t, \ kappa)

u_x \xrightarrow{\mathcal{F}} i k \hat u

\hat u_t = - \alpha^2 k^2 \hat u

Here are the solutions in python and matlab.

and here is my take on it

using FFTW, DifferentialEquations

function heat()
α = 1.0
L = 100.0
N = 1000
dx = L / N
domain = range(-L/2, stop=L/2-dx, length=N)

u₀ = collect(0 * domain)
u₀[400:600] .= 1
u₀ = fft(u₀, 1)

kappa = (2π / L) * (-N/2:N/2-1)
kappa = fftshift(kappa, 1)

tspan = (0, 10)
t = tspan[1]:0.1:tspan[2]
params = (α, kappa)

function rhs(dû, û, p, t)
α, k = p
dû = -α^2 * (k.^2)' .* û
end
prob = ODEProblem(rhs, u₀, tspan, params)
û = solve(prob, saveat=t)
u = zeros(Complex, size(û)...)

for k in 1:length(t)
u[:, k] = ifft(û[:, k])
end
domain, real(u)
end


My results don’t match the other solutions. I have verified that up to calculating the kappa, I am consistent with the other solutions. So I am skeptical of the ODE solution or the ifft step.

1 Like

Compare your solution to a known analytic solution for the same \alpha and initial and boundary conditions. E.g. a simple sinusoid initial condition should decay exponentially.

My solution is wrong, I couldn’t make any sense of the output

But the python and matlab code solutions match the analytical solution

Why are you doing an ifft in time instead of space? It seems you messed up your indexing here.

1 Like

Probably, I’ll flip it and try

Switching the Indies fixed the shape of the solution indeed

# fixed ifft indices
#...
û = solve(prob)
# ...

for k in 1:length(eachcol(u))
u[:, k] = ifft(û[:, k])
end


But the solution doesn’t get time-stepped (i.e., the final the distribution is the same as the initial condition)

You’re not updating.

dû .= -α^2 * (k.^2)' .* û

1 Like

The solution is unstable do I need a specific solver to handle Complex numbers?

function rhs(dû, û, p, t)
α, k = p
dû .= α^2 * (k.^2) .* û
end
prob = ODEProblem(rhs, u₀, tspan, params)
û = solve(prob, Rodas4P(autodiff=false))


Aren’t you missing a -?

3 Likes

exactly

Thank you for y

Thank you for your precious time!

1 Like