Double click here and enter the student numbers for all minigroup members. Don't enter any names.

1. First student number
2. Second student number
3. Third student number.

Before you start work on the project, click on this link to read the MATH0011 project instructions.

Project 6 - The Gibbs Phenomenon

This project is about a strange occurence in the theory of Fourier series. You don't need to have any previous knowledge of Fourier series (which are covered in MATH0016 Methods 3) to complete this project.

In your analysis modules you looked at power series representations for functions - for example, you can represent the function $f(x) = e^x$ by a series $$ \sum_{i=0}^\infty \frac{x^i}{i!}.$$ A power series for a function $f$ expresses $f(x)$ as a sum of multiples of powers of $x$. Fourier series are similar, but they represent functions as sums of multiples of functions of the form $\sin(nx)$ and $\cos(mx)$.

The first cell below imports some functions that you will need for plotting.

import math
import matplotlib.pyplot as plt
import numpy as np

Part 1

Write a function partial_sum(x, n) whose output is

$$ \sum_{i=0}^n \frac{\sin((2i+1)x)}{2i+1}. $$

One way to do this is with a list comprehension producing

$$[ \sin(x), \sin(3x)/3, \sin(5x)/5, \dots ] $$

together with the sum function, but you can use any method you want.

def partial_sum(x,n):
    s=[math.sin((2*i+1)*x)/(2*i+1) for i in range(n+1)]
    return sum(s)

Your partial_sum function is an example of a partial sum of a Fourier series. Specifically, this is the Fourier series representation of the square wave function which you will plot in the next exercise.

Part 2

The function square_wave(x) is defined by $$ \mathtt{square\_wave}(x) = \begin{cases} -\pi/4 & x < 0 \\ \pi/4 & 0 \leq x \leq \pi \\ -\pi/4 & \pi < x \leq 2\pi \\ \pi/4 & x > 2\pi\end{cases}$$ Plot the graph of square_wave(x) for $-0.2 \leq x \leq 2\pi + 0.2$ using plt.plot. On the same axes plot partial_sum(x, n) for n = 3, and n=15.

#plot square_wave(x)
plt.plot([-0.2,0,0,math.pi,math.pi,2*math.pi,2*math.pi,2*math.pi+0.2],[-math.pi/4,-math.pi/4,math.pi/4,math.pi/4,-math.pi/4,-math.pi/4,math.pi/4,math.pi/4],'b',label="square_wave(x)")

#plot partial_sum(x,3) and partial_sum(x,15)
xs=np.linspace(-0.2,2*math.pi+0.2,500)
ys=[partial_sum(x,3) for x in xs]
zs=[partial_sum(x,15) for x in xs]
plt.plot(xs,ys,'r',label="partial_sum(x,n) for n=3")
plt.plot(xs,zs,'g',label="partial_sum(x,n) fro x=15")

#label
plt.xlabel("x-axis")
plt.ylabel("y-axis")
plt.grid()
plt.legend()

<matplotlib.legend.Legend at 0x7fe1cbdcb760>

You are familiar with the fact that power series don't always converge - that's why you spend so much time in Analysis 1 thinking about radii of convergence. Fourier series also have interesting convergence behaviour. Your plots should show a special behaviour of Fourier series called the Gibbs phenomenon, where the partial sums are larger than the value of square_wave(x) for values of $x$ near to where square_wave(x) jumps from positive to negative. In the next exercise you will investigate exactly how large is the difference between the partial sums and the square wave.

Part 3

Compute an approximation to $I = \int_0 ^\pi \operatorname{sinc}(t) \mathrm{d}t$, where

$$ \operatorname{sinc}(x) = \begin{cases} 1 & x = 0 \\ \frac{\sin(x)}{x} & x \neq 0 \end{cases}.$$

You can either do this by using the code for the rectangle rule from our weekly notebooks, or by importing a command from one of the modules included with Anaconda - try searching the online documentation for numpy and scipy.

Then, on the same axes, plot partial_sum(x, 100), square_wave(x), and the lines $y = \pm \frac I2$ for $-0.2 \leq x \leq 2\pi + 0.2$. You should see that the straight lines $y = \pm \frac I2$ show exactly how much the partial sums of the Fourier series differ from the square wave near to the points where the square wave changes value.

