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
  • 5
    $\begingroup$ Something suspicious about applying Fourier transform to nonlinear operator, that should lead to a convolution. $\endgroup$ May 3, 2021 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, 2021 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, 2021 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, 2021 at 15:45

2 Answers 2

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$
0
$\begingroup$

Thanks for the code but when I try the code with T = 1000, it is not stable. When I add

u_hat[:,j+1] = (1 / nx) * np.fft.fftshift(np.fft.fft(u[:,j+1]))

which removes the accumulated error in the imaginary part, I can get stable results for T=1000

The code becomes

import numpy as np
import matplotlib.pyplot as plt

from matplotlib.pyplot import cm

nu = 1
L = 100 
nx = 1024

t0 = 0 
tN = 1000
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])))

  # clean the imaginary part contribution in u_hat
  u_hat[:,j+1] = (1 / nx) * np.fft.fftshift(np.fft.fft(u[:,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}")

enter image description here

$\endgroup$
2
  • $\begingroup$ It's not the answer to the question. It wasn't meant to be stable for any $t$ and OP didn't want to run it for $t > 200$. $\endgroup$ Jun 13, 2022 at 21:15
  • $\begingroup$ As it’s currently written, your answer is unclear. Please edit to add additional details that will help others understand how this addresses the question asked. You can find more information on how to write good answers in the help center. $\endgroup$
    – Community Bot
    Jun 13, 2022 at 21:15

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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