# Solving numerically the 1D Kuramoto-Sivashinsky Equation using spectral methods

Lately, I've been trying to solve numerically the 1D Kuramoto-Sivashinsky Equation using spectral methods.

Let $$\nu$$ be the viscosity and $$[0,L]$$ the domain. The 1D equation is, $$u_t + uu_x + u_{xx} + \nu u_{xxxx} = 0 \iff u_t = - uu_x - (u_{xx} + \nu u_{xxxx})$$

Introducing the linear operator $$\mathcal{L}u = - (u_{xx} + \nu u_{xxxx})$$ and the non-linear operator $$\mathcal{N}u = - u u_x = - \frac{1}{2} u^2_{x}$$ We can write this PDE as a sum of these two operators $$u_t = \mathcal{L}u + \mathcal{N}u$$

We can solve the linear part with a Crank-Nicholson scheme $$\begin{equation} \frac{u(x,t + \Delta t) - u(x,t)}{\Delta t} = \frac{1}{2} \left( \mathcal{L}u(x,t + \Delta t) + \mathcal{L}u(x,t) \right) \end{equation}$$

And we can solve the non-linear part with an Adam Bashforth scheme $$\begin{equation} \frac{u(x,t + \Delta t) - u(x,t)}{\Delta t} = \frac{3}{2} \mathcal{N}u(x,t) - \frac{1}{2} \mathcal{N} u(x,t - \Delta t) \end{equation}$$

After some algebra and going into Fourier space, we find $$\begin{equation} \hat{u}(k, t + \Delta t) = \frac{1}{1 - \frac{\Delta t}{2} \mathcal{F}[\mathcal{L}]} \left( \left[1 + \frac{\Delta t}{2} \mathcal{F}[\mathcal{L}] \right] \hat{u}(k,t) + \left[ \frac{3}{2} \mathcal{F}[\mathcal{N}] (\hat{u}(k,t))^2 - \frac{1}{2} \mathcal{F} [\mathcal{N}] (\hat{u}(k,t - \Delta t))^2 \right] \Delta t \right) \end{equation}$$

where the Fourier transform of the linear and non-linear operator is given by \begin{align*} \mathcal{F}[\mathcal{L}] &= - \left(i \frac{2\pi}{L} k \right)^2 - \nu \left(i \frac{2\pi}{L} k \right)^4 = \left( \frac{2\pi}{L} k \right)^2 - \nu \left( \frac{2\pi}{L} k \right)^4 \\ \mathcal{F}[\mathcal{N}u(x,t)] &= - \frac{1}{2} \left(i \frac{2\pi}{L} k \right) \hat{u}^2(k,t) \end{align*}

I'm trying to "show" that for $$L = 100$$ we have chaos. Taking $$\nu = 1$$, $$x \in [0, L]$$ where $$L = 100$$, $$t \in [0, 200]$$ and the initial condition, $$\begin{equation} u(x,0) = \cos(2\pi x / L) + 0.1 \cos(4\pi x / L) \end{equation}$$

We're also discretizing the domain in $$N = 1024$$ points and taking a step time $$\Delta t = 0.05$$.

The 1D plot I'm trying to obtain is the following, However, the one I obtain is different. I guess I'm not too far away of obtaining the good one but the chaos is not very much visible in mine. And here is the code:

import numpy as np
import matplotlib.pyplot as plt

from matplotlib.pyplot import cm

nu = 1
L = 100
nx = 1024

t0 = 0
tN = 200
dt = 0.05
nt = int((tN - t0) / 0.05)

# wave number mesh
k = np.arange(-nx/2, nx/2, 1)

t = np.linspace(start=t0, stop=tN, num=nt)
x = np.linspace(start=0, stop=L, num=nx)

# solution mesh in real space
u = np.ones((nx, nt))
# solution mesh in Fourier space
u_hat = np.ones((nx, nt), dtype=complex)

# initial condition
u0 = np.cos((2 * np.pi * x) / L) + 0.1 * np.cos((4 * np.pi * x) / L)

# Fourier transform of initial condition
u0_hat = (1 / nx) * np.fft.fftshift(np.fft.fft(u0))

# set initial condition in real and Fourier mesh
u[:,0] = u0
u_hat[:,0] = u0_hat

# Fourier Transform of the linear operator
FL = (((2 * np.pi) / L) * k) ** 2 - nu * (((2 * np.pi) / L) * k) ** 4
# Fourier Transform of the non-linear operator
FN = - (1 / 2) * ((1j) * ((2 * np.pi) / L) * k)

# resolve EDP in Fourier space
for j in range(0,nt-1):
uhat_current = u_hat[:,j]
if j == 0:
uhat_last = u_hat[:,0]
else:
uhat_last = u_hat[:,j-1]

