Ejemplo 8. Uso de matrices: Serie trigonométrica de Fourier.

A partir de este ejemplo empiezan otro tipo de cosas. Esta entrada va a ser un poco larga.

Haremos un programa que calcule una serie de Fourier para funciones definidas a trozos.

Empecemos por lo básico: matrices como estructura de datos.

Podemos ver a una matriz como un arreglo de vectores, en otras palabras, un vector de vectores 😫. ¿Cómo?

Imagina un vector simple de 3 elementos, representaré cada elemento así [ ]. Entonces tenemos lo siguiente:

[ ] -> Elemento 1

[ ] -> Elemento 2

[ ] -> Elemento 3

Vector negro

Ahora supongamos que creo otros tres vectores de 3 elementos:

[ ][ ][ ] -> Vector verde

[ ][ ][ ] -> Vector azul

[ ][ ][ ] -> Vector rojo

Puedo colocar cada uno de los tres nuevos vectores en cada uno de los elementos del vector negro (el inicial), es decir, el vector verde en el elemento 1 del vector negro, el vector azul en el elemento 2 del vector negro y el vector rojo en el elemento 3 del vector negro:

[ ][ ][ ]

Como puedes observar el elemento 1 del vector negro contendrá el mismo valor que el elemento 1 del vector verde, el elemento 2 del vector negro el mismo valor que el elemento 1 del vector azul y el elemento 3 del vector negro el valor del elemento 1 del vector rojo.

De esta manera tenemos una matriz, la cual consta de filas (renglones, horizontal) y columnas (vertical).

Por ejemplo, la siguiente es una matriz de 4 filas y 5 columnas a la que normalmente llamamos matriz de 4x5:

[ ][ ][ ][ ][ ] -> Fila 1

[ ][ ][ ][ ][ ] -> Fila 2

[ ][ ][ ][ ][ ] -> Fila 3

[ ][ ][ ][ ][ ] -> Fila 4

Podemos identificar cada elemento con un par ordenado, en este el primer numero indica la fila y el segundo indica la columna. Para la matriz anterior tenemos lo siguiente:

[ 1,1 ][ 1,2 ][ 1,3 ][ 1,4 ][ 1,5 ]

[ 2,1 ][ 2,2 ][ 2,3 ][ 2,4 ][ 2,5 ]

[ 3,1 ][ 3,2 ][ 3,3 ][ 3,4 ][ 3,5 ]

[ 4,1 ][ 4,2 ][ 4,3 ][ 4,4 ][ 4,5 ]

Por ejemplo, para acceder al elemento de la esquina superior derecha usaremos el índice 1,5, para el elemento de la esquina inferior izquierda el índice 4,1.

Ahora veamos unos comandos que serán útiles para realizar el ejemplo.

La dimensión de una matriz será un par ordenado que indica el número total de filas y columnas de la matriz. El comando:

dim(matriz)

regresa una lista con ese par ordenado:

Si queremos conocer únicamente la dimensión en columnas existe el comando:

colDim(matriz)

que regresa el número total de columnas.

Para la dimensión en Filas se usa el comando:

rowDim(matriz)

que de igual manera, regresa el número total de filas.

Para acceder a una fila, se hace lo que muestra la siguiente imagen:

Se escribe un vector con un único elemento justo delante de la matriz. El elemento contiene el número de la fila.

Para acceder a un elemento, se hace algo bastante parecido:

Se escribe un vector de dos elementos justo delante de la matriz. El primer elemento es el número de la fila y el segundo el número de la columna.

Algo que también vamos a usar es la siguiente forma de operar:
Tenemos una variable a con el valor de 3.

Puedo hacer lo siguiente:
A la variable a, le asignaré la suma de el valor de la variable a y el valor de la variable a.
La suma nos da un resultado de 6, y ese valor se asigna a a.

Si vuelvo a hacerlo, se hará la suma del nuevo valor de a, que es 6, y deberá dar 12.

Ahora a=12. Hago de nuevo la suma, debe resultar 24.

De esta manera "se actualiza" el valor de la variable, relacionando el valor que actualmente tiene. Mira esto:

a comienza con valor de 3, le sumo 1 y asigno a la variable a. Ahora a=4.
Sumo 1 y lo asigno a a. Como a=4, entonces 4+1=5 y lo asigno a a. a=5.
Sumo 1 y de nuevo lo asigno a a. a=6.

Veamos un poco de teoría de las series de Fourier. Una serie de Fourier es una serie infinita que sirve para aproximar una función periódica y continua a una función periódica y definida por partes.

