# 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
# definitions we made.
> 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);
>
