Fourier Series

You can use the following commands to calculate the $n\mbox{th}$ partial sum of the Fourier series of the expression $f$ on the interval $[-L,L]$ .

syms x k L n
evalin(symengine,'assume(k,Type::Integer)');
a = @(f,x,k,L) int(f*cos(k*sym('pi')*x/L)/L,x,-L,L);
b = @(f,x,k,L) int(f*sin(k*sym('pi')*x/L)/L,x,-L,L);
fs = @(f,x,n,L) a(f,x,0,L)/2 + ...
    symsum(a(f,x,k,L)*cos(k*pi*x/L) + b(f,x,k,L)*sin(k*pi*x/L),k,1,n);
pfs = @(f,x,n,L) pretty(simple(fs(f,x,n,L)));

I'll do it for the function $f(x)=|x|$ on the interval $[-1,1]$ .

f = abs(x)

 
f =
 
abs(x)

The $n\mbox{th}$ partial sum is

fs(f,x,n,1)

 
ans =
 
1/2 - (2*sum(((exp(-pi*k*x*i)/2 + exp(pi*k*x*i)/2)*(2*((exp(-(pi*k*i)/2)*i)/2 - (exp((pi*k*i)/2)*i)/2)^2 - pi*k*((exp(-pi*k*i)*i)/2 - (exp(pi*k*i)*i)/2)))/k^2, k == 1..n))/pi^2

Here's a more readable form.

pfs(f,x,n,1)

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

The Fourier series is

pfs(f,x,inf,1)

           Inf              /      / pi k \2                  \
           ---  cos(pi k x) | 2 sin| ---- |  - pi k sin(pi k) |
           \                \      \  2   /                   /
        2  /    -----------------------------------------------
           ---                         2
          k = 1                       k
  1/2 - -------------------------------------------------------
                                    2
                                  pi

Here are the plots of the partial sums for $n=2,5,10$ . The plot also shows the function $f$ .

ezplot(fs(f,x,2,1),-1,1)
hold on
ezplot(f,-1,1)
hold off
title('Partial sum with n=2')

ezplot(fs(f,x,5,1),-1,1)
hold on
ezplot(f,-1,1)
hold off
title('Partial sum with n=5')

ezplot(fs(f,x,10,1),-1,1)
hold on
ezplot(f,-1,1)
hold off
title('Partial sum with n=10')

Now I do it with the function $f(x) = x$ on $[-1,1]$ .

f = x

fs(f,x,n,1)

 
f =
 
x
 
 
ans =
 
-(2*sum((((exp(-pi*k*x*i)*i)/2 - (exp(pi*k*x*i)*i)/2)*(- (exp(-pi*k*i)*i)/2 + (exp(pi*k*i)*i)/2 + pi*k*(exp(-pi*k*i)/2 + exp(pi*k*i)/2)))/k^2, k == 1..n))/pi^2

pfs(f,x,n,1)

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

Here are plots of the partial sums for $n=2,5,10,20,50$ .

ezplot(fs(f,x,2,1),-1,1)
hold on
ezplot(f,-1,1)
hold off
title('Partial sum with n=2')

ezplot(fs(f,x,5,1),-1,1)
hold on
ezplot(f,-1,1)
hold off
title('Partial sum with n=5')

ezplot(fs(f,x,10,1),-1,1)
hold on
ezplot(f,-1,1)
hold off
title('Partial sum with n=10')

ezplot(fs(f,x,20,1),-1,1)
hold on
ezplot(f,-1,1)
hold off
title('Partial sum with n=20')

ezplot(fs(f,x,50,1),-1,1)
hold on
ezplot(f,-1,1)
hold off
title('Partial sum with n=50')

If you zoom in here, you notice that the graph of the partial sum of the Fourier series looks very jagged. That is because ezplot does not plot enough points compared to the number of oscillations in the functions in the partial sum. We can fix this problem using plot. This also allows us to use a different colors for the plot of the original function and the plot of the partial sum. To use plot, we need to turn the partial sum into an inline vectorized function and specify the points where it will be evaluated.

g = inline(vectorize(fs(f,x,50,1)));
X = -1:.001:1;
plot(X,g(X),'r')
hold on
ezplot(f,-1,1)
hold off
title('Partial sum with n=50')

g = inline(vectorize(fs(f,x,100,1)));
X = -1:.0005:1;
plot(X,g(X),'r')
hold on
ezplot(f,-1,1)
hold off
title('Partial sum with n=100')

Notice that no matter how many terms we take, the partial sum always overshoots the function $x$ at the end points of the interval. If we do our two plots on the interval $[-2,2]$ , we see that outside $[-1,1]$ , the partial sum doesn't have much to do with the function.

ezplot(fs(f,x,20,1),-2,2)
hold on
ezplot(f,-2,2)
hold off

A Fourier series on $[-L,L]$ is $2L$ periodic, and so are all its partial sums. So, what we are really doing when we compute the Fourier series of a function $f$ on the interval $(-L,L)$ is computing the Fourier series of a $2L$ periodic extension of f. To do that in MATLAB, we have to make use of the unit step function u(x), which is 0 if $x < 0$ and 1 if $x \ge 0$ . It is known as the Heaviside function, and the MATLAB command for it is heaviside. The function $h(x) = u(x-a)u(b-x)$ is 1 on the interval [a,b] and 0 outside the interval. So the formula

f = heaviside(x+3)*heaviside(-1-x)*(x+2) + heaviside(x+1)*heaviside(1-x)*x ...
    + heaviside(x-1)*heaviside(3-x)*(x-2);

extends $f(x) = x$ to be periodic on $[-3,3]$ , with period 2. To check that we've extended it correctly, we plot it.

ezplot(f,-3,3)

Here is a plot of the function and its Fourier series.

ezplot(fs(f,x,20,1),-3,3)
hold on
ezplot(f,-3,3)
hold off

It is somewhat clearer if the plots aren't the same color, so I'll use plot for the partial sum.

X = -3:0.001:3;
g = inline(vectorize(fs(f,x,20,1)));

ezplot(f,-3,3)
hold on
plot(X,g(X),'r')
hold off

The overshoots are at the discontinuities of the 2 periodic extension of f. This is called the Gibbs Phenomenon.