Solving 1D heat equation using FFT

Abdulrhmn_Ghanem · March 11, 2024, 12:15pm

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.

John_Gibson · March 11, 2024, 12:29pm

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.

Abdulrhmn_Ghanem · March 11, 2024, 12:43pm

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

But the python and matlab code solutions match the analytical solution

ChrisRackauckas · March 11, 2024, 2:30pm

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

Abdulrhmn_Ghanem · March 11, 2024, 2:56pm

Probably, I’ll flip it and try

Abdulrhmn_Ghanem · March 11, 2024, 3:18pm

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)

ChrisRackauckas · March 11, 2024, 4:26pm

You’re not updating.

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

Abdulrhmn_Ghanem · March 11, 2024, 9:55pm

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

ChrisRackauckas · March 12, 2024, 12:02pm

Aren’t you missing a -?

photor · March 12, 2024, 12:56pm

exactly

Abdulrhmn_Ghanem · March 12, 2024, 3:32pm

Thank you for y

Thank you for your precious time!