Hello Julia community,
I am currently going through a toy-problem trying to solve the heat PDE using the FFT (I’m following this series by Steve Brunton Solving PDEs with the FFT [Matlab] - YouTube). This is the code I came up with based on the video
using FFTW, DifferentialEquations
##
#variables
a=1 #thermal diffusivity constant
L=100 #length of domain
N=1000 #number of discrete points
dx=L/N
x=collect(-L/2:dx:L/2) #x domain
##
#define discrete wavenumbers
kappa=2π/L*collect(-N/2:N/2)
kappa=fftshift(kappa)
#parameters tuple
parameters=(kappa,a)
#initial conditions
u0=zeros(Float64,length(x))
u0[Int((L/2-L/10)/dx):Int((L/2+L/10)/dx)].=1
##
#define ODEs problem
t=(0.0,20.0)
function heatfft!(du,uhat,params,t)
du=-params[2]^2*(params[1].^2).*uhat;
end
problem=ODEProblem(heatfft!,fft(u0),t,parameters)
##
#solve
solution=solve(problem)
##inverse transform
timeee=0.0:0.1:20.0
uhat=real.(ifft.(solution(timeee)));
While the solution returns “retcode: success” and gives 8 points of the solutions that look allright, when I call the interpolation solution(timeee)
the resulting array of arrays is full of NaNs.
Am I missing something or did something incorrect ?
My Julia is at version 1.5.3 and all the packages used should be updated.