3
$\begingroup$

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,

chaos goal

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.

choas fail

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 ?

$\endgroup$
4
  • 4
    $\begingroup$ Something suspicious about applying Fourier transform to nonlinear operator, that should lead to a convolution. $\endgroup$ May 3 at 17:38
  • 1
    $\begingroup$ 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. $\endgroup$ May 3 at 18:30
  • $\begingroup$ 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. $\endgroup$
    – bpachev
    May 3 at 19:34
  • $\begingroup$ 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}$$ $\endgroup$ May 5 at 15:45
5
$\begingroup$

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
  # Cranck-Nicholson + Adam 
  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}")

And I got similar answer:

enter image description here

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.