Spectral method for 1D heat equation OrdinaryDiffEq.jl

General Usage
, ,
1

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.

I would appreciate your help.

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