Closed
Description
Implementing a numerical scheme of ETDRK4, published in SIAM J.S.C in 2005, cited over 600 times. Trefethen is the God in spectral method. So the algorithm is in high trust. Matlab code is directly in the paper. So there should be no issue with the matlab code. Python code is what I copied from matlab and did some grammar correction. It is also proved to be match perfect. I did numerous sanity checks.
Matlab gives me stable result.
Scipy.fftpack gives me stable result.
Numpy.fft gives me blow up.
See the post here (full code is contained in that post ):
https://www.linkedin.com/pulse/weird-behavior-numpyfft-solving-chaotic-system-shaowu-shawn-pan/
The code with postprocessing is given below. Very short. Trust me.
#######################################################################
#
# Showing the difference between numpy.fft and scipy.fft
#
# Author: Shaowu Pan
#
# Date: 02/11/2018
#
#######################################################################
# import numpy and pyplot
import numpy as np
from matplotlib import pyplot as plt
import sys
# try to import numpyr or scipy fft based on argv
if sys.argv[1] == 'numpy':
## choose numpy
from numpy.fft import fft,ifft
else:
# choose scipy
from scipy.fftpack import fft, ifft
# define constant
pi = np.pi
# viscosity
nu = 1
# mesh
mesh = 128
# time restriction
tot = 30
# distributed x point
x = np.linspace(0, 2.0 * pi, mesh, endpoint=False)
# force IC
u0 = np.cos(x)*(1.0 + np.sin(x))
# set up mesh grid and k array
N = mesh
k_array_noshift = np.fft.fftfreq(mesh)*mesh
# set up initial condition
u = u0
v = fft(u)
# time step
h = 0.01
## follow them, set -N/2 to zer
k_array_noshift[mesh / 2] = 0
## Linear part in fourier space
L = 1.0/16.0*k_array_noshift**2 - nu*1.0/16.0**3 * k_array_noshift**4
## Fourier mulitplier
E = np.exp(h * L)
E2 = np.exp(h * L / 2.0)
## select number of points on the circle
M = 16
## choose radius 1, since k^2-k^4 ranges a lot... so 1 is enough to make sure only one singular point
r = np.exp(1j * np.pi * (np.arange(1, M + 1) - 0.5) / M)
r = r.reshape(1, -1)
r_on_circle = np.repeat(r, mesh, axis=0)
## define hL on M copies
LR = h * L
LR = LR.reshape(-1, 1)
LR = np.repeat(LR, M, axis=1)
## obtain hL on circle
LR = LR + r_on_circle
## important quantites used in ETDRK4
# g = -0.5*i*k
g = -0.5j * k_array_noshift
# averaged Q, f1,f2,f3
Q = h*np.real(np.mean( (np.exp(LR/2.0)-1)/LR, axis=1 ))
f1 = h*np.real(np.mean( (-4.0-LR + np.exp(LR)*(4.0-3.0*LR+LR**2))/LR**3, axis=1 ))
f2 = h*np.real(np.mean( (2.0+LR + np.exp(LR)*(-2.0 + LR))/LR**3, axis=1 ))
f3 = h*np.real(np.mean( (-4.0-3.0*LR - LR**2 + np.exp(LR)*(4.0 - LR))/LR**3, axis=1 ))
def compute_u2k_zeropad_dealiased(uk_):
"""
3/2 zero padding de-aliasing
"""
# u2k = fft(np.real(ifft(uk_))**2)
# three over two law
N = uk_.size
# map uk to uk_fine
uk_fine = np.hstack((uk_[0:N / 2], np.zeros((N / 2)), uk_[-N / 2:])) * 3.0 / 2.0
# convert uk_fine to physical space
u_fine = np.real(ifft(uk_fine))
# compute square
u2_fine = np.square(u_fine)
# compute fft on u2_fine
u2k_fine = fft(u2_fine)
# convert u2k_fine to u2k
u2k = np.hstack((u2k_fine[0:N / 2], u2k_fine[-N / 2:])) / 3.0 * 2.0
return u2k
print 'dt =', h
# np.linalg.norm(np.fft.ifft(uk_0)-u0) # very good
ntsnap = int(tot/h)
isnap = 0
tsnap = np.linspace(0, tot, ntsnap)
usnap = np.zeros((mesh, ntsnap))
uksnap = np.zeros((mesh, ntsnap),dtype=complex)
aasnap = np.zeros((mesh, ntsnap),dtype=complex)
bbsnap = np.zeros((mesh, ntsnap),dtype=complex)
ccsnap = np.zeros((mesh, ntsnap),dtype=complex)
Nvsnap = np.zeros((mesh, ntsnap),dtype=complex)
Nasnap = np.zeros((mesh, ntsnap),dtype=complex)
Nbsnap = np.zeros((mesh, ntsnap),dtype=complex)
Ncsnap = np.zeros((mesh, ntsnap),dtype=complex)
tcur = 0.0
## main loop time stepping
while tcur <= tot and isnap < ntsnap:
# print tcur
# record snap
# if abs(tcur - tsnap[isnap]) < 1e-2:
# print ' current progress =', tcur / tsnap[-1] * 100, ' % '
u = np.real(ifft(v))
usnap[:, isnap] = u
uksnap[:, isnap] = v
# Nv = g * fft(np.real(ifft(v))**2)
Nv = g * compute_u2k_zeropad_dealiased(v)
a = E2 * v + Q * Nv
# Na = g * fft(np.real(ifft(a))**2)
Na = g * compute_u2k_zeropad_dealiased(a)
b = E2 * v + Q * Na
# Nb = g * fft(np.real(ifft(b))**2)
Nb = g * compute_u2k_zeropad_dealiased(b)
c = E2 * a + Q * (2.0*Nb - Nv)
# Nc = g * fft(np.real(ifft(c))**2)
Nc = g * compute_u2k_zeropad_dealiased(c)
v = E*v + Nv*f1 + 2.0*(Na + Nb)*f2 + Nc*f3
aasnap[:, isnap] = a
bbsnap[:, isnap] = b
ccsnap[:, isnap] = c
Nvsnap[:, isnap] = Nv
Nasnap[:, isnap] = Na
Nbsnap[:, isnap] = Nb
Ncsnap[:, isnap] = Nc
# record time
tcur = tcur + h
isnap = isnap + 1
# save uksnap
np.savez('output_uk',uksnap=uksnap,f1=f1,f2=f2,f3=f3,Q=Q,LR=LR,L=L,g=g, Nasnap=Nasnap, Nvsnap=Nvsnap, Nbsnap=Nbsnap, Ncsnap=Ncsnap,
aasnap=aasnap, bbsnap=bbsnap, ccsnap=ccsnap, usnap=usnap)
# plot contours
from matplotlib import cm
fig = plt.figure(figsize=(8, 8))
X, Y = np.meshgrid(x, tsnap)
Z = usnap.transpose()
V = np.linspace(-10, 10, 500)
surf = plt.contourf(X, Y, Z, V, cmap=cm.seismic)
plt.savefig('contour.png')