En nuestro caso usaremos la Serie Trigonométrica de Fourier. Su característica principal se encuentra en su nombre, pues la serie se compone principalmente de las funciones Seno y Coseno.

Algunos de los mejores ejemplos se encuentran en análisis de circuitos multifrecuencia, aunque claro, en toda la física puede utilizarse, por ejemplo, en el análisis de ondas en una cuerda y todo lo que esté relacionado con movimiento armónico.

La Serie Trigonométrica de Fourier (desde ahora STF) tiene la siguiente forma:

$$f(x)=\frac{a_0}{2}+\sum_{n=1}^{\infty}a_n\cdot cos(n\cdot \omega \cdot x)+\sum_{n=1}^{\infty}b_n\cdot sin(n\cdot \omega \cdot x)$$

Donde:

$\omega=\frac{2\pi}{T}=2\pi \cdot f$

$a_0$, $a_n$ y $b_n$ son los coeficientes de Fourier:

$a_0=\frac{2}{T}\int_{0}^{T}f(x) dx$

$a_n=\frac{2}{T}\int_{0}^{T}f(x)\cdot cos(\omega\cdot n \cdot x) dx$

$b_n=\frac{2}{T}\int_{0}^{T}f(x)\cdot sin(\omega\cdot n \cdot x) dx$

Para $n=1,2,3,...$

Como sabes, debe trabajarse con una función definida a trozos y periódica. Aquí hay unos ejemplos:

Obviamente la onda senoidal no aplica, pues es continua. Pero la onda cuadrada, la triangular y diente de sierra se trata de funciones definidas por partes.

Aunque ya tenemos las fórmulas que necesitamos, resolvamos un ejemplo sencillo para obtener el algoritmo. Vamos a aproximar la siguiente onda cuadrada con una STF.

Primero, el periodo T=4, ya que para completar ese ciclo comienza en x=0 y termina en x=4, entonces 4 - 0 = 4, es decir $T=x_{final}-x_{inicial}$

Calculemos los coeficientes:

$a_0=\frac{2}{T}\int_{0}^{T}f(x) dx = \frac{1}{4}\int_{0}^{T}f(x) dx $

Como $f(x)$ está definida a trozos, la integral se "separa" en todas esas partes a manera de cubrir el periodo de la siguiente forma:

$a_0=\frac{2}{4}\cdot (\int_{0}^{2}0\cdot dx +\int_{2}^{4} 1\cdot dx )$

La primera integral corresponde a la parte de la función definida en el intervalo $x\in (0,2]$. La segunda corresponde a la parte de la función definida en $x\in (2,4]$, de ahí "salen" los límites de integración y ese cero y uno respectivamente.

La primer integral resulta cero. Resolvemos la segunda:

$a_0=\frac{2}{4}\cdot (\int_{2}^{4} 1\cdot dx =1)$

Ahora el coeficiente $a_n$:

$a_n=\frac{2}{4}\cdot(\int_{0}^{2}0\cdot cos(\omega\cdot n \cdot x) dx +\int_{2}^{4}1\cdot cos(\omega\cdot n \cdot x) dx)$

$\omega=\frac{2\pi}{4}=\frac{\pi}{2}$

$a_n=\frac{2}{4}\cdot(\int_{2}^{4}1\cdot cos(\frac{\pi}{2}\cdot n \cdot x) dx)=\frac{2sin(2n\pi)}{n\pi}-\frac{2sin(n\pi)}{n\pi}$

*$sin(2n\pi)|_{n=1,2,3,...}=0$

*$sin(n\pi)|_{n=1,2,3,...}=0$

Por lo tanto $a_n=0$

Por último el coeficiente $b_n$:

$b_n=\frac{2}{4}\cdot(\int_{0}^{2}0\cdot sin(\omega\cdot n \cdot x) dx+\int_{2}^{4}1\cdot sin(\omega\cdot n \cdot x) dx)$

$\omega=\frac{\pi}{2}$

$b_n=\frac{2}{4}\cdot(\int_{2}^{4}1\cdot sin(\omega\cdot n \cdot x) dx)=\frac{-cos(2n\pi)-cos(n\pi)}{n\pi}$

Ahora sólo falta "escribir" la serie de Fourier:

$$f(x)=\frac{1}{2}+\sum_{n=1}^{\infty}0\cdot cos(n\cdot \omega \cdot x)+\sum_{n=1}^{\infty}(\frac{-cos(2n\pi)-cos(n\pi)}{n\pi}\cdot sin(n\omega \cdot x))$$

