1

I'd like to differentiate a real, periodic function on (0,2*pi) which is also symmetric about x=pi using a discrete Fourier transform. I have written a Python code which does this using a FFT/IFFT but this does not take into account the symmetry of the function and so is a bit wasteful.

(The overall aim is to make a pseudospectral fluid flow solver and the periodicity and symmetry in one of the directions should allow me to expand the variables in that direction using only the cosine part of a Fourier series)

I know I need to use a Discrete Cosine Transform (DCT) to do this but cannot work out what needs to be changed about my domain (x), wavenumber vector (k) and implementation of the DCT/IDCT save that the former two should be half the length.

import sympy as sp
import numpy as np
import matplotlib.pylab as plt
from scipy.fftpack import fft, ifft

# Number of grid points
N = 2**5 

# Test function to check results with (using SymPy)
w = 3.; X = sp.Symbol('x'); Y=sp.cos(w*X)

# Domain of regularly spaced points in [0,2pi)
x=(2*np.pi/N)*np.arange(0,N)

# Calc exact derivatives using SymPy then turn into functions
dY      = Y.diff(X)
d2Y     = dY.diff(X)
d3Y     = d2Y.diff(X)
f       = sp.lambdify(X, Y,'numpy')
df_ex   = sp.lambdify(X, dY, 'numpy')
d2f_ex  = sp.lambdify(X, d2Y, 'numpy')
d3f_ex  = sp.lambdify(X, d3Y, 'numpy')

# Wavenumber vector
k=np.hstack(( np.arange(0,N/2), 0, np.arange(-N/2+1,0) )); 
k2=k**2; k3=k**3;

# Trans. to Fourier domain, diff, then return to phyical space
F   = fft(f(x))
df  = np.real(ifft(1j*k*F))
d2f = np.real(ifft( -k2*F))
d3f = np.real(ifft(-1j*k3*F))

# Plot result
fh=plt.figure(figsize=(8,4)); ah=fh.add_subplot(111)
plt.plot(x,f(x),'b-',x,df_ex(x), 'r-',x,d2f_ex(x),'g-',x,d3f_ex(x),'k-')
plt.plot(x,df,'ro',x,d2f,'go',x,d3f,'ko')
plt.xlim([0,2*np.pi])
anothermh
  • 9,815
  • 3
  • 33
  • 52
EDK
  • 63
  • 5
  • *"I have written a Python code which does this using a FFT/IFFT..."* Just FYI: you can also use [`scipy.fftpack.diff`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.fftpack.diff.html) to do this. For an example, see http://scipy-cookbook.readthedocs.io/items/KdV.html – Warren Weckesser Apr 15 '17 at 22:08
  • I do not understand what the question is. – user3344003 Apr 16 '17 at 01:42
  • Hello @user3344033. The question is, how to modify the code in order to take advantage of the symmetry of the function about x=pi so that the full FFT need not be used. As I said, I think it requires the DCT but have not been able to implement this. – EDK Apr 16 '17 at 08:36
  • @Warren Weckesser: thanks for that, I'm surprised that I didn't come across that in my lengthy searches. It looks very useful. For others interested, the following online notes also give a bit of background on it for pseudospectral purposes [link](http://homes.ioc.ee/lints/pubs/r312_15.pdf). – EDK Apr 16 '17 at 08:45

0 Answers0