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

I'm having trouble with the following implementation of the KS model (see below) found on

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

When tN > 300 an overflow occurs on the line:

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

I tried using a workaround with the convolution theorem with no use, any help would be much apreciated.

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}")$$$$

• Are the inputs finite? Dec 6, 2021 at 13:37
• Yes they are, time goes from t0 to tN, and x goes from 0 to L with periodic boundary conditions Dec 7, 2021 at 1:05
• I meant the u0 variable in your code, or the u variable in the one line you show above. If u0 is already very (unreasonably) large, then it may not surprise that the FFT overflows. Dec 7, 2021 at 17:25
• It is not going to be larger than 1024 elements, u0 is of size nx`, which in this case is 1024, but I also tried with only 64 elements and also end up having the overflow Dec 7, 2021 at 23:07
• No, I was wondering about the size of the elements of that vector, not its length. Dec 8, 2021 at 5:09