$$f(x)=\frac{1}{2}+\sum_{n=1}^{\infty}(\frac{-cos(2n\pi)-cos(n\pi)}{n\pi}\cdot sin(n \frac{\pi}{2} \cdot x))$$

Y listo, tal vez es un proceso largo, pero aprovecharemos las funciones de la calculadora para hacer integrales definidas.

El algoritmo también se complica, pero al final obtener una serie de Fourier será tan rápido, incluso para funciones por partes muy grandes.

Tómate un respiro, porque vamos por ese algoritmo.

Veamos como vamos a introducir nuestra función a trozos. Yo propongo lo siguiente:

Se crea una matriz de n filas y 3 columnas:

[ ][ ][ ] -> f1

[ ][ ][ ] -> f2

[ ][ ][ ] -> f3

...

[ ][ ][ ] -> fn

La primer columna tendrá la función definida en un intervalo.

La segunda el límite inferior del intervalo y la tercera el límite superior.

Para la función anterior:

La matriz queda de la siguiente forma:

[ 0 ][ 0 ][ 2 ] -> 0 para $0<x\leq2$

[ 1 ][ 2 ][ 4 ] -> 1 para $2<x\leq4$

Otro ejemplo:

La matriz queda:

[ 0 ][ 0 ][ 1 ] -> 0 para $0<x\leq1$

[ 1 ][ 1 ][ 2 ] -> 1 para $1<x\leq2$

[ 2 ][ 2 ][ 3 ] -> 2 para $2<x\leq3$

[ 3 ][ 3 ][ 4 ] -> 3 para $3<x\leq4$

El número de filas corresponde con el número de funciones que conforman f(x), en este caso 4 funciones -> 4 filas. SIEMPRE TENDRÁ 3 COLUMNAS (límite inferior y límite superior).

Comenzamos agregando un editor de programas:

Vamos a declarar el argumento del programa:

Recibe una matriz.
Recibe un número j de "aproximaciones".

Las variables las iremos declarando como se necesiten, pero por ahora podemos decir que tenemos las siguientes:

t -> Almacena el periodo.
d -> Almacena la dimensión en Filas de la matriz.
w -> Almacena la frecuencia angular (no escribo $\omega$ para ahorrar tiempo).
a0 -> Coeficiente de Fourier independiente.
an -> Coeficiente de Fourier.
bn -> Coeficiente de Fourier.

Comencemos obteniendo la dimensión en Filas de matriz, la cual se almacena en la variable d. Usamos el comando rowDim().

Ahora hay que calcular el periodo $T$. Como vimos arriba, está dado por:

$$T=x_{final}-x_{inicial}$$

Imaginemos la siguiente matriz 4x3:

[ f1 ][ x0 ][ x1 ]

[ f2 ][ x1 ][ x2 ]

[ f3 ][ x2 ][ x3 ]

[ f4 ][ x3 ][ x4 ]

$x_{inicial}$ está en el elemento 1,2. $x_{final}$ está en el elemento 4,3.

Ahora otra matriz 6x3:

[ f1 ][ x0 ][ x1 ]

[ f2 ][ x1 ][ x2 ]

[ f3 ][ x2 ][ x3 ]

[ f4 ][ x3 ][ x4 ]

[ f5 ][ x4 ][ x5 ]

[ f6 ][ x5 ][ x6 ]

$x_{inicial}$ está en 1,2. $x_{final}$ en 6,3.

Como se observa, $x_{inicial}$ siempre estará en el elemento 1,2, mientras que $x_{final}$ en la columna 3 de la última fila, es decir, si guardamos la dimensión en filas en la variable d, el último valor de x estará en el elemento d,3. Calculamos $T$:

Mostremos ese resultado en pantalla para comprobar:

Para esa función el periodo es de 2, como muestra el resultado. Continuemos.

Calculamos la frecuencia angular $\omega=\frac{2\pi}{T}$, y la mostramos:

Ahora a calcular el primer coeficiente $a_0=\frac{2}{T}\int_{0}^{T}f(x) dx$. Como vimos la integral se va a separar en más integrales, dependiendo del número de funciones definidas en los intervalos. Imagina la siguiente matriz 2x3:

[ f1 ][ x0 ][ x1 ]

[ f2 ][ x1 ][ x2 ]

Son dos funciones, por lo tanto la integral se descompone en dos integrales:

$$a_0=\frac{2}{T}\cdot(\int_{x0}^{x1}f1 dx+\int_{x1}^{x2}f2 dx)$$

Expresemos f1, f2, x0, x1 y x2 como elementos de la matriz:

$$a_0=\frac{2}{T}\cdot(\int_{[1,2]}^{[1,3]}[1,1] dx+\int_{[2,1]}^{[2,3]}[2,1] dx)$$

Para la matriz 3x3:

[ f1 ][ x0 ][ x1 ]

[ f2 ][ x1 ][ x2 ]

[ f3 ][ x2 ][ x3 ]

Nuestra integral queda:

$$a_0=\frac{2}{T}\cdot(\int_{[1,2]}^{[1,3]}[1,1] dx+\int_{[2,1]}^{[2,3]}[2,1] dx+\int_{[3,2]}^{[3,3]}[3,1] dx)$$

En una matriz de d funciones:

[ f1 ][ x0 ][ x1 ]

[ f2 ][ x1 ][ x2 ]

[ f3 ][ x2 ][ x3 ]

...

[ fd ][ x(d-1) ][ xd ]

$$a_0=\frac{2}{T}\cdot(\int_{[1,2]}^{[1,3]}[1,1] dx+\int_{[2,1]}^{[2,3]}[2,1] dx+...+\int_{[d,2]}^{[d,3]}[d,1] dx)$$

Esto nos lleva a tener que hacer iteraciones igual al número de funciones, el cual está almacenado en la variable d.

Haremos uso de la forma de operar:

$$variable = variable + valor$$$

Al final, el coeficiente $a_0$ será la suma de cada una de las integrales en las que se descompuso la integral principal. Por lo tanto haremos un ciclo que haga esa suma. Necesitamos una variable iteradora, la declaramos en la parte del comando Local.

Inicializamos la variable a0, asignando el valor de cero.

NOTA: Esto se hace porque el cero es el "neutro" de la suma, es decir, si yo hago a+0=a. No afecta nuestra suma y es lo que necesitamos para nuestra suma iterada.

Codificaré el ciclo, que inicia en 1 y termina con el número de ecuaciones o filas d:

Finalmente, la suma de cada integral de la función definida a trozos:

Como se ve en la fórmula, el resultado de la parte de las integrales se multiplica por $\frac{2}{T}$, entonces:

Hecho esto ya tenemos el coeficiente $a_0$. Mostremoslo para comprobar.

Vamos por el coeficiente $a_n=\frac{2}{T}\int_{0}^{T}f(x)\cdot cos(\omega\cdot n \cdot x) dx$

De nuevo, hay que ponerlo en términos de elementos de una matriz:

[ f1 ][ x0 ][ x1 ]

[ f2 ][ x1 ][ x2 ]

[ f3 ][ x2 ][ x3 ]

...

[ fd ][ x(d-1) ][ xd ]

$$a_n=\frac{2}{T}\cdot(\int_{[1,2]}^{[1,3]}[1,1]\cdot cos(\omega\cdot n \cdot x) dx+\int_{[2,2]}^{[2,3]}[2,1]\cdot cos(\omega\cdot n \cdot x) dx)+...+\int_{[d,2]}^{[d,3]}[d,1]\cdot cos(\omega\cdot n \cdot x) dx$$

Es prácticamente lo mismo que el coeficiente $a_0$, sólo se agrega $cos(\omega \cdot n\cdot n)$:

De nuevo, se inicializa an con el valor de cero y usamos el mismo ciclo, agregando la parte indiciada a la integral:

También se multiplica el resultado de la suma por $\frac{2}{T}$:

Listo, mostremos el resultado para probar lo que llevamos.

Hagamos lo mismo para el coeficiente $b_n$, sólo copiemos y peguemos el código de $a_n$ y cambiemos $cos(\omega n x)$ por $sin(\omega n x)$, y obviamente reemplazar ab por bn:

Probemos:

Parece que todo va bien.

Una vez que se tienen los coeficientes se procede a escribir la serie.

$$f(x)=\frac{a_0}{2}+\sum_{n=1}^{\infty}a_n\cdot cos(n\cdot \omega \cdot x)+\sum_{n=1}^{\infty}b_n\cdot sin(n\cdot \omega \cdot x)$$

Entonces declaramos una variable f que almacenará dicha serie:

Ahora a escribir la serie y mostrarla:

Se ve diferente la imagen porque la fórmula no cabía en la pantalla normal. Observa que las sumas (notación sigma) comienzan en 1 y terminan en j. Por definición la serie de Fourier es infinita, pero computacionalmente no podemos hacer eso. Se opta por elegir un número muy grande, que en nuestro caso será j.

Por cuestiones de potencia de la calculadora j no puede ser muy grande, el proceso se ralentiza e incluso nos veremos en la necesidad de reiniciar la calculadora.

A probar con nuestro ejemplo: