# Lesson 25. Fourier Series
# --------------------------
#
# Let f be a function defined on the interval [-Pi, Pi]. We often consider it extended to a
# periodic function on the whole real line. The Fourier series of f is the series
> a[0]/2 + Sum(a[k]*cos(k*t)+b[k]*sin(k*t), k=1..infinity);
/infinity \
| ----- |
| \ |
1/2 a[0] + | ) (a[k] cos(k t) + b[k] sin(k t))|
| / |
| ----- |
\ k = 1 /
# where
> a[k]=1/Pi*Int(f(t)*cos(k*t), t=-Pi..Pi);
Pi
/
|
| f(t) cos(k t) dt
|
/
- Pi
a[k] = ----------------------
Pi
> b[k]=1/Pi*Int(f(t)*sin(k*t), t=-Pi..Pi);
Pi
/
|
| f(t) sin(k t) dt
|
/
- Pi
b[k] = ----------------------
Pi
# If f is piecewise continuous and piecewise continuously differentiable, then the Fourier
# series converges to f(t) wherever f is continuous and to (f(t+)+f(t-))/2 where it isn't.
# Let's try an example.
> f:= x -> x^3 + Pi*x^2 - Pi^2*x;
3 2 2
f := x -> x + Pi x - Pi x
# I chose this so f(Pi) = f(-Pi). How do we make f periodic? The following function will
# help.
> saw:= x -> x - 2*Pi*floor(x/(2*Pi)+1/2);
x
saw := x -> x - 2 Pi floor(1/2 ---- + 1/2)
Pi
> plot(saw(x), x=-10 .. 10, discont=true);
> fp:= f @ saw;
fp := f@saw
> plot({f(x), fp(x)}, x=-10 .. 10, y = -10 .. 40);
> a:= unapply( 1/Pi*int(f(t)*cos(k*t), t=-Pi..Pi), k);
a := k -> (
3 3 2 2
k sin(k Pi) Pi + 4 k cos(k Pi) Pi - 8 k sin(k Pi) Pi - 6 cos(k Pi)
----------------------------------------------------------------------
4
k
3 3
k sin(k Pi) Pi + 6 cos(k Pi) + 4 k sin(k Pi) Pi
+ -------------------------------------------------)/Pi
4
k
> simplify(a(k));
2 2
k sin(k Pi) Pi + 2 k cos(k Pi) Pi - 2 sin(k Pi)
2 -------------------------------------------------
3
k
> a(0);
Error, (in a) division by zero
> a(0):= 1/Pi*int(f(t), t=-Pi .. Pi);
3
a(0) := 2/3 Pi
> a(3);
- 4/9 Pi
> assume(k,integer);
> sin(k*Pi);
sin(k~ Pi)
> k:= 'k';
k := k
> sin(k*Pi):= 0; cos(k*Pi):= (-1)^k;
sin(k Pi) := 0
k
cos(k Pi) := (-1)
> sin(j*Pi);
sin(j Pi)
> sin((k+1)*Pi);
sin((k + 1) Pi)
> expand(");
0
> a(k);
2 k 2 k k
4 k (-1) Pi - 6 (-1) (-1)
------------------------ + 6 -----
4 4
k k
----------------------------------
Pi
> normal(");
k
(-1) Pi
4 --------
2
k
> a:= unapply(", k);
k
(-1) Pi
a := k -> 4 --------
2
k
> a(0);
Error, (in a) division by zero
> a(0):= 1/Pi*int(f(t), t=-Pi .. Pi);
3
a(0) := 2/3 Pi
> 1/Pi*int(f(t)*sin(k*t), t=-Pi .. Pi);
k 2 2 k 2 2
(-1) Pi (Pi k - 8) (-1) Pi (Pi k + 4)
- --------------------- + ---------------------
3 3
k k
-----------------------------------------------
Pi
> normal(");
k
(-1)
12 -----
3
k
> b:= unapply(", k);
k
(-1)
b := k -> 12 -----
3
k
>
> a(0)/2 + sum(a(k)*cos(k*t)+b(k)*sin(k*t), k=1..infinity);
/infinity \
| ----- / k k \|
3 | \ | (-1) Pi cos(k t) (-1) sin(k t)||
1/3 Pi + | ) |4 ----------------- + 12 --------------||
| / | 2 3 ||
| ----- \ k k /|
\ k = 1 /
> PSum:= n -> a(0)/2 + sum(a(k)*cos(k*t)+b(k)*sin(k*t), k=1..n);
/ n \
|----- |
| \ |
PSum := n -> 1/2 a(0) + | ) (a(k) cos(k t) + b(k) sin(k t))|
| / |
|----- |
\k = 1 /
> PSum(5);
3
1/3 Pi - 4 Pi cos(t) - 12 sin(t) + Pi cos(2 t) + 3/2 sin(2 t) - 4/9 Pi cos(3 t)
- 4/9 sin(3 t) + 1/4 Pi cos(4 t) + 3/16 sin(4 t) - 4/25 Pi cos(5 t)
12
- --- sin(5 t)
125
> plot(PSum(5), t=-10 .. 10);
> plot({fp(t), PSum(5)}, t=-Pi .. Pi);
> with(plots,display);
[display]
> display( [ seq( plot({fp(t), PSum(nn)}, t=-Pi .. Pi),
> nn = 1 .. 20) ], insequence=true);
> display( [ seq( plot(fp(t) - PSum(nn), t= Pi-1 .. Pi+1), nn=1 .. 20)], insequence=true);
> subs(t=Pi, PSum(n));
/ n \
|----- / k k \|
3 | \ | (-1) Pi cos(k Pi) (-1) sin(k Pi)||
1/3 Pi + | ) |4 ------------------ + 12 ---------------||
| / | 2 3 ||
|----- \ k k /|
\k = 1 /
> ";
3
Pi - 4 Pi Psi(1, n + 1)
> asympt(", n);
3 Pi Pi Pi Pi 1
Pi - 4 ---- + 2 ---- - 2/3 ---- + 2/15 ---- + O(----)
n 2 3 5 6
n n n n
> fp(Pi);
3
Pi
> fp(0);
0
> subs(t=0,PSum(n));
/ n \
|----- / k k \|
3 | \ | (-1) Pi cos(0) (-1) sin(0)||
1/3 Pi + | ) |4 --------------- + 12 ------------||
| / | 2 3 ||
|----- \ k k /|
\k = 1 /
> ";
3
1/3 Pi - 4 Pi hypergeom([1, 1, 1], [2, 2], -1)
(n + 1)
(-1) Pi hypergeom([1, n + 1, n + 1], [2 + n, 2 + n], -1)
- 4 ---------------------------------------------------------------
2
(n + 1)
> asympt(",n);
Error, (in asympt) unable to compute series
>
