f =

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.

41 views (last 30 days)

Show older comments

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
on 31 Mar 2024

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.

Paul
on 31 Mar 2024

Edited: Paul
on 1 Apr 2024

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_n = 2/T*int(f*cos(2*sym(pi)/T*n*t),t,0,T)

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

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)

%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)

b_n = simplify(b_n)

The reconstruction yields the expected result

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

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)

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

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)

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

a_n = release(a_n(1:5))

b_n = release(b_n(1:5))

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

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!