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
# Cranck-Nicholson + Adam
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 ?
