DifferentialEquations.jl: ODE returns "retcode: success" but interpolation gives "NaN"

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.

You used the in-place form but created a new array instead of mutating.

function heatfft!(du,uhat,params,t)
    du .= -params[2]^2*(params[1].^2).*uhat;
end
2 Likes

Thank you Chris for the speedy reply and your awesome work!
For posterity here is the documentation link for in-place vs out-of-place Problem Interface · DifferentialEquations.jl .

1 Like