Skip to content
Open in molab View on GitHub Download .py

Symbolic evaluation of Fourier series coefficients

from sympy import cos, init_printing, integrate, pi, plot, sin, symbols
init_printing(pretty_print=True,use_latex=True)
import matplotlib.pyplot as plt

f,t,T,G = symbols('f t T G')

f_1 = G * t
print('f = ', f_1)

f =  G*t

plot((f_1.subs(G, 25), (t, 0, 1)), ((f_1 - 25).subs(G, 25), (t, 1, 2)), ((f_1 - 50).subs(G, 25), (t, 2, 3)), ylabel='f = 25 t [V]', xlabel='t [sec]', xlim=(0, 3))

<sympy.plotting.backends.matplotlibbackend.matplotlib.MatplotlibBackend object at 0x7f19040fb380>
c0 = integrate(f_1, (t, 0, T)) * (1 / T)
print(c0)
print(c0.subs([(G, 25), (T, 1)]))

G*T/2
25/2

a1 = 2 / T * integrate(f_1 * sin(2 * pi * t), (t, 0, T))

print(a1)
print(a1.subs([(T,1),(pi, 3.14),(G,25)]))

2*G*(-T*cos(2*pi*T)/(2*pi) + sin(2*pi*T)/(4*pi**2))/T
-7.96178343949045

a, b = ([], [])
for _n in range(10):
    a.append(2 / T * integrate(f_1 * sin(_n * pi * t), (t, 0, T)))
    b.append(2 / T * integrate(f_1 * cos(_n * pi * t), (t, 0, T)))

a[9], b[9]
  • 0: ||[ \frac{2 G \left(- \frac{T \cos{\left(9 \pi T \right)}}{9 \pi} + \frac{\sin{\left(9 \pi T \right)}}{81 \pi^{2}}\right)}{T} ||]
  • 1: ||[ \frac{2 \left(G \left(\frac{T \sin{\left(9 \pi T \right)}}{9 \pi} + \frac{\cos{\left(9 \pi T \right)}}{81 \pi^{2}}\right) - \frac{G}{81 \pi^{2}}\right)}{T} ||]
c = []
for _n in range(1, 10):
    print(a[_n])
    expr = a[_n] ** 2 + b[_n] ** 2
    expr = expr.subs([(T, 1), (pi, 3.14), (G, 25)])
    c.append(expr)  # print expr  # val = lambdify(t,expr,'numpy')

2*G*(-T*cos(pi*T)/pi + sin(pi*T)/pi**2)/T
2*G*(-T*cos(2*pi*T)/(2*pi) + sin(2*pi*T)/(4*pi**2))/T
2*G*(-T*cos(3*pi*T)/(3*pi) + sin(3*pi*T)/(9*pi**2))/T
2*G*(-T*cos(4*pi*T)/(4*pi) + sin(4*pi*T)/(16*pi**2))/T
2*G*(-T*cos(5*pi*T)/(5*pi) + sin(5*pi*T)/(25*pi**2))/T
2*G*(-T*cos(6*pi*T)/(6*pi) + sin(6*pi*T)/(36*pi**2))/T
2*G*(-T*cos(7*pi*T)/(7*pi) + sin(7*pi*T)/(49*pi**2))/T
2*G*(-T*cos(8*pi*T)/(8*pi) + sin(8*pi*T)/(64*pi**2))/T
2*G*(-T*cos(9*pi*T)/(9*pi) + sin(9*pi*T)/(81*pi**2))/T

plt.figure()
plt.bar(range(1,10),c)
plt.xlabel('n')
plt.ylabel('Fourier coefficients $a_n^2 + b_n^2$')