Generalized fourier series and orthogonality

Let \(f:\mathbb{R}\rightarrow\mathbb{C}\), \(g:\mathbb{R}\rightarrow\mathbb{C}\), and \(w:\mathbb{R}\rightarrow\mathbb{C}\) be complex functions. For square-integrable¹ \(f\), \(g\), and \(w\), the inner product of \(f\) and \(g\) with weight function \(w\) over the interval \([a,b] \subseteq \mathbb{R}\) is² \[\begin{aligned} \langle f,g \rangle_w &= \int_a^b f(x) \overline{g}(x) w(x) \diff x \end{aligned}\] where \(\overline{g}\) denotes the complex conjugate of \(g\). The inner product of functions can be considered analogous to the inner (or dot) product of vectors.

The fourier series components can be found by a special property of the \(\sin\) and \(\cos\) functions called orthogonality. In general, functions \(f\) and \(g\) from above are orthogonal over the interval \([a,b]\) iff \[\begin{aligned} \langle f,g \rangle_w = 0 \end{aligned}\] for weight function \(w\). Similar to how a set of orthogonal vectors can be a basis for a vector space, a set of orthogonal functions can be a basis for a function space: a vector space of functions from one set to another (with certain caveats).

In addition to some sets of sinusoids, there are several other important sets of functions that are orthogonal. For instance, sets of legendre polynomials (Kreyszig 2011, sec. 5.2) and bessel functions (Kreyszig 2011, sec. 5.4) are orthogonal.

As with sinusoids, the orthogonality of some sets of functions allows us to compute their series components. Let functions \(f_0, f_1, \cdots\) be orthogonal with respect to weight function \(w\) on interval \([a,b]\) and let \(a_0, a_1, \cdots\) be real constants. A generalized fourier series is (Kreyszig 2011, sec. 11.6) \[\begin{aligned} \label{eq:gen_fourier_synthesis} f(x) = \sum_{m=0}^\infty a_m f_m(x) \end{aligned}\] and represents a function \(f\) as a convergent series. It can be shown that the Fourier components \(a_m\) can be computed from $$\begin{align} a_m = \frac{\langle f, f_m \rangle_w} {\langle f_m, f_m \rangle_w}.\label{eq:gen-fourier-analysis} \end{align}$$ In keeping with our previous terminology for fourier series, eq. ¿eq:gen_fourier_synthesis? and eq. ¿eq:gen-fourier-analysis? are called general fourier synthesis and analysis, respectively.

For the aforementioned legendre and bessel functions, the generalized fourier series are called fourier-legendre and fourier-bessel series (Kreyszig 2011, sec. 11.6). These and the standard fourier series (section 6.1) are of particular interest for the solution of partial differential equations (chapter 7).

Example 6.4

The Bessel functions of the first kind are series solutions to the Bessel differential equation [@kreyszig2011,§ 5.4], which models boundary value problems associated with circular membranes in polar coordinates. The nth Bessel function is defined as $$ J_\nu(x) = x^\nu \sum_{m=0}^\infty \frac{(-1)^m x^{2 m}}{2^{2 m + \nu} m! \Gamma(\nu + m + 1)}, $$ where we have used the following definitions:

• ν is a real number called the order of the Bessel function
• Γ is the gamma function defined as $\Gamma(\nu+1) = \int_0^\infty e^{-t} t^\nu \diff t$ for ν >  − 1. The gamma function is a generalization of the factorial function to real numbers.

The Bessel functions J_n(k_n, mx), for nonnegative integer n, are orthogonal with respect to the weight function w(x) = x on the interval [0,R] for R > 0 [@kreyzig2011,§ 11.6]. Here $$\begin{align} k_{n,m} = \frac{\alpha_{n,m}}{R}, \label{eq:bessel-orthogonality-k} \end{align}$$ where α_n, m are the zeros of the Bessel function J_n. This orthogonality allows us to compute the Fourier-Bessel series components of a function f(x) on [0,R] from [@eq:gen-fourier-analysis].

In this example, we will explore the definition of the Bessel function by plotting J₀, J₁, J_1/2, and J_−1/2. Furthermore, we will compute the first several Fourier-Bessel series components of the function f(x) = 1 − x³ on the interval x ∈ [0,1] and plot partial sums of the series.

We proceed in Python. First, load packages.

import numpy as np
import sympy as sp
import scipy
from scipy.special import jn_zeros
import matplotlib.pyplot as plt

Define the Bessel functions of the first kind of orders 0 and 1.

x = sp.symbols('x', nonnegative=True)
J_0 = sp.besselj(0, x)
J_1 = sp.besselj(1, x)

Plot these functions.

