Calculating a Fourier series with MATLAB manually problem

Hi there, I'm trying to transfer this academic terms learnings on the fourier series to MATLAB, and running into difficulties.

I want to be able to plot a fourier series approximation against the original function but when I try and calculate the fourier series to 5 steps of n in the final line I get a division by zero error.

I can't see where this division by zero occurs, the only division I am using is dividing by the period, T, which is just pi.

I understand this is probably a bad way to do things! And I don't think I understand the symsum enough, but any help gratefully recieved.

clear all
close all hidden
clc
syms t % time
syms n % number of terms to calculate the sum to
syms T % period
syms a_n(t, n, T) % fourier co-effecients
syms b_n(t, n, T)
syms a_0(t, n, T)
syms f(t)

f = sin(t)^2

tvec = -3 : 0.01 : 3;
xvec = sin(tvec).^2;

figure
plot(tvec, xvec)
title('Plot of f(t)')

T = sym(pi)

a_0 = (2/T)*int(f, t, 0, T)/2

a_n = (2/T)*int(f*cos(((2*pi)/T)*n*t), t, 0, T)
b_n = (2/T)*int(f*sin(((2*pi)/T)*n*t), t, 0, T)

exp = a_0 + symsum(a_n*cos(((2*pi)/T)*n*t) + b_n*sin(((2*pi)/T)*n*t), n, 1, 5)
% division by zero error

[Torsten]

Look at a_n - you get a division by zero if n = 1 because you divide be (n^2-1).

It's a bug of the symbolic toolbox that the same case distinction that is made for b_n is not done for a_n.

[VBBV]

clear all
close all hidden
clc
syms t % time
syms n % number of terms to calculate the sum to
syms T % period
syms a_n(t, n) % fourier co-effecients
syms b_n(t, n)
syms a_0(t, n)
syms f(t)

f = sin(t)^2

f = 

tvec = -3 : 0.01 : 3;
xvec = sin(tvec).^2;

figure
plot(tvec, xvec)
title('Plot of f(t)')

T = sym(pi)

T = 

π

a_0 = (2/T)*int(f, t, 0, T)/2

a_0 = 

a_n = (2/T)*int(f*cos(((2*pi)/T)*n*t), t, 0, T)

a_n = 

b_n = (2/T)*int(f*sin(((2*pi)/T)*n*t), t, 0, T)

b_n = 

exp = a_0 + symsum(a_n*cos(((2*pi)/T)*n*t) + b_n*sin(((2*pi)/T)*n*t), n, 2, 5) % change n to 2

exp = 

% division by zero error

[VBBV]

It is indeed strange when n varied from 1 to 5 it produces division by zero error, but when you change n to 2, it doesnt.

[VBBV]

@Torsten, you are right, but interestingly, when sym(pi) is changed to pi , the summation behavior changes

[Torsten]

... but the division by zero remains.

pi is the numerical pi - thus you will get differences due to the precision of the representation of pi by a fracture.

From MATLAB documentation:

For example, the command a = sym('pi') creates a symbolic variable named pi and assigns it to the workspace variable a . To create a symbolic number representing the constant π, use a = sym(pi) instead.

[Owen]

Thanks, I see it, I'm not sure where I am going wrong as I have written code to do the same thing in python / sympy and it works. I know the series should expand to 1/2 - cos(2*t)/2. In my MATLAB code the average value is correct at 1/2 (a_0/2). I reset the limits to be about -T/2 and T/2 and now b_n evaluates to 0 as expected. When n is 1 a_n should evaluate to -1/2. I get this in my python code and of course 1/2 - cos(2*t)/2 is just an identity for sin(t)^2 so that all makes sense.

I have gone wrong somewhere in my MATLAB code..

[Owen]

Thanks, yes, I see why we try to divide by zero when subbing in n=1, subbing in n=2 avoids this but a_n evaluates to 0 for all other n, meaning we just end up with a straight line at the average value, 1/2! When n=1 a_n should evaluate to -1/2, which in the complete fourier series equation would give us 1/2 - cos(2t)/2. For some reason I am getting an unexpected result for my formula for a_n.

[Owen]

>> syms t

>> n = 1;

>> ans = 2/sym(pi) * int(sin(t)^2 * cos(2*pi/T * n * t), t, -T/2, T/2)

ans =

-1/2

Ok so if I set n to 1 BEFORE I integrate I get the expected answer. This is how I was doing it in the python code as well. However for simpler fourier series expansions I would calculate a_n as a function in terms of n and then plug in the n's to get the fourier co-effecients. I'm interested in why for this specific function we get an a_n that is invalid for n=1.

[Torsten]

I'm interested in why for this specific function we get an a_n that is invalid for n=1.

But I alrady answered it - it's a MATLAB bug.

[Owen]

Thank you, I had misread that it was specifically a bug that b_n is not also invalid, but I see you meant the other way round, thanks!

Answer 1

[Paul]

I don't think anything is wrong with the code in the Question, at least based on visual inspection (though it's always best to use sym(pi) rather than pi in symbolic expressions). Rather, int often doesn't recognize special cases when evaluating parameterized, Fourier integrals. The user needs to identify those cases and manually process them.

syms t % time
syms n % number of terms to calculate the sum to

f = sin(t)^2;
T = sym(pi);

a_0 = 2/T*int(f,t,0,T)/2

a_0 = 

a_n = 2/T*int(f*cos(2*sym(pi)/T*n*t),t,0,T)

a_n = 

b_n = 2/T*int(f*sin(2*sym(pi)/T*n*t),t,0,T)

b_n = 

I don't know why int() captures the special cases of n for b_n but does not for a_n. So we have to handle the special case for a_n manually (keeping in mind that we only care about n >= 1)

a_n = piecewise(n==1,limit(a_n,n,1),a_n)

a_n = 

%b_n = piecewise(n==1,limit(b_n,n,1),b_n)

We can simplify things by placing appropriate assumptions on n. This step isn't really necessary, but adds some clarity IMO.

assume(n,'integer')
assumeAlso(n,'positive'); % n >= 1 for the sin/cos Fourier series
a_n = simplify(a_n)

a_n = 

b_n = simplify(b_n)

b_n = 

0

The reconstruction yields the expected result

fr = a_0 + symsum(a_n*cos(2*sym(pi)/T*n*t),n,1,5)

fr = 

However, an incorrect result for a_n ensues if evaluating the integrals with the assumptions on n.

a_n = 2/T*int(f*cos(2*sym(pi)/T*n*t),t,0,T)

a_n = 

0

b_n = 2/T*int(f*sin(2*sym(pi)/T*n*t),t,0,T)

b_n = 

0

That looks like a bug.

Alternatively, we can Hold the integrals and then release() them after evaluating them at specified values of n (the @doc functionality isn't finding symbolic release() ?)

a_n(n) = 2/T*int(f*cos(2*sym(pi)/T*n*t),t,0,T,'Hold',true)

a_n(n) = 

b_n(n) = 2/T*int(f*sin(2*sym(pi)/T*n*t),t,0,T,'Hold',true)

b_n(n) = 

a_n = release(a_n(1:5))

a_n = 

b_n = release(b_n(1:5))

b_n = 

fr = a_0 + sum(a_n.*cos(2*sym(pi)/T*(1:5)*t))

fr = 

[Owen]

Thanks Paul this has cleared everything up and thank you for the work-arounds, great to have some suggestions of new functrionality for me to learn and glad to know that my code was not the problem!