Solving 1D heat equation using FFT

Specific Domains
1

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(dû, û, p, t)
		α, k = p 
		dû = -α^2 * (k.^2)' .* û
	end
	prob = ODEProblem(rhs, u₀, tspan, params)
	û = solve(prob, saveat=t)
	u = zeros(Complex, size(û)...)
	
	for k in 1:length(t)
		u[:, k] = ifft(û[:, 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.

2

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.

3

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

But the python and matlab code solutions match the analytical solution

image
image959×687 107 KB

4

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

5

Probably, I’ll flip it and try

6

Switching the Indies fixed the shape of the solution indeed

# fixed ifft indices
#...
û = solve(prob)
# ...
	
for k in 1:length(eachcol(u))
	u[:, k] = ifft(û[:, k])
end

image
image1105×657 16.4 KB

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

image
image1062×626 14.5 KB

7

You’re not updating.

dû .= -α^2 * (k.^2)' .* û

8

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

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

image
image1090×172 27.4 KB

9

Aren’t you missing a -?

10

exactly

11

Thank you for y

:man_facepalming:

Thank you for your precious time!