# compute solution in Fourier space through a finite difference method
u_hat[:,j+1] = (1 / (1 - (dt / 2) * FL)) * ( (1 + (dt / 2) * FL) * uhat_current + ( ((3 / 2) * FN) * (uhat_current ** 2) - ((1 / 2) * FN) * (uhat_last ** 2) ) * dt )
# go back in real space
u[:,j+1] = np.real(nx * np.fft.ifft(np.fft.ifftshift(u_hat[:,j+1])))

# plot the result
fig, ax = plt.subplots(figsize=(10,8))

xx, tt = np.meshgrid(x, t)
levels = np.arange(-3, 3, 0.01)
cs = ax.contourf(xx, tt, u.T, cmap=cm.jet)
fig.colorbar(cs)

ax.set_xlabel("x")
ax.set_ylabel("t")
ax.set_title(f"Kuramoto-Sivashinsky: L = {L}, nu = {nu}")


After checking multiple times computations and code, I don't see where I'm wrong ? Maybe there's some non-linear effect I should manage ?

• Something suspicious about applying Fourier transform to nonlinear operator, that should lead to a convolution. May 3 at 17:38
• Indeed. I applied a pseudo-spectral method to deal with non-linearities. No I have the good plot. Gonna edit my post with this solution. May 3 at 18:30
• Do you just need a KSE solver or are you writing your own one for the learning experience? Because github.com/jswhit/pyks is another pseudospectral solver I've used in research. May 3 at 19:34
• Well I believe @MaximUmansky is right here and the Fourier transform of your nonlinear term is indeed wrong. I think you need to know that: $$\mathcal{F} [ (f(x))^{2} ] \neq (\mathcal{F} [ f(x) ])^{2}$$ May 5 at 15:45

Indeed the problem was that you were trying to calculate your nonlinear term incorrectly and you forgot that:

$$\mathcal{F}[(f(x))^{2}] \neq (\mathcal{F}[f(x)])^{2}$$

I changed your code to this:

import numpy as np
import matplotlib.pyplot as plt

from matplotlib.pyplot import cm

nu = 1
L = 100
nx = 1024

t0 = 0
tN = 200
dt = 0.05
nt = int((tN - t0) / 0.05)

# wave number mesh
k = np.arange(-nx/2, nx/2, 1)

t = np.linspace(start=t0, stop=tN, num=nt)
x = np.linspace(start=0, stop=L, num=nx)

# solution mesh in real space
u = np.ones((nx, nt))
# solution mesh in Fourier space
u_hat = np.ones((nx, nt), dtype=complex)

u_hat2 = np.ones((nx, nt), dtype=complex)

# initial condition
u0 = np.cos((2 * np.pi * x) / L) + 0.1 * np.cos((4 * np.pi * x) / L)

# Fourier transform of initial condition
u0_hat = (1 / nx) * np.fft.fftshift(np.fft.fft(u0))

u0_hat2 = (1 / nx) * np.fft.fftshift(np.fft.fft(u0**2))

# set initial condition in real and Fourier mesh
u[:,0] = u0
u_hat[:,0] = u0_hat

u_hat2[:,0] = u0_hat2

# Fourier Transform of the linear operator
FL = (((2 * np.pi) / L) * k) ** 2 - nu * (((2 * np.pi) / L) * k) ** 4
# Fourier Transform of the non-linear operator
FN = - (1 / 2) * ((1j) * ((2 * np.pi) / L) * k)

# resolve EDP in Fourier space
for j in range(0,nt-1):
uhat_current = u_hat[:,j]
uhat_current2 = u_hat2[:,j]
if j == 0:
uhat_last = u_hat[:,0]
uhat_last2 = u_hat2[:,0]
else:
uhat_last = u_hat[:,j-1]
uhat_last2 = u_hat2[:,j-1]

# compute solution in Fourier space through a finite difference method
u_hat[:,j+1] = (1 / (1 - (dt / 2) * FL)) * ( (1 + (dt / 2) * FL) * uhat_current + ( ((3 / 2) * FN) * (uhat_current2) - ((1 / 2) * FN) * (uhat_last2) ) * dt )
# go back in real space
u[:,j+1] = np.real(nx * np.fft.ifft(np.fft.ifftshift(u_hat[:,j+1])))
u_hat2[:,j+1] = (1 / nx) * np.fft.fftshift(np.fft.fft(u[:,j+1]**2))

# plot the result
fig, ax = plt.subplots(figsize=(10,8))

xx, tt = np.meshgrid(x, t)
levels = np.arange(-3, 3, 0.01)
cs = ax.contourf(xx, tt, u.T, cmap=cm.jet)
fig.colorbar(cs)

ax.set_xlabel("x")
ax.set_ylabel("t")
ax.set_title(f"Kuramoto-Sivashinsky: L = {L}, nu = {nu}") 