# Issue with ODE solver for solving Burgers equation using spectral method

Hi all,

I am trying to solve the 1D Burgers equation Burgers equation wiki using Fourier transform and the ODE package.
Things works smoothly in Matlab, but not in Julia. I’m just wondering where I went wrong.

Matlab version:

close all; clear; clc
%% Parameters
nu = 0.02; % drag coeff.
nT = 1001; % total timesteps
dt = 0.02; % timestep
L  = 16;   %
N  = 256;  % number of grid in x
x  = L/N*(-N/2:N/2-1); % discretized grid location
dx = x(2) - x(1);      % grid interval
kx = (2*pi/L)*[0:N/2-1 -N/2:-1]';
u  = exp(-(x+3).^2);
ut = fft(u);
t  = (1:nT)*dt;

%% Main
[t,utsol] = ode45(@(t,ut) burgers_rhs(t,ut,nu,kx), t, ut);

usol = ifft(utsol,[],2);

%%
figure('unit','normalized','position',[0.1,0.1,0.4,0.5],...
'DefaultAxesFontSize',12);

for it=1:floor(nT/5):nT
if it<=1
plot(x,usol(it,:),'r:','LineWidth',2); hold on
else
plot(x,usol(it,:),'b','LineWidth',2); hold on
end
[ym,idx] = max(usol(it,:));
text(x(idx),ym+0.05,['t=',num2str(t(it))])
end

title(['Burgers FFT, \nu=',num2str(nu),', L=',num2str(L),', dx=', ...
num2str(dx),', dt=',num2str(dt),', nt=',num2str(nT)])
xlim([min(x),max(x)]); ylim([0,1.1])
xlabel('x'); ylabel('u')
print(gcf,'-dpng',['burgers_fft_nu=',num2str(nu),',L=',num2str(L),...
',dx=',num2str(dx),',dt=',num2str(dt),...
',nt=',num2str(nT),'.png'])

%%
function dut=burgers_rhs(t,ut,nu,kx)

u = ifft(ut);
dut = -nu*kx.^2.*ut - fft(u.*ifft(1i*kx.*ut));

end


The above Matlab code works fine, both with ode45 and ode23.

Julia version:

using PyPlot, FFTW, DifferentialEquations

"Solve Burgers equation with spectral method."
function solveBurgers_FT()

# Parameters
nu = 0.02 # drag coeff.
nt = 1001 # total timesteps
dt = 0.02 # timestep
L  = 16   #
N  = 256  # number of grid in x
x  = L/N*(-N/2:N/2-1) # discretized grid location
dx = x[2] - x[1]      # grid interval
kx = (2*pi/L)*[0:N/2-1... -N/2:-1...]
u  = exp.(-(x.+3).^2)
ut = fft(u)
#t  = (1:nt)*dt
tspan = (0.0, nt*dt)
#tspan = (dt, 2*dt) # for debug
p = [nu, kx]

# Main
prob = ODEProblem(burgers_rhs,ut,tspan,p)

# I need to find out why none of these works!
#utsol = solve(prob,DP5(),saveat=dt,reltol=1e-4,abstol=1e-4,dt=dt)
utsol = solve(prob,DP5(),saveat=dt,reltol=1e-4,abstol=1e-4,dt=0.0002)
#utsol = solve(prob,Tsit5(),saveat=dt)
#utsol = solve(prob,Rodas5(),saveat=dt)
#utsol = solve(prob,alg_hints=[:stiff],saveat=dt)

usol = ifft(hcat(utsol.u...),1)

# Visualization
figure()

for it=1:nt

if mod(it,floor(Int,nt/5))==1
if it<=1
plot(x,usol[:,it],"r:",LineWidth=2)
else
plot(x,usol[:,it],"b",LineWidth=2)
end
uabsmax = findmax(abs.(usol[:,it]))
umax = usol[uabsmax[2],it]
text(x[uabsmax[2]],umax+0.05,"t=$((it-1)*dt)") end end title("Burgers, \nu=$(nu), L= $(L), dx=$(dx), dt= $(dt), nt=$(nt)")
xlim([x[1],x[end]])
ylim([0,1.1])
xlabel("x")
ylabel("u")
end

function burgers_rhs(dut,ut,p,t)

nu,kx = p[1], p[2]
u = ifft(ut)
dut = -nu*kx.^2 .* ut - fft(u.*ifft(1im*kx.*ut))

end


However, when I transformed the code to Julia, it either shows stability issue or the solution does not evolve at all. What am I doing wrong?

as a starter:

function burgers_rhs(dut,ut,p,t)
# carefull here, kx has size (256,1) if you dont flatten it
nu,kx = p[1], vec(p[2])
u = ifft(ut)
# you need to put .= in order to modify dut inplace
dut .= -nu*kx.^2 .* ut - real.(fft(u.*ifft(1im*kx.*ut)))

end

2 Likes