Text provided under a Creative Commons Attribution license, CC-BY. All code is made available under the FSF-approved MIT license. (c) Marc Spiegelman



In [ ]:

    
import numpy as np
import scipy.linalg as la
import matplotlib.pyplot as plt
% matplotlib inline

Fun with Fourier Series

GOAL: visualize some basic behavior of Fourier Series

Fourier Sine Series of 1

Here we will calculate the Fourier Sine series on the interval $[0,1]$ as

$$ 1 = \sum_{n-odd}^\infty \frac{4}{n\pi} \sin(n\pi x) $$



In [ ]:

    
x = np.linspace(0,1,1000)

# small python function to define the partial sum for the truncated Fourier Series f_N
def f_N(x,N):
    na = range(1,N+1,2)
    f_N = np.zeros(x.shape)
    for n in na:
        f_N += 4./(n*np.pi)*np.sin(n*np.pi*x) 
    return f_N

# And make a figure showing f_N for increasing values of N
Na = [ 1, 3, 5, 11, 101]

plt.figure()
for N in Na:
    plt.plot(x,f_N(x,N),label='N={}'.format(N))
plt.plot(x,np.ones(x.shape),'k',label='$f(x)=1$')
plt.xlabel('x')
plt.legend(loc='best')
plt.grid()
plt.show()

Fourier Sine Series of $x(1-x)$

The previous series converges extremely slowly due to the discontinuities at the end-points and that $f(x)=1$ does not share the same boundary conditions as the eigenfunctions $\phi_n(x)=\sin(n\pi x)$ which are zero at the boundaries. For $C^2$ functions that share the same boundary conditions, however, Fourier Series can converge quite quickly. Here we will calculate the Fourier Sine series of the parabola $f(x) = x(1-x)$ on the interval $[0,1]$ which satisfies these conditions.

The Fourier coefficients for this function are

$$ a_n = 2\int_0^1 (x - x^2) sin(n\pi x) dx = \frac{8}{n^3\pi^3} $$

for $n$ odd. These can be found relatively easily by successive integration by parts (Homework).

So the Fourier Sine series of $f$ is $$ x(1-x) = \sum_{n-odd}^\infty \frac{8}{(n\pi)^3} \sin(n\pi x) $$



In [ ]:

    
x = np.linspace(0,1,1000)

# small python function to define the partial sum for the truncated Fourier Series f_N
def f_N(x,N):
    na = range(1,N+1,2)
    f_N = np.zeros(x.shape)
    for n in na:
        f_N += 8./(n*np.pi)**3*np.sin(n*np.pi*x) 
    return f_N

# And make a figure showing f_N for increasing values of N
Na = [ 1, 3, 5, 11, 101]

plt.figure()
for N in Na:
    plt.plot(x,f_N(x,N),label='N={}'.format(N))
plt.plot(x,x*(1-x),'k',label='$f(x)=x(1-x)$')
plt.xlabel('x')
plt.legend(loc='best')
plt.grid()
plt.show()