Differentiation of continuous functions using the FFT

Differentiation of continuous functions using the FFT#

In this notebook, we use the FFTArray to compute the derivatives of two periodic test functions. Afterwards we will evaluate these functions on a one-dimensional grid and compare the results to their analytic counterparts.

To that end, we briefly recall the definitions for differentiation via Fourier transform from (Wikipedia):

Suppose \(g(x)\) with position \(x {\in} {\mathbb{R}}\) is an absolutely continuous differentiable function, and both \(g\) and its derivative \(g′\) are integrable. Then, the Fourier transform of the derivative is given by: \begin{equation*} \mathcal{F}\left\{ \frac{d}{dx} g(x)\right\} = i 2\pi f \ \hat{g}(f). \end{equation*} Note, that the differential operator transformed into a multiplication with the frequencies \(f {\in} {\mathbb{R}}\).

This in principle allows to compute the derivative of any function \(g(x)\) by first computing its Fourier transformation and then inversely transforming the product: \begin{equation*} \frac{d^n}{d^n x} g(x) = \mathcal{F^{-1}}{ \left( \mathcal{F}{ \left( \frac{d^n}{d^n x} g(x) \right)} \right)} = \mathcal{F^{-1}}{ \left( (i 2\pi f)^n \ \hat{g}(f) \right)}. \end{equation*}

However, finding the analytic expressions for the Fourier transform and its inverse by computing double integrals is often not straightforward or even impossible. Instead, numerically sampling \(g(x)\) allows us to apply the FFT to solve this problem.

Note: To avoid errors when evaluating the above equations numerically due to sampling on a finite grid, we require that \(g(x) \rightarrow 0\) sufficiently fast as \(x\) approaches the grid boundaries. This can be ensured, for example, via sufficiently large grids. Additionally the sampling step size is small enough to properly represent all frequencies contained in the function.

Definition and implementation of periodic test function#

We define the following test function \(g(x)\) and compute its derivativites \(g^{(n)}(x) = \frac{\mathrm{d}^{(n)}}{\mathrm{d}x^{(n)}}\):

\begin{align} g(x) &= \cos{(x)}\, e^{-(x-a)^{2}/25} \\ g'(x) &= - [\sin{(x)} + \frac{2}{25}(x-a) \cos{(x)}]\, e^{-(x-a)^{2}/25} \\ g''(x) &= \{[ \frac{1}{625} (4(x-a)^2-50)-1 ]\cos{(x)} + \frac{4}{25} (x-a) \sin{(x)} \}\, e^{-(x-a)^{2}/25} \\ \end{align}

The analytical Fourier transform of \(g(x)\) is given by: \begin{align*} G(f) &=\int_{-\infty}^{\infty}dx \ g(x)\ e^{- 2 \pi i f x},\quad \forall\ x \in \mathbb R . \end{align*}

Definitions and numerical sampling#

We proceed to implement the test function and its analytic derivative which we will later compare to the result of the differentiation via the FFT:

[1]:

# required packages
import numpy as np
from bokeh.io import output_notebook

import fftarray as fa

# plotting functions
from helpers import plt_deriv_comparison, plt_deriv_sampling
output_notebook(hide_banner=True)


# Analytical definition of test function and its derivatives.
# test function in position space
def g_x(x, a):
    return fa.cos(x)*fa.exp((-(x - a)**2)/25.)
# first derivative of the test function
def g_d1_x(x, a):
    return -(fa.sin(x)+ (2*(x-a)*fa.cos(x))/(25.))*fa.exp(-(x-a)**2/25.)
# second derivative of the test function
def g_d2_x(x, a):
    return (((1./625.)*(4*(x-a)**2-50)-1)*fa.cos(x)+(4./25.)*(x-a)*fa.sin(x))*fa.exp(-(x - a)**2/25.)

# We sample the test function on a finite one-dimensional and plot it using FFTarray.


# initialize finite coordinate grid and backend of FFTarray.
x_dim = fa.dim_from_constraints(
    "x",              # dimension name
    n=2048,           # number of grid points
    pos_middle=0.,    # center of position grid
    pos_extent=80.,   # extent of position grid
    freq_middle=0.,   # center of frequency grid
)

# get Array from Dimension in position space.
x = fa.coords_from_dim(x_dim, "pos")
# numerically sampled test function in position space.
a=1.25
g_x_analytic = g_x(a=a, x=x)
g_d1_x_analytic = g_d1_x(a=a, x=x)
g_d2_x_analytic = g_d2_x(a=a, x=x)

# plot analytic test functions. For plot parameters see helpers.py.
plt_deriv_sampling("Plotting analytic test functions",g_x_analytic,g_d1_x_analytic,g_d2_x_analytic)

Numerical differentiation and residual analysis#

Finally, we implement the formula for computing the derivative via the FFT and compare the results with the analytic solutions:

[2]:

# Implementation of the RHS of the second euqation of this notebook
def derivative_pos(dim: fa.Dimension, arr: fa.Array, order: int) -> fa.Array:
    """Takes the derivative of order `order` along Dimension `dim` of the
    Array `arr` in position space.
    """
    kernel = (1.j*2*np.pi*fa.coords_from_dim(dim,"freq"))**order
    return (kernel*arr.into_space("freq")).into_space("pos")

# first order derivative
g_d1_x_numeric = derivative_pos(x_dim, g_x_analytic, 1) # numerical result

# second order derivative
g_d2_x_numeric = derivative_pos(x_dim, g_x_analytic, 2) # numerical result

# plot comparison as well as residual differences. For plot parameters see helpers.py
plt_deriv_comparison("First Order Derivative",g_d1_x_analytic,"g'(x)",g_d1_x_numeric,"g'_num(x)")
plt_deriv_comparison("Second Order Derivative",g_d2_x_analytic,"g''(x)",g_d2_x_numeric,"g''_num(x)")

# machine confirmation that results are almost equal to analytic solutions.
np.testing.assert_array_almost_equal(g_d1_x_analytic.values("pos", xp=np), g_d1_x_numeric.values("pos", xp=np), decimal=11)
np.testing.assert_array_almost_equal(g_d2_x_analytic.values("pos", xp=np), g_d2_x_numeric.values("pos", xp=np), decimal=9)

Summary#

The two comparisons show good agreement between the analytical test functions and our numerical results. We have thus demonstrated in this notebook how straightforward it is to numerically compute the derivatives of continous functions using FFTArray.

Contents