5
$\begingroup$

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()

enter image description here

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.

enter image description here

The effect is even worse for double derivative.

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

enter image description here

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.

$\endgroup$
1

1 Answer 1

8
$\begingroup$

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.

$\endgroup$
10
  • 1
    $\begingroup$ @OverLordGoldDragon The answer is quite correct. The discrete Fourier transform is a trigonometric interpolation an thus suffers, like any other interpolation method, from Gibb's oscillations. Padding will improve the approximation a little, however won't fix it. $\endgroup$
    – ConvexHull
    Commented Aug 11, 2023 at 8:30
  • 1
    $\begingroup$ @OverLordGoldDragon I think it's pretty simple. Instead of using a trigonometric expansion for the second example, you might rather choose a hierarchical polynomial expansion (e.g. Chebyshev polynomials), which completely avoids the oscillations, since the approximation is able to represent the solution exactly. All depends on your chosen basis $\endgroup$
    – ConvexHull
    Commented Aug 11, 2023 at 10:15
  • 1
    $\begingroup$ @OverLordGoldDragon proper padding would indeed fix this problem, but applying a Fourier series truncation to a function that is discontinuous across it's periodic boundary is very well-studied and exhibits this phenomena, even with infinite precision. $\endgroup$
    – whpowell96
    Commented Aug 11, 2023 at 11:21
  • 1
    $\begingroup$ Note that the name of the person after whom Gibbs phenomenon was called is Gibbs, not "Gibb". $\endgroup$
    – Ruslan
    Commented Aug 11, 2023 at 14:53
  • 3
    $\begingroup$ @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. $\endgroup$
    – ConvexHull
    Commented Aug 11, 2023 at 15:12

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.