#the code for the rectangle rule from our weekly notebooks
def rectangle_rule(f, a, b, N):
    total = 0
    for i in range(N):
        total = total + f(a + (b-a)*i / N)
    return ((b-a) / N) * total

#define a function sinc(x)
def sinc(x):
    if x==0:
        return 1
    return math.sin(x)/x

#compute an approximation to I
I=rectangle_rule(sinc,0,math.pi,100000)
I_str=str(rectangle_rule(sinc,0,math.pi,100000))
print("I="+I_str)

#plot y=+-I/2
plt.plot([-0.2,2*math.pi+0.2],[I/2,I/2],'r',label="y=I/2")
plt.plot([-0.2,2*math.pi+0.2],[-I/2,-I/2],'r',label="y=-I/2")

#plot square_wave(x)
plt.plot([-0.2,0,0,math.pi,math.pi,2*math.pi,2*math.pi,2*math.pi+0.2],[-math.pi/4,-math.pi/4,math.pi/4,math.pi/4,-math.pi/4,-math.pi/4,math.pi/4,math.pi/4],'b',label="square_wave(x)")

#plot partial_sum(x,100)
xs=np.linspace(-0.2,2*math.pi+0.2,500)
ws=[partial_sum(x,100) for x in xs]
plt.plot(xs,ws,'g',label="partial_sum(x,n) for n=100")

#label
plt.xlabel("x-axis")
plt.ylabel("y-axis")
plt.grid()
plt.legend()

I=1.8519527599195542

<matplotlib.legend.Legend at 0x7fe1c9b3da30>

Part 4

Multiplying the $i$th term in the partial sum

$$\sum_{i=0}^n \frac{\sin((2i+1)x)}{2i+1}$$

by the Lanczos sigma factor $\operatorname{sinc}\left(\frac{(2i+1)\pi}{4n}\right)$ gives a new function

$$\mathtt{lanczos\_partial\_sum}(x,n) = \sum_{i=0}^n \operatorname{sinc}\left(\frac{(2i+1)\pi}{4n}\right) \frac{\sin((2i+1)x)}{2i+1}.$$

On the same axes, plot lanczos_partial_sum(x, 100), square_wave(x), and the lines $y = \pm \frac I2$. You will see that the Gibbs phenomenon is reduced, but not eliminated.

def lanczos_partial_sum(x,n):
    S=[(sinc((2*i+1)*math.pi/(4*n)))*(math.sin((2*i+1)*x)/(2*i+1)) for i in range(n+1)]
    return sum(S)

#plot y=+-I/2
plt.plot([-0.2,2*math.pi+0.2],[I/2,I/2],'r',label="y=I/2")
plt.plot([-0.2,2*math.pi+0.2],[-I/2,-I/2],'r',label="y=-I/2")

#plot square_wave(x)
plt.plot([-0.2,0,0,math.pi,math.pi,2*math.pi,2*math.pi,2*math.pi+0.2],[-math.pi/4,-math.pi/4,math.pi/4,math.pi/4,-math.pi/4,-math.pi/4,math.pi/4,math.pi/4],'b',label="square_wave(x)")

#plot lanczos_partial_sum(x,100)
xs=np.linspace(-0.2,2*math.pi+0.2,500)
ls=[lanczos_partial_sum(x,100) for x in xs]
plt.plot(xs,ls,'g',label="lanczos_partial_sum(x,n) for n=100")

#label
plt.xlabel("x-axis")
plt.ylabel("y-axis")
plt.grid()
plt.legend()

<matplotlib.legend.Legend at 0x7fe1c9b3d6a0>

Submitting your project

Make sure you have done all of the following things.

0. Included all minigroup members' student numbers at the top of this notebook.
1. Read through every part of the project to check you have answered all of it.
2. Carefully read and followed all of the MATH0011 project instructions.
3. Checked that all of your code works correctly.

If you have, you're ready to submit. One of the minigroup members should download the completed notebook (in CoCalc, click the File menu next to the green Save button, then click Download) and submit it on the MATH0011 Moodle. Please submit one .ipynb file per minigroup.