


# Series de Fourier en Forma Compleja

[](https://colab.research.google.com/github/ivan-jgr/computacion-cientifica/blob/main/Laboratorios/Laboratorio-5.ipynb)

Definimos la serie de Fourier de una función periodica $f(t)$ de la siguiente forma:

$$
\hat{f}(t) = \sum_{n = -\infty}^\infty c_n e^{i \omega_n t},
$$

donde $c_n = \displaystyle \dfrac{1}{2L} \int_{-L}^L f(t) e^{-i \omega_n t} dt$.

En este laboratorio utilizaremos la librería `cmath` para manipular numeros complejos, calcular la serie de Fourier en forma compleja y utilizarla para resolver ecuaciones diferenciales ordinarias.


```python
import cmath # Libreria para manipular complejos
```

Para crear un número complejo $z = a + jb$ (note que se usa $j$ para representar la unidad imaginaria) hacemos `z = complex(a, b)`

### Aritmética de Complejos


```python
z = complex(3, 2) # creamos el complejo 3 + 2j
w = complex(1, -2) # creamos el complejo 1 - 2j

# puede obtener a la parte real de un complejo z haciendo z.real
print(f"La parte real de z es: {z.real}") 
# puede obtener a la parte imaginaria de un complejo z haciendo z.imag
print(f"La parte imaginaria de z es: {z.imag}")
# Puede calcular el módulo de un complejo w haciendo abs(w)
print(f"El módulo de w es: {abs(w)}")
# Puede calcular el argumento (en radianes) de un complejo z haciendo phase(z)
print(f"El argumento de z es: {cmath.phase(z)}")
# puede sumar, restar, multiplicar  dividir complejo
print(f"El resultado de z + w es: {z + w}")
print(f"El resultado de zw es: {z*w}")
```

    La parte real de z es: 3.0
    La parte imaginaria de z es: 2.0
    El módulo de w es: 2.23606797749979
    El argumento de z es: 0.5880026035475675
    El resultado de z + w es: (4+0j)
    El resultado de zw es: (7-4j)



```python
cn = complex(3, 2) # Creamos el complejo 3 + 2j
cn_conjugado = complex(cn.real, -cn.imag) # su conjugado
print('El conjugado de  {} es {}'.format(cn, cn_conjugado))

# También podemos crear complejos de la siguiente forma
j = complex(0, 1)
dn = 3 + 2*j # cuidado con redefinir j

print('Complejo dn = {}'.format(dn))
```

    El conjugado de  (3+2j) es (3-2j)
    Complejo dn = (3+2j)


### Serie de Fourier en Forma Compleja

#### <ins> Problema 1 </ins>

1. Implemente la función `cn(f, L, n)` que calcula en coeficiente de la serie Compleja de Fourier:
\begin{equation}c_n = \dfrac{1}{2L} \int_{-L}^L f(t) e^{-i \omega_n t} dt. \label{cn} \tag{1} \end{equation}

 Recuerde que si $f(t)$ es compleja, i.e. es de la forma $f(t) = u(t) + i v(t)$, entonces: $$\int f(t) dt = \int \operatorname{Re}(f(t)) dt + i \int \operatorname{Im}(f(t)) dt$$

 Es decir que puede integrar la parte real del integrado, luego la parte imaginaria y luego formar el complejo resulte. El método `quad` que hemos usado para integrar **no** puede manejar terminos complejos, por lo que deberá integrar por separado las partes real e imaginaria de \eqref{cn}.


2. Implemente la función `Sn(f, L, N)` que calcula la $N$-ésima suma parcial de Fourier:
$$S_N (t) = \sum_{n = -N}^N c_n e^{i \omega_n t}.$$


```python
# la función exponencial la puede encontrar como np.exp y sí acepta valores complejos
import numpy as np
from scipy.integrate import quad

def cn(f, L, n):
    """
    Calcula el coeficiente cn para la serie de Fourier compleja.
    
    :param f: es la una función
    :param L: es el semiperíodo de f
    :param n: es un entero (positivo o negativo)
    
    returns: el coeficiente cn de Fourier descrito por la formula anterior.
    """
    pass
```


```python
def Sn(f, L, N):
    """
    Devuelve la suma parcial de Fourier.
    """
    pass
```


```python
# Verifiquemos que Sn nos devuelva complejos con parte imaginaria igual a cero
f_test = lambda x: abs(x)

f_test_serie = Sn(f_test, 1, 10)

t = np.linspace(-1, 1)

f_test_serie(t).imag
```




    array([-8.20588570e-18, -1.65326325e-03, -5.27721283e-04,  1.48481473e-03,
            1.00167383e-03, -1.16508065e-03, -1.37356733e-03,  7.26638248e-04,
            1.60551011e-03, -2.14159758e-04, -1.67386986e-03, -3.20139149e-04,
            1.57168151e-03,  8.21819557e-04, -1.30935688e-03, -1.23976600e-03,
            9.13623833e-04,  1.53139453e-03, -4.24803019e-04, -1.66699157e-03,
           -1.07300339e-04,  1.63274132e-03,  6.28471027e-04, -1.43213350e-03,
           -1.08560774e-03,  1.08560774e-03,  1.43213350e-03, -6.28471027e-04,
           -1.63274132e-03,  1.07300339e-04,  1.66699157e-03,  4.24803019e-04,
           -1.53139453e-03, -9.13623833e-04,  1.23976600e-03,  1.30935688e-03,
           -8.21819557e-04, -1.57168151e-03,  3.20139149e-04,  1.67386986e-03,
            2.14159758e-04, -1.60551011e-03, -7.26638248e-04,  1.37356733e-03,
            1.16508065e-03, -1.00167383e-03, -1.48481473e-03,  5.27721283e-04,
            1.65326325e-03,  8.20588570e-18])



Note que la parte imaginaria de todos los valores valuados en la serie es cero (o muy cercana a cero por la precisión del metodo `quad`) ¿Por qué?


```python
import matplotlib.pyplot as plt
from matplotlib import rcParams

rcParams['figure.figsize'] = 10, 6
plt.style.use('seaborn-whitegrid')

t = np.linspace(-2, 2, 500)
plt.plot(t, f_test_serie(t).real) # Note que graficamos valores reales
plt.show()
```




    [<matplotlib.lines.Line2D at 0x7fc67ececbd0>]



## Serie Compleja de Fourier para resolver EDOs

Nuestro objetivo es resolver la ecuación diferencial (encontrar la solución al estado estable para la carga $q(t)$) para un circuito RLC con una fuente de voltaje descrito por una función periódica $v(t)$:



La carga en el capacitor está modelada por la ecuación:
\begin{equation} L \frac{d^2q}{dt^2} + R\frac{dq}{dt} + \frac{1}{C} q = v(t). \label{eq:edo} \tag{1}\end{equation}

Podemos expandir $v(t)$ en su serie compleja de Fourier: $v(t) = \displaystyle\sum_{n = -\infty}^\infty c_n e^{i \omega_n t}$ y al usar el principio de superposición resolver
$$ L \frac{d^2q}{dt^2} + R\frac{dq}{dt} + \frac{1}{C} q = e^{i \omega_n t},$$

que podemos resolver usando coeficientes indeterminados:

- Proponemos $q_n(t) = A_n e^{i \omega_n t}$ y obtendremos que $A_n = \dfrac{C}{(1-LC\omega_n^2) + RC\omega_n i}$

Por simplicidad asumamos que $f(t)$ tiene período $T = 2 \pi$ y entonces:

$$A_n = \dfrac{C}{(1-LCn^2) + RCn i}.$$

La solución a la ecuación $\eqref{eq:edo}$ es entonces:

$$q(t) = \sum_{n = -\infty}^\infty c_n A_n e^{i \omega_n t} = \sum_{n = - \infty}^\infty c_n \dfrac{C}{(1-LCn^2) + RCn i} e^{i n t}.$$


```python
# valores de L, R, C
L, R, C = 1, 1, 1

def An(n):
    # denominador de An
    den = complex(1-L*C*n**2, R*C*n)
    return C / den

print(f"El coeficiente A3 = {An(3)}")
```

    El coeficiente A3 = (-0.1095890410958904-0.0410958904109589j)



```python
def f(t):
    if t > 0:
        return 1
    else: 
        return -1

# Terminos de la serie
N = 30
    
q_sol = lambda x: sum([cn(f, np.pi, k)*An(k)*np.exp(-j*k*x) for k in range(-N-1, N+1)])
```


```python
rcParams['figure.figsize'] = 10, 6

v = Sn(f, np.pi, 30)

t = np.linspace(0, 6*np.pi, 200)
plt.plot(t, v(t).real, label='$v(t)$', color='green')

plt.plot(t, q_sol(t).real, label='$q(t)$', color='red')

plt.xlabel('Tiempo en segundos')
plt.legend()
plt.show()
```


```python
# Veamos como evoluciona la solución al variar el tiempo
from ipywidgets import interact, FloatSlider
rcParams['figure.figsize'] = 10, 6

def plot_sol(tmax):
    t = np.linspace(0, tmax, 200)
    plt.plot(t, q_sol(t).real, label='$q(t)$', color='red')
    plt.plot(t, v(t).real, label='$v(t)$', color='green')
    plt.scatter(tmax, q_sol(tmax).real, color='red')

    plt.ylim(-2, 2)
    plt.xlabel('Tiempo en segundos')
    plt.legend(loc='upper left')
    plt.show()
    
interact(plot_sol, tmax= FloatSlider(min = 0.01, max = 50, step=0.5));
```


    interactive(children=(FloatSlider(value=0.01, description='tmax', max=50.0, min=0.01, step=0.5), Output()), _d…



```python
from matplotlib import rcParams
rcParams['figure.figsize'] = 12, 8

def plot_sol(tmax):
    t = np.linspace(0, tmax, 200)
    
    ###########################################
    # fila 1 columna 1: Gráfica de la entrada
    plt.subplot(4, 2, 1)
    plt.plot(t, v(t).real, label='$v(t)$', color='green')
    
    plt.ylim(-2, 2)
    plt.legend(loc='upper left')
    
    # fila 1 columna 2: Gráfica de la salida
    plt.subplot(4, 2, 2)
    plt.plot(t, q_sol(t).real, label='$q(t)$', color='red')
    
    plt.ylim(-2, 2)
    plt.legend(loc='upper left')
    
    ###########################################
    # fila 2 columna 1: Primer armónico
    plt.subplot(4, 2, 3)
    
    v1 = cn(f, np.pi, 1)*np.exp(-j*1*t) + cn(f, np.pi, -1)*np.exp(-j*(-1)*t) 
    plt.plot(t, v1.real, label='primer armónico', color='blue')
    
    plt.ylim(-2, 2)
    plt.legend(loc='upper left')
    
     # fila 2 columna 2: Solucion al primer armónico
    plt.subplot(4, 2, 4)
    
    q1 = cn(f, np.pi, 1)*An(1)*np.exp(-j*1*t) + cn(f, np.pi, -1)*An(-1)*np.exp(-j*(-1)*t) 
    plt.plot(t, q1.real, label='$q_1(t)$', color='orange')
    
    plt.ylim(-2, 2)
    plt.legend(loc='upper left')
    
    
    ###########################################
     # fila 3 columna 2: Tercer armónico
    plt.subplot(4, 2, 5)
    
    v3 = cn(f, np.pi, 3)*np.exp(-j*3*t) + cn(f, np.pi, -3)*np.exp(-j*(-3)*t) 
    plt.plot(t, v3.real, label='tercer armónico', color='black')
    
    plt.ylim(-2, 2)
    plt.legend(loc='upper left')
    
    # fila 3 columna 2: Solucion al tercer armónico
    plt.subplot(4, 2, 6)
    
    q3 = cn(f, np.pi, 3)*An(3)*np.exp(-j*3*t) + cn(f, np.pi, -3)*An(-3)*np.exp(-j*(-3)*t) 
    plt.plot(t, q3.real, label='$q_5(t)$', color='gray')
    
    plt.ylim(-2, 2)
    plt.legend(loc='upper left')
    
    ###########################################
    # fila 4 columna 1: Quinto armonico
    plt.subplot(4, 2, 7)
    
    v5 = cn(f, np.pi, 5)*np.exp(-j*5*t) + cn(f, np.pi, -5)*np.exp(-j*(-5)*t) 
    plt.plot(t, v5.real, label='quinto armónico', color='y')
    
    plt.ylim(-2, 2)
    plt.legend(loc='upper left')
    
    # fila 4 columna 2: Solucion al quinto armonico
    plt.subplot(4, 2, 8)
    
    q5 = cn(f, np.pi, 5)*An(5)*np.exp(-j*5*t) + cn(f, np.pi, -5)*An(-5)*np.exp(-j*(-5)*t)
    plt.plot(t, q5.real, label='$q_5(t)$', color='c')
    
    plt.ylim(-2, 2)
    plt.legend(loc='upper left')
   
   
    
interact(plot_sol, tmax= FloatSlider(min = 0.01, max = 50, step=0.5, value=22));
```


    interactive(children=(FloatSlider(value=22.0, description='tmax', max=50.0, min=0.01, step=0.5), Output()), _d…


Las gráficas anteriores muestran la respuesta a $v(t)$ en la primera fila y luego las respuestas al primer, tercer y quinto armónico. La serie de Fourier nos permite ver como responde el circuito a los componentes de frecuencia específicos en la función de entrada como si fuera un prima que nos permite descomponer la solución $q(t)$ en las respuestas a cada componente de frecuencia presente en $v(t)$.

Además, note la diferencia en amplitud para $q_1(t)$ con respecto a los otros armónicos. ¿Puede explicar por qué la amplitud de este armónico es mayor que el resto?

#### <ins> Problema 2 </ins>

Considere el siguiente sistema masa-resorte:



se aplica una fuerza externa periódica dada por $f(t) = \left| \sin t \right|$ con período $T = 2 \pi$.
1. Encuentre la respuesta (solución al estado estable) para la posición de la masa $x(t)$ usando superposición con $\operatorname{N} = 30$ terminos para la serie.
 - $m$ es la mada del bloque
 - $b$ es la constante de amortigüamiento viscoso
 - $k$ es la constante del resorte

2. Encuentre una expresión en forma de serie para $v(t) = \dfrac{dx}{dt} = \dfrac{d}{dt} \displaystyle \left( \sum_{n = -\infty}^\infty c_n A_n e^{i \omega_n t} \right)$, i.e. derívela término por término.

3. Genere gráficas para $x(t), v(t)$ y $f(t)$.


4. Grafique la solución para el segundo y cuarto armónico.


```python
# Defina la función den entrada f(t)
def f(t):
    pass

# Encuentre la expansion en serie de f(t) y guardela en la variable f_hat
f_hat = None

# Utilíce coeficientes indeterminados y encuentre la solucion al estado estable para la posición de la masa
# Encuentre el coeficiente An como en el ejemplo mostrado
def An(n):
    pass


N = 30
# Escriba la solución como usa serie compleja (desde -N hasta N)
x_sol = None

# Escriba v(t) en forma de serie
v_sol = None
```


```python
# Grafique f(t) junto con x_sol(t) y v_sol(t) en el intervalo [0, 50]

```


```python
# Grafique la respuesta para el segundo y cuarto armónico
```

<hr>

# **Instrucciones Generales**
1. Este laboratorio será realizado de manera *individual*.
2. **Fecha de entrega**: Viernes 17 de Septiembre de 2021. Su solución debe subirla en un archivo ZIP enviado por GES y debe contener el archivo .ipynb con sus respuesta a cada inciso solicitado en cada uno de los *Problemas* indicados en la parte superior.

11f9c6d291304f765234d94f392893b33c8d5ea1
ipynb
Jupyter Notebook
Laboratorios/Laboratorio-5.ipynb
ivan-jgr/computacion-cientifica
a235b537cded383e6962720ed4b55d1d75c6c472
MIT
Laboratorios/Laboratorio-5.ipynb
ivan-jgr/computacion-cientifica
a235b537cded383e6962720ed4b55d1d75c6c472
MIT
Laboratorios/Laboratorio-5.ipynb
ivan-jgr/computacion-cientifica
a235b537cded383e6962720ed4b55d1d75c6c472
MIT
Qwen/Qwen-72B
1. YES
2. YES