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

**URL:** https://discourse.julialang.org/t/issue-with-ode-solver-for-solving-burgers-equation-using-spectral-method/26224
**Category:** Modelling & Simulations
**Tags:** diffeq
**Created:** [July 10, 2019, 3:03pm UTC](https://discourse.julialang.org/t/issue-with-ode-solver-for-solving-burgers-equation-using-spectral-method/26224 "2019-07-10T15:03:42Z")
**Posts on this page:** 2
**Page:** 1

<div class="post-metadata">

### Author: ![henry2004y](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/henry2004y/32/9284_2.png) [@henry2004y](https://discourse.julialang.org/u/henry2004y)
#### Post date: [July 10, 2019, 3:03pm UTC](https://discourse.julialang.org/t/issue-with-ode-solver-for-solving-burgers-equation-using-spectral-method/26224/1 "2019-07-10T15:03:42Z")

</div>

Hi all,

I am trying to solve the 1D Burgers equation [Burgers equation wiki](https://en.wikipedia.org/wiki/Burgers%27_equation) 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:

```julia
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:

```julia
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?

---

<div class="post-metadata">

### Author: ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)
#### Post date: [July 10, 2019, 3:26pm UTC](https://discourse.julialang.org/t/issue-with-ode-solver-for-solving-burgers-equation-using-spectral-method/26224/2 "2019-07-10T15:26:04Z")

</div>

as a starter:

```julia
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

```