#
#
#
# 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
--------------------------------------------------------------------------------
