# Taking derivative using FFT

I would like to calculate derivative of a given function ( a 1D array) using Array. Here is the code

import numpy as np
import matplotlib.pyplot as plt

Lleft=-np.pi
Lright=np.pi
L=Lright-Lleft
nx = 128
dx=L/nx

x = np.linspace(start=Lleft, stop=Lright-dx, num=nx)

u0=np.sin(x) #initial_profile
uhat = (1 / nx) * np.fft.fftshift(np.fft.fft(u0))  #forward_fourier_transform
u0transformed = np.real(nx * np.fft.ifft(np.fft.ifftshift(uhat)))   #transforming back to initial condition in real space

k = np.arange(-nx/2, nx/2, 1)

uhat= 1j*(k)*(2*np.pi/L) * uhat  #taking the derivative
utransformed = np.real(nx * np.fft.ifft(np.fft.ifftshift(uhat)))   #transforming back to solution in real space

plt.plot(x,u0transformed, "b",  label = " u0 = Sin(x)")
plt.plot(x,u0,"k--" , label = " u0 backtrasformed ")

plt.plot(x,np.cos(x),"y", label = "derivative of u0 = Cos(x)")
plt.plot(x,utransformed,"r--",label="derivative of u0 using FFT = Cos(x)")

plt.legend()


For particular example of trigonometric functions, the code works well but if I choose some periodic polynomial function, what I have observed that, resultant derivatives have "Runge like phenomenon" at the boundaries. Near boundaries of the domain the function fluctuates a lot.

The effect is even worse for double derivative.

uhat= -(k)**2*(2*np.pi/L)**2 * uhat


I am not sure what is the reason behind this. This also makes me doubt my understanding of all the FFT implementation I've encountered for simulating ODEs.

I believe this stems from the fact that your function $$f(x) = x^2$$ does not have continuous derivatives once it is extended periodically like $$\tilde{f}(x) = f(x \ \mathrm{mod} \ 12 -6),$$ which manifests as the Fourier coefficients of $$f$$ being proportional to $$n^{-2}$$, $$f'$$ to $$n^{-1}$$ and $$f''$$ to $$1$$.
In your plots of $$f'(x) = 2x$$, we see the typical Gibbs phenomena, where there are oscillations near the discontinuity at $$x + 6 = 0 \ \mathrm{mod} \ 12$$. Also notice that your Fourier reconstruction of $$f'$$ tends downwards at $$x=6$$ and upwards at $$x=-6$$. With better resolution, we would see that both of these converge to $$\frac{1}{2}(f'(-6) + f'(6)) = 0.$$
As for $$f''(x)$$, we are now differentiating a discontinuous function at $$x = \pm 6$$, so there is no surprise that the Fourier series has trouble representing this. The big takeaway here is that Fourier (and indeed, most globally supported modal methods with smooth modes) have trouble representing discontinuities within the domain. It just happens to be on the boundary of a domain, but the Fourier series of a function is equivalent to the Fourier series of a periodic extension of the function, which can be discontinuous if the values on the boundary of the interval are not equal.
• @OverLordGoldDragon It would be nice, if you could show me an oscillation free approximation of $x^2$ especially for the first and second derivative with a global Fourier expansion. Commented Aug 11, 2023 at 15:12