# Spectral method for 1D heat equation OrdinaryDiffEq.jl

**URL:** <https://discourse.julialang.org/t/spectral-method-for-1d-heat-equation-ordinarydiffeq-jl/97586>\
**Category:** General Usage\
**Tags:** question, fftw, diffeq\
**Created:** [April 17, 2023, 7:25pm UTC](https://discourse.julialang.org/t/spectral-method-for-1d-heat-equation-ordinarydiffeq-jl/97586 "2023-04-17T19:25:29Z")\
**Posts on this page:** 1\
**Page:** 1

<div class="post-metadata">

**Author:** ![mleprovost](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mleprovost/32/7166_2.png) [@mleprovost](https://discourse.julialang.org/u/mleprovost)\
**Post date:** [April 17, 2023, 7:25pm UTC](https://discourse.julialang.org/t/spectral-method-for-1d-heat-equation-ordinarydiffeq-jl/97586/1 "2023-04-17T19:25:29Z")

</div>

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](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.

I would appreciate your help.

```julia
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))

```