x_plt = np.linspace(0, 20, 1001)
J_0_fun = sp.lambdify(x, J_0, modules=['numpy', 'scipy'])
J_1_fun = sp.lambdify(x, J_1, modules=['numpy', 'scipy'])
fig, ax = plt.subplots()
ax.plot(x_plt, J_0_fun(x_plt), label='$J_0(x)$')
ax.plot(x_plt, J_1_fun(x_plt), label='$J_1(x)$')
ax.spines['bottom'].set_position('zero')
ax.spines['left'].set_position('zero')
ax.set_xlabel('$x$')
ax.set_ylabel('$J_n(x)$')
ax.legend()
plt.draw()

The plot shows the Bessel functions of the first kind of orders 0 and 1. The functions are oscillatory and decay to zero as x increases. They resemble sinusoids with decreasing amplitude.

Next, we will plot the Bessel functions of the first kind of orders 1/2 and  − 1/2.

J_half = sp.besselj(1/2, x)
J_half_n = sp.besselj(-1/2, x)

Plot these functions.

J_half_fun = sp.lambdify(x, J_half, modules=['numpy', 'scipy'])
J_half_n_fun = sp.lambdify(x, J_half_n, modules=['numpy', 'scipy'])
fig, ax = plt.subplots()
ax.plot(x_plt, J_half_fun(x_plt), label='$J_{1/2}(x)$')
ax.plot(x_plt, J_half_n_fun(x_plt), label='$J_{-1/2}(x)$')
ax.set_ylim(-1, 1)
ax.spines['bottom'].set_position('zero')
ax.spines['left'].set_position('zero')
ax.set_xlabel('$x$')
ax.set_ylabel('$J_n(x)$')
ax.legend()
plt.draw()

These functions are similar to the Bessel functions of orders 0 and 1. However, J_−1/2(x) has a singularity at x = 0.

Note that, like sinusoids, all these Bessel functions produce zeros at certain values of x. The zeros are crucial to the orthogonality of the Bessel functions, as the introduction to this example explains. Therefore, we must be able to compute these zeros. kreyszig2011 provides a table for the first several zeros of the Bessel functions. Alternatively, we can use the scipy package’s special module to compute the zeros, as follows for J₀(x).

n_zeros = 5
alpha_0 = jn_zeros(0, n_zeros)

Turning to the Fourier-Bessel series, we will consider the function f(x) = 1 − x³ on the interval [0,1]. The Fourier-Bessel series analysis of this function computes the coefficients a_m from [@eq:gen-fourier-analysis]. The coefficients are given by $$ a_m = \frac{\langle f, f_m\rangle_w}{\langle f_m, f_m\rangle_w}, $$ where f_m = J₀(k_0, mx) and k_0, m are given by [@eq:bessel-orthogonality-k] and w(x) = x is the weight function. The inner products are computed as $$\begin{align} \langle f, f_m\rangle_w &= \int_0^R f(x) f_m(x) w(x) \, dx, \\ \langle f_m, f_m\rangle_w &= \int_0^R f_m(x)^2 w(x) \, dx. \end{align}$$ Computing these integrals, we obtain the following:

R = 1
m = sp.symbols('m', integer=True, positive=True)
k_0_m = sp.symbols('k_0_m', real=True, positive=True)
f = 1 - x**3
f_m = J_0.subs(x, k_0_m * x)
w = x
inner_f_f_m = sp.integrate(f * f_m * w, (x, 0, R))
inner_f_m_f_m = sp.integrate(f_m**2 * w, (x, 0, R))
a_m = inner_f_f_m / inner_f_m_f_m

Convert the symbolic expression for a_m to a numerical function of k_0, m and compute the coefficients for the first several values of k_0, m.

a_m_fun = sp.lambdify(k_0_m, a_m, modules=['numpy', 'mpmath'])
a_m_ = np.zeros_like(alpha_0)
for i, alpha in enumerate(alpha_0):  # Vectorization fails; workaround
    a_m_[i] = a_m_fun(alpha)
print(a_m_)

[ 1.27215302 -0.33478735  0.09593875 -0.04835341  0.02549602]

Next, we compute the partial sum of the Fourier-Bessel series up to m = 5. The partial sum is given by $$ f_N(x) = \sum_{m=1}^N a_m J_0(k_{0,m} x). $$ We will evaluate this sum numerically.

N = n_zeros
f_N = 0
for i in range(N):
    f_N += a_m_[i] * sp.besselj(0, alpha_0[i] * x)
f_N_fun = sp.lambdify(x, f_N, modules=['numpy', 'scipy'])

Plot the function f(x) and the partial sum f_N(x).

x_plt = np.linspace(0, R, 101)
fig, ax = plt.subplots()
ax.plot(x_plt, 1 - x_plt**3, label='$f(x)$')
ax.plot(x_plt, f_N_fun(x_plt), 
        label='Partial sum $f_N(x)$', linestyle='--')
ax.spines['bottom'].set_position('zero')
ax.spines['left'].set_position('zero')
ax.set_xlabel('$x$')
ax.set_ylabel('$f(x)$')
ax.legend()
plt.show()

The plot shows the function f(x) = 1 − x³ and the partial sum of the Fourier-Bessel series up to m = 5. The partial sum approximates the function quite well, capturing its general shape.