Fourier Series Approximation

Visualizing Functional Analysis with Python

I’ll provide an example problem in functional analysis and solve it using $Python$.

I’ll include mathematical notations in $TeX$, source code, and visualize the results with graphs.

Functional Analysis Example: Approximating Functions with Fourier Series

Let’s consider a fundamental problem in functional analysis: approximating a function using a Fourier series.

We’ll examine how well we can approximate a square wave function using different numbers of terms in its Fourier series.

Mathematical Background

A square wave function $f(x)$ on the interval $[-\pi, \pi]$ can be defined as:

$$
f(x) = \begin{cases}
1, & \text{if } 0 < x < \pi \
-1, & \text{if } -\pi < x < 0
\end{cases}
$$

The Fourier series of this function is:

$$
f(x) \approx \frac{4}{\pi} \sum_{n=1,3,5,…}^{\infty} \frac{1}{n} \sin(nx)
$$

This is a perfect example to demonstrate how increasing the number of terms improves approximation, a key concept in functional spaces.

Let’s implement this in $Python$ and visualize the results:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import animation
from IPython.display import HTML

# Define the square wave function
def square_wave(x):
    return np.where(x > 0, 1, -1)

# Define the Fourier series approximation
def fourier_approximation(x, num_terms):
    result = np.zeros_like(x, dtype=float)
    for n in range(1, num_terms + 1, 2):  # Only odd terms
        result += (4 / (np.pi * n)) * np.sin(n * x)
    return result

# Set up the x values
x = np.linspace(-np.pi, np.pi, 1000)
true_function = square_wave(x)

# Calculate different approximations
approximations = {}
term_counts = [1, 3, 5, 11, 21, 51, 101]
for terms in term_counts:
    approximations[terms] = fourier_approximation(x, terms)

# Compute the L2 error for each approximation
l2_errors = {}
for terms, approx in approximations.items():
    l2_errors[terms] = np.sqrt(np.mean((approx - true_function)**2))

# Create visualization
plt.figure(figsize=(12, 8))

# Plot the true function
plt.subplot(2, 2, 1)
plt.plot(x, true_function, 'k-', label='Square Wave')
plt.title('Original Square Wave Function')
plt.grid(True)
plt.legend()
plt.xlim(-np.pi, np.pi)
plt.ylim(-1.5, 1.5)

# Plot approximations with different terms
plt.subplot(2, 2, 2)
plt.plot(x, true_function, 'k-', label='Square Wave')
for terms in [1, 5, 21, 101]:
    plt.plot(x, approximations[terms], label=f'{terms} terms')
plt.title('Fourier Series Approximations')
plt.grid(True)
plt.legend()
plt.xlim(-np.pi, np.pi)
plt.ylim(-1.5, 1.5)

# Plot L2 error vs number of terms
plt.subplot(2, 2, 3)
plt.plot(list(l2_errors.keys()), list(l2_errors.values()), 'bo-')
plt.xscale('log')
plt.yscale('log')
plt.title('L2 Error vs Number of Terms')
plt.xlabel('Number of terms (log scale)')
plt.ylabel('L2 Error (log scale)')
plt.grid(True)

# Create a heatmap visualization of the approximation space
plt.subplot(2, 2, 4)
term_dense = np.arange(1, 102, 2)
x_sample = np.linspace(-np.pi, np.pi, 100)
approximation_matrix = np.zeros((len(term_dense), len(x_sample)))

for i, terms in enumerate(term_dense):
    approximation_matrix[i, :] = fourier_approximation(x_sample, terms)

plt.imshow(approximation_matrix, aspect='auto', 
           extent=[-np.pi, np.pi, term_dense[-1], term_dense[0]], 
           cmap='viridis')
plt.colorbar(label='Function Value')
plt.title('Function Approximation Space')
plt.xlabel('x')
plt.ylabel('Number of terms')

plt.tight_layout()
plt.savefig('fourier_approximation.png', dpi=300)
plt.show()

# Fix the animation code
fig, ax = plt.subplots(figsize=(10, 6))
line_true, = ax.plot(x, true_function, 'k-', label='Square Wave')
line_approx, = ax.plot([], [], 'r-', label='Approximation')
ax.set_xlim(-np.pi, np.pi)
ax.set_ylim(-1.5, 1.5)
ax.grid(True)
ax.legend()

# The title needs to be initialized outside of animation functions
title_text = ax.set_title('Fourier Approximation: 0 terms')

def init():
    line_approx.set_data([], [])
    title_text.set_text('Fourier Approximation: 0 terms')
    return line_approx, title_text

