# 1.26 Lesson 26

```# Lesson 26.  Fourier Series and Gibbs Phenomenon
# ===============================================
#
# Last time we looked at the Fourier series for the
# function f(x) = x^3 + Pi*x^2 - Pi^2*x.  Here are the
> restart; with(plots,display):
> f:= x -> x^3 + Pi*x^2 - Pi^2*x:
> saw:= x -> x - 2*Pi*floor(x/(2*Pi)+1/2):
> fp:= f @ saw:
> sin(k*Pi):= 0: cos(k*Pi):= (-1)^k:
> a:= unapply(normal(1/Pi*int(f(t)*cos(k*t), t=-Pi .. Pi)), k);

k
(-1)  Pi
a := k -> 4 --------
2
k
> a(0):= 1/Pi*int(f(t),t=-Pi..Pi);

3
a(0) := 2/3 Pi
> b:= unapply(normal(1/Pi*int(f(t)*sin(k*t),t=-Pi..Pi)),k);

k
(-1)
b := k -> 12 -----
3
k
# It may be helpful to have two "partial sum" functions,
# one using "sum" that tries to do the sum symbolically,
# and one using "Sum" that does not.
> 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:= subs(sum=Sum,eval(Psum));

/  n                                  \
|-----                                |
| \                                   |
PSum := n -> 1/2 a(0) + |  )   (a(k) cos(k t) + b(k) sin(k t))|
| /                                   |
|-----                                |
\k = 1                                /
# We say that the difference between f and PSum(n) at
# t=Pi was of order 1/n.  What about at t=0?
> eval(subs(t=0,PSum(n)));

/  n             \
|-----       k   |
3   | \      (-1)  Pi|
1/3 Pi  + |  )   4 --------|
| /          2   |
|-----      k    |
\k = 1           /
# This is an alternating series.  The size of the tail
# is therefore at most the absolute value of the next
# term, and thus of order 1/n^2.
#
# The next example is a "sawtooth" function which has
# a jump at t=Pi.
> f:= x -> x;

f := x -> x
> fp(x);

x
x - 2 Pi floor(1/2 ---- + 1/2)
Pi
> a:= unapply(normal(1/Pi*int(f(t)*cos(k*t), t=-Pi .. Pi)), k);

a := 0
# We should have expected this: f is an odd function,
# so there are no cosine terms.
> b:= unapply(normal(1/Pi*int(f(t)*sin(k*t), t=-Pi .. Pi)), k);

k
(-1)
b := k -> - 2 -----
k
> Psum(5);

2 sin(t) - sin(2 t) + 2/3 sin(3 t) - 1/2 sin(4 t) + 2/5 sin(5 t)
> PSum(n);

n
-----         k
\        (-1)  sin(k t)
)   - 2 --------------
/               k
-----
k = 1
> plot({PSum(15),fp(t)},t= 0 .. 2*Pi);

> display([seq(plot({PSum(nn),fp(t)},t=0 .. 2*Pi), nn=1 .. 15)],insequence=true);

# If you look closely, you'll see that the highest wiggle is always some distance above the
# graph of f.
# As n -> infinity, the "overshoot" doesn't disappear.
# This is called the Gibbs phenomenon.
#
# The highest point in the graph of PSum(n) should be a point where the derivative is 0.

> ds:=diff(PSum(n),t);

n
-----
\            k
ds :=   )   - 2 (-1)  cos(k t)
/
-----
k = 1
> ds:=value(");

(n + 1)
(n + 1)                  (cos(t) - 1) (-1)        sin((n + 1) t)
ds := (-1)        cos((n + 1) t) - --------------------------------------- + 1
sin(t)
# One way for ds to be 0 is for cos((n+1) t) to be (-1)^n and sin((n+1) t) = 0.  That will
# happen at
# t = Pi - Pi/(n+1).  Is this the point we want?
> plot(subs(n=10, ds), t=Pi - Pi/11 .. Pi);

> eval(subs(t=Pi-Pi/(n+1), ds));

(n + 1)             /       Pi \
(-1)        cos((n + 1) |Pi - -----|)
\     n + 1/

/        Pi      \     (n + 1)             /       Pi \
|- cos(-----) - 1| (-1)        sin((n + 1) |Pi - -----|)
\      n + 1     /                         \     n + 1/
- -------------------------------------------------------- + 1
Pi
sin(-----)
n + 1
> simplify(");

n                 Pi         n                 Pi         n
- ((-1)  cos(Pi n) sin(-----) + (-1)  sin(Pi n) cos(-----) + (-1)  sin(Pi n)
n + 1                        n + 1

Pi      /       Pi
- sin(-----))  /  sin(-----)
n + 1   /       n + 1
> subs(cos(Pi*n)=(-1)^n, sin(Pi*n)=0,");

n 2       Pi           Pi
((-1) )  sin(-----) - sin(-----)
n + 1        n + 1
- --------------------------------
Pi
sin(-----)
n + 1
> s1:= subs(t=Pi-Pi/m,PSum(m-1));

m - 1         k       /      Pi \
-----     (-1)  sin(k |Pi - ----|)
\                    \       m /
s1 :=   )   - 2 ------------------------
/                    k
-----
k = 1
> simplify(expand(s1));

/m - 1     (2 k)     k Pi \
|----- (-1)      sin(----)|
| \                    m  |
2 |  )   -------------------|
| /             k         |
|-----                    |
\k = 1                    /
> s1:=subs((-1)^(2*k)=1, ");

/m - 1     k Pi \
|----- sin(----)|
| \          m  |
s1 := 2 |  )   ---------|
| /        k    |
|-----          |
\k = 1          /
> limit(s1,m=infinity);

/m - 1     k Pi \
|----- sin(----)|
| \          m  |
limit         2 |  )   ---------|
m -> infinity   | /        k    |
|-----          |
\k = 1          /
> g:= x -> sin(Pi*x)/x;

sin(Pi x)
g := x -> ---------
x
> s1 =  2*Sum(1/m * g(k/m), k=1..m);

/m - 1     k Pi \     /  m       k Pi \
|----- sin(----)|     |----- sin(----)|
| \          m  |     | \          m  |
2 |  )   ---------| = 2 |  )   ---------|
| /        k    |     | /        k    |
|-----          |     |-----          |
\k = 1          /     \k = 1          /
> lims:= 2*int(g(x), x=0..1);

lims := 2 Si(Pi)
> evalf(");

3.703874104
> gibbs:= ("-Pi)/(2*Pi);

3.703874104 - Pi
gibbs := 1/2 ----------------
Pi
> evalf(");

.08948987215
# The next example is an illustration of the use of
# Fourier series in a physical problem: the one-dimensional heat equation.  The
# equation describes the temperature u(x,t) at position x and time t in a uniform bar
# of material with insulated surface:
> heateqn:= diff(u(x,t),x\$2) = diff(u(x,t),t);

2
d               d
heateqn := ----- u(x, t) = ---- u(x, t)
2             dt
dx
# We'll suppose the temperature at the ends of the bar (x=0 and x=Pi) is kept at 0.
> boundary:= { u(0,t) = 0 , u(Pi, t) = 0 };

boundary := {u(0, t) = 0, u(Pi, t) = 0}
# We also specify the initial condition u(x,0) = f(x):
> f:= x -> x*cos(2*x);

f := x -> x cos(2 x)
> plot(f(x), x=0..Pi);
# Note that u(x,0) doesn't satisfy the boundary conditions.  But that's not a problem.
# Now it turns out that there are solutions to the differential equation sin(a x) exp(-b t)
> v:= (x,t) -> sin(a*x)*exp(-b*t);

v := (x,t) -> sin(a x) exp(- b t)
> subs(u=v, heateqn);

2
d               d
----- v(x, t) = ---- v(x, t)
2             dt
dx
> ";

2
- sin(a x) a  exp(- b t) = - sin(a x) b exp(- b t)
# So this works if b = a^2.  It would also work with "cos" instead of "sin".  But for our
# boundary condition u(0,t)=0 we need "sin", and for u(Pi,t)=0 we need "a" to be an integer.
# A linear combination of solutions is also a solution.  So we could look for a solution of the
# form
> Sum(c[k]*sin(k*x)*exp(-k^2*t), k=1..infinity);

infinity
-----
\                           2
)     c[k] sin(k x) exp(- k  t)
/
-----
k = 1
# At t=0 we should have
> f(x) = eval(subs(t=0, "));

infinity
-----
\
x cos(2 x) =    )     c[k] sin(k x)
/
-----
k = 1
# So we want to expand f(x) for 0 <= x <= Pi in a kind of Fourier series, but just involving
# sin and not cos.  Well, as we've seen, for a Fourier series on [-Pi..Pi] the cos terms will be
# 0 if the function is odd.   It happens that our x*cos(2*x) is odd: if it were not, we'd redefine
# it on [-Pi .. 0] to make it so.
> b:= unapply(normal(1/Pi*int(f(x)*sin(k*x), x=-Pi .. Pi)), k);

k
(-1)  k
b := k -> - 2 -----------------
(2 + k) (- 2 + k)
# Hmm, k=2 would seem to be a problem here.
> b(2):= 1/Pi*int(f(x)*sin(2*x), x=-Pi .. Pi);

b(2) := -1/4
# For an approximate solution, we'll take a partial sum of the Fourier series.  I'll use 5
# terms (note the terms for higher n are damped down very rapidly).
> sum(b(k)*sin(k*x)*exp(-k^2*t), k=1..5);
Error, (in sum) division by zero

> Sum(b(k)*sin(k*x)*exp(-k^2*t), k=1..5);

5
-----         k                   2
\        (-1)  k sin(k x) exp(- k  t)
)   - 2 ----------------------------
/              (2 + k) (- 2 + k)
-----
k = 1
# There's a problem here: in the sum (or the Sum), the b(k) is evaluated first with k a
# symbolic variable, so the formula is used.  This is what's known as "premature
# evaluation".  We need to ensure that b(k) is not evaluated before numbers are substituted
# for k, so Maple will use our special value for b(2).  This can be ensured
# by putting quotes around 'b(k)'.
> sum('b(k)'*sin(k*x)*exp(-k^2*t), k=1..5);

- 2/3 sin(x) exp(- t) - 1/4 sin(2 x) exp(- 4 t) + 6/5 sin(3 x) exp(- 9 t)

10
- 2/3 sin(4 x) exp(- 16 t) + ---- sin(5 x) exp(- 25 t)
21
> u:= unapply(",(x,t));

u := (x,t) -> - 2/3 sin(x) exp(- t) - 1/4 sin(2 x) exp(- 4 t)

+ 6/5 sin(3 x) exp(- 9 t) - 2/3 sin(4 x) exp(- 16 t)

+ 10/21 sin(5 x) exp(- 25 t)
> with(plots,animate):
> animate(u(x,t), x=0..Pi, t = 0 .. 1);
>

```