# 2.11 Fourier series cont.

```#
#
#
#                                               WORKSHEET # 22
#
#
#                                                Fourier series cont
#
# .
> restart;
--------------------------------------------------------------------------------
# In this worksheet we examine the Fourier series for the function
# which equals x on the interval [-Pi,Pi) and has period 2*Pi.
--------------------------------------------------------------------------------
> f:=x->x;

f := x -> x
--------------------------------------------------------------------------------
> ?floor
--------------------------------------------------------------------------------
# floor(x) is the greatest integer less than or equal to x.
--------------------------------------------------------------------------------
> fper:=x->f(x-2*Pi*floor((x/(2*Pi))+1/2));

x
fper := x -> f(x - 2 Pi floor(1/2 ---- + 1/2))
Pi
--------------------------------------------------------------------------------
> with(plots):
--------------------------------------------------------------------------------
> Fper:=plot(fper(x),x=-3*Pi..3*Pi):
--------------------------------------------------------------------------------
> with(plots):
--------------------------------------------------------------------------------
> display(Fper,title= `the periodic extension of x`);
--------------------------------------------------------------------------------
#

** Maple V Graphics **

--------------------------------------------------------------------------------
# Notice that the vertical lines are not really part of the graph.  The
# function fper(x) has a jump discontinuity at all odd multiples of Pi.
--------------------------------------------------------------------------------
> fper(Pi);

- Pi
--------------------------------------------------------------------------------
> limit(fper(x),x=Pi,left);fper(Pi);limit(fper(x),x=Pi,right);

Pi

- Pi

- Pi
--------------------------------------------------------------------------------
> for k from -5 to 5 do print(k,limit(fper(x),x=k*Pi,left),fper(k*Pi),limit(fper(x),x=k*Pi,right)) od;k:='k':

-5, Pi, - Pi, - Pi

-4, 0, 0, 0

-3, Pi, - Pi, - Pi

-2, 0, 0, 0

-1, Pi, - Pi, - Pi

0, 0, 0, 0

1, Pi, - Pi, - Pi

2, 0, 0, 0

3, Pi, - Pi, - Pi

4, 0, 0, 0

5, Pi, - Pi, - Pi
--------------------------------------------------------------------------------
# Recall that for a function f(x) defined on the interval [0,2*Pi] the
# Fourier series is a linear combination of the periodic functions
# cos(n*x), sin(n*x) and a constant term:
#
#
# f(x)=(1/2)*a[0]+Sum(a[n]*cos(n*x)+b[n]*sin(n*x),n=1..infinity);
#
# The coefficients in this formula are given by the following Fourier
# integral formulas.  The formula for a[n] is valid for n=0,1,... and that
# for b[n] is valid for n=1,2,...
#
#         a[n]=(1/Pi)*Int(f(x)*cos(n*x),x=0..2*Pi);
#         b[n]=(1/Pi)*Int(f(x)*sin(n*x),x=0..2*Pi)
#
--------------------------------------------------------------------------------
# For functions that are periodic we can choose any interval of length
# 2*Pi for the Fourier integrals.  Choosing the interval [-Pi,Pi] gives:
--------------------------------------------------------------------------------
> A:=int(x*cos(n*x),x=-Pi..Pi);

A := 0
--------------------------------------------------------------------------------
# Thus all the cosine coefficients are zero.  This should have been
# obvious since the periodic extension of f(x)=x is an odd function.
--------------------------------------------------------------------------------
> B:=int(x*sin(n*x),x=-Pi..Pi);

- sin(n Pi) + n cos(n Pi) Pi
B := - 2 ----------------------------
2
n
--------------------------------------------------------------------------------
> weknow:=sin(n*Pi)=0,cos(n*Pi)=(-1)^n;

n
weknow := sin(n Pi) = 0, cos(n Pi) = (-1)
--------------------------------------------------------------------------------
> subs(weknow,B);

n
(-1)  Pi
- 2 --------
n
--------------------------------------------------------------------------------
> b:='b':
--------------------------------------------------------------------------------
> b:=n->-2*(-1)^n/n;

n
(-1)
b := n -> - 2 -----
n
--------------------------------------------------------------------------------
> S:=n->sum(b(k)*sin(k*x),k=1..n);

n
-----
\
S := n ->   )   b(k) sin(k x)
/
-----
k = 1
--------------------------------------------------------------------------------
# Next we compare the the graph of f(x)=x with the partial sums of the
# Fourier series.
--------------------------------------------------------------------------------
> for n from 1 to 30 do P.n:=plot(S(n),x=-3*Pi..3*Pi) od:n:='n':
--------------------------------------------------------------------------------
> with(plots):
--------------------------------------------------------------------------------
> display([Fper,P1,P2,P3,P4,P5],title=`first 5 partial sums for f(x)=x` );
--------------------------------------------------------------------------------
#

** Maple V Graphics **

--------------------------------------------------------------------------------
#
> display([Fper,P1,P5],title=`1st and 5th partial sums for  f(x)=x` );\
display([Fper,P10],title=`10th partial sum for f(x)=x` );\
display([Fper,P20 ],title=` 20th partial sum for f(x)=x` );\
display([Fper,P20 ],title=` 30th partial sum for f(x)=x` );
--------------------------------------------------------------------------------
#

** Maple V Graphics **

--------------------------------------------------------------------------------
#

** Maple V Graphics **

--------------------------------------------------------------------------------
#

** Maple V Graphics **

--------------------------------------------------------------------------------
#

** Maple V Graphics **

--------------------------------------------------------------------------------
# These graphs nicely approximate the saw tooth wave fper(x) except at
# odd multiples of Pi; but this is to be expected since fper(x) has jump
# discontinuities at the odd multiples of Pi.  Also notice that the nth
# partial sum S(n) seems to overshoot the saw tooth wave as it
# approaches the jump discontinuities from the left.  In fact the nth
# partial sum does overshoot the graph at x=Pi by about 18%.  This is
# known as Gibb's phenomenum.  The overshoot actually occurs at the
# local maximum of S(n).  This we find by solving Diff(S(n))=0.
--------------------------------------------------------------------------------
> S(n);

n
-----         k
\        (-1)  sin(k x)
)   - 2 --------------
/               k
-----
k = 1
--------------------------------------------------------------------------------
> Diff(S(n),x)=sum(-2*(-1)^k*cos(k*x),k=1..n);

/  n                     \
|-----         k         |              (n + 1)
d  | \        (-1)  sin(k x)|   sin(x) (-1)        sin((n + 1) x)
---- |  )   - 2 --------------| = ---------------------------------
dx  | /               k      |               1 + cos(x)
|-----                   |
\k = 1                   /

2
(n + 1)                    sin(x)
+ (-1)        cos((n + 1) x) + ---------- + cos(x)
1 + cos(x)
--------------------------------------------------------------------------------
# This is correct.  In other words the following  is an identity.
--------------------------------------------------------------------------------
> Sum(-2*(-1)^k*cos(k*x),k=1..n)=sum(-2*(-1)^k*cos(k*x),k=1..n);

n
-----                                 (n + 1)
\            k            sin(x) (-1)        sin((n + 1) x)
)   - 2 (-1)  cos(k x) = ---------------------------------
/                                     1 + cos(x)
-----
k = 1

2
(n + 1)                    sin(x)
+ (-1)        cos((n + 1) x) + ---------- + cos(x)
1 + cos(x)
--------------------------------------------------------------------------------
> exp(I*x)=cos(x)+I*sin(x);exp(-I*x)=cos(x)-I*sin(x);

exp(I x) = cos(x) + I sin(x)

exp(- I x) = cos(x) - I sin(x)
--------------------------------------------------------------------------------
# Therefore
--------------------------------------------------------------------------------
> cos(x)=(exp(I*x)+exp(-I*x))/2;sin(x)=(exp(I*x)-exp(-I*x))/2*I;

cos(x) = 1/2 exp(I x) + 1/2 exp(- I x)

sin(x) = 1/2 I (exp(I x) - exp(- I x))
--------------------------------------------------------------------------------
# Substituting in cos(kx)=1/2(exp(kIx)+exp(-kIx)), summing the
# geometric series involved and then simplifying leads to the above
# identity.
--------------------------------------------------------------------------------
# Now we return to the problem of finding where S(n) achieves its
# maximum value.
--------------------------------------------------------------------------------
> Diff(S(n),x)=sum(-2*(-1)^k*cos(k*x),k=1..n);

/  n                     \
|-----         k         |              (n + 1)
d  | \        (-1)  sin(k x)|   sin(x) (-1)        sin((n + 1) x)
---- |  )   - 2 --------------| = ---------------------------------
dx  | /               k      |               1 + cos(x)
|-----                   |
\k = 1                   /

2
(n + 1)                    sin(x)
+ (-1)        cos((n + 1) x) + ---------- + cos(x)
1 + cos(x)
--------------------------------------------------------------------------------
# Near x=Pi the last 3 terms in this expression are approximately equal
# to 2, so we should be considering the function which is  just the sum
# of the first term and 2.  Multiplying the first term, top and bottom, by
# 1-cos(x) and then simplifying leads to the function h(x) defined as
# follows.
--------------------------------------------------------------------------------
> h(x):=(-1)^(n+1)*(1-cos(x))*sin((n+1)*x)/sin(x)+2;

(n + 1)
(-1)        (1 - cos(x)) sin((n + 1) x)
h(x) := --------------------------------------- + 2
sin(x)
--------------------------------------------------------------------------------
# The critical point of S(n)  just  to the left of Pi will approximate the
# solution of h(x)=0 near Pi.  We can calculate the limit of h(x) as x
# tends to Pi by using L'Hospital's rule.
--------------------------------------------------------------------------------
> Limit(h(x),x=Pi)=-2*n;

(n + 1)
(-1)        (1 - cos(x)) sin((n + 1) x)
Limit   --------------------------------------- + 2 = - 2 n
x -> Pi                  sin(x)
--------------------------------------------------------------------------------
> weknow:=sin(n*Pi)=0;

weknow := sin(n Pi) = 0
--------------------------------------------------------------------------------
# Next we calculate h(x) at x=n*Pi/(n+1).
--------------------------------------------------------------------------------
> subs(weknow,subs(x=n*Pi/(n+1),h(x)));

2
--------------------------------------------------------------------------------
# In other words in the range n*Pi/(n+1)<x<Pi, h(x) has gone from the
# small positive number 2 down to the large negative number -2*n.
# There is a solution of h(x)=0 in this range, and it should be close to
# x=n*Pi/(n+1).  Therefore the maximum value of S(n) would appear to
# occur close to x=n*Pi/(n+1).  In any case we will test the amount of
# overshooting by calculating S(n) at x=n*Pi/(n+1).
--------------------------------------------------------------------------------
> G:=n->evalf(subs(x=(n*Pi)/(n+1),S(n)));

n Pi
G := n -> evalf(subs(x = -----, S(n)))
n + 1
--------------------------------------------------------------------------------
#
> for n from 1 to 10 do print(n,G(n),evalf(G(n)/Pi)) od;n:='n':

1, 2., .6366197722

2, 2.598076210, .8269933425

3, 2.885618083, .9185207633

4, 3.054557324, .9722957939

5, 3.165704771, 1.007675125

6, 3.244375368, 1.032716754

7, 3.302985531, 1.051372948

8, 3.348338914, 1.065809378

9, 3.384475470, 1.077312001

10, 3.413945202, 1.086692508
--------------------------------------------------------------------------------
> for n from 10 by 10 to 100 do print(n,G(n),evalf(G(n)/Pi)) od;n:='n':

10, 3.413945202, 1.086692508

20, 3.553086981, 1.130982712

30, 3.601987524, 1.146548238

40, 3.626938404, 1.154490350

50, 3.642072922, 1.159307817

60, 3.652231869, 1.162541510

70, 3.659522442, 1.164862172

80, 3.665009202, 1.166608662

90, 3.669287874, 1.167970605

100, 3.672717881, 1.169062409
--------------------------------------------------------------------------------
> for n from 100 by 10 to 200 do print(n,G(n),evalf(G(n)/Pi)) od;n:='n':

100, 3.672717878, 1.169062409

110, 3.675528964, 1.169957206

120, 3.677874768, 1.170703898

130, 3.679861972, 1.171336445

140, 3.681566993, 1.171879170

150, 3.683045873, 1.172349912

160, 3.684340875, 1.172762124

170, 3.685484303, 1.173126089

180, 3.686501263, 1.173449797

190, 3.687411610, 1.173739570

200, 3.688231329, 1.174000494
--------------------------------------------------------------------------------
# The third column in these calculations is suggesting that the
# overshoot is approximately 17%.  Actually it is closer to 18%.
--------------------------------------------------------------------------------
> evalf(G(1000)/Pi);

1.177980559
--------------------------------------------------------------------------------
# In fact it can be shown that the limit of G(n)/Pi as n tends to infinity
# exists and equals 1.178979744.  the limit is exactly given by the
# following formula.
--------------------------------------------------------------------------------
> Limit(G(n),n=infinity)=(2/Pi)*Int(sin(x)/x,x=0..Pi);

Pi
/
|   sin(x)
|   ------ dx
n           k     k n Pi       |      x
-----     (-1)  sin(------)     /
\                   n + 1      0
Limit           )   - 2 ----------------- = 2 --------------
n -> infinity  /                k                   Pi
-----
k = 1
--------------------------------------------------------------------------------
> evalf((2/Pi)*Int(sin(x)/x,x=0..Pi));\

1.178979744
--------------------------------------------------------------------------------

```