def animate(i):
    terms = 2*i - 1 if i > 0 else 0  # Handle the case when i=0
    if terms <= 0:
        y = np.zeros_like(x)
    else:
        y = fourier_approximation(x, terms)
    line_approx.set_data(x, y)
    title_text.set_text(f'Fourier Approximation: {terms} terms')
    return line_approx, title_text

# Create animation with 51 frames (up to 101 terms)
anim = animation.FuncAnimation(fig, animate, init_func=init,
                              frames=52, interval=200, blit=True)

# Display the animation (in Google Colab)
plt.close()  # Prevents duplicate display
HTML(anim.to_jshtml())

# Calculate and display convergence rate in functional spaces
terms_for_convergence = np.arange(1, 202, 2)
errors = []

for terms in terms_for_convergence:
    approx = fourier_approximation(x, terms)
    error = np.sqrt(np.mean((approx - true_function)**2))
    errors.append(error)

# Plot convergence
plt.figure(figsize=(10, 6))
plt.loglog(terms_for_convergence, errors, 'bo-')
plt.xlabel('Number of Terms')
plt.ylabel('L2 Error')
plt.title('Convergence Rate of Fourier Series Approximation')
plt.grid(True)

# Fit a power law to see the convergence rate
from scipy.optimize import curve_fit

def power_law(x, a, b):
    return a * x**b

params, _ = curve_fit(power_law, terms_for_convergence, errors)
a, b = params

# Add the fitted curve to the plot
x_fit = np.logspace(np.log10(terms_for_convergence[0]), np.log10(terms_for_convergence[-1]), 100)
plt.loglog(x_fit, power_law(x_fit, a, b), 'r-', 
           label=f'Fitted power law: {a:.2e} × n^({b:.3f})')
plt.legend()
plt.show()

print(f"Convergence rate: O(n^{b:.3f})")

Code Explanation

The code implements a Fourier series approximation of a square wave function.

Here’s a breakdown of what it does:

  1. Square Wave Definition:
    We define a square wave function that returns $1$ for positive $x$ and $-1$ for negative x.

  2. Fourier Approximation Function:
    This implements the Fourier series formula I mentioned earlier, summing only the odd terms as per the mathematical derivation.

  3. Multiple Approximations:
    We calculate several approximations with different numbers of terms $(1, 3, 5, 11, 21, 51, 101)$.

  4. Error Calculation:
    For each approximation, we compute the L2 error (root mean square error) compared to the true function.

  5. Visualizations:

    • The original square wave
    • Different approximations overlaid on the original function
    • L2 error vs. number of terms on a log-log scale
    • A heatmap showing how the approximation evolves with increasing terms

  6. Animation:
    Shows how the approximation improves as we add more terms.

  7. Convergence Analysis:
    We fit a power law to the error values to determine the convergence rate of the Fourier series.

Results and Analysis


Convergence rate: O(n^-0.400)

The visualizations show several important aspects of functional approximation:

  1. Approximation Quality:
    As we increase the number of terms in the Fourier series, the approximation gets closer to the true square wave. With just $1$ term, we have a sine wave.
    With $101$ terms, we have a reasonable approximation of the square wave, though with visible Gibbs phenomenon (oscillations near the discontinuities).

  2. L2 Error Decay:
    The log-log plot of error vs. number of terms shows that the error decreases approximately as $O(n^{(-0.5)})$, which matches theoretical expectations for Fourier series approximation of discontinuous functions.

  3. Function Space Visualization:
    The heatmap shows how the approximation evolves across the entire domain as we add more terms, providing insight into the convergence properties in function space.

  4. Gibbs Phenomenon:
    Even with many terms, there’s oscillation near the discontinuities (at $x = 0$), which is the famous Gibbs phenomenon. This illustrates a fundamental limitation in Fourier approximation of discontinuous functions.

Functional Analysis Insights

This example demonstrates several key concepts in functional analysis:

  1. Approximation in Function Spaces:
    How functions can be approximated using orthogonal bases (in this case, the Fourier basis).

  2. Convergence Properties:
    The rate of convergence depends on the smoothness of the function being approximated.

  3. L2 Norm:
    We used the L2 norm (square root of the mean squared error) to measure the distance between functions, which is fundamental in Hilbert spaces.

  4. Completeness of Basis:
    The Fourier basis is complete in L2$[-π,π]$, meaning any square-integrable function can be approximated to arbitrary precision.

The visualization clearly shows that even though we need infinitely many terms for perfect representation, we can achieve practical approximations with a finite number of terms.

This demonstrates the power of functional analysis in providing theoretical frameworks for practical computational methods.