Periodización de funciones. Coeficientes de Fourier. Aproximación por polinomios de Fourier.

Francisco Javier Pérez González
Departamento de Análisis Matemático
Universidad de Granada

Inicialización

Unprotect[Plot];
Plot[f_, {x_, lims__ /; Length[{lims}] > 2}, 
        opts___] := 
Block[{intervalos, discont, graf, 
          lineas, us}, 
intervalos = 
({#[[1]] +  (10. ^(-15)), #[[2]] -  
                    (10. ^(-15))} &) /@ 
            Partition[{lims}, 2, 
              1]; 
graf = (Plot[f, 
                  Evaluate[
                    Prepend[#, 
                      x]], 
DisplayFunction -> 
                    Identity, opts] &) /@ 
            intervalos; 
discont = (Drop[
                  Drop[{lims}, 1], -1] /. 
                us_?NumberQ -> {us - (10. ^(-15)), 
                    us + (10. ^(-15))}) // 
            Flatten; 
discont = 
          discont /. 
            us_?NumberQ ->  {us, 
                f /. x -> us}; 
lineas = 
          Line /@ Partition[discont, 2]; 
Show[{graf, 
            Graphics[{Dashing[{ .01,  .01}], 
                lineas}]}, 
DisplayFunction -> $DisplayFunction]]
Protect[Plot];

Introducción

En esta práctica vamos a usar Mathematica para realizar cálculos y representaciones gráficas que nos ayudarán a comprender mejor los conceptos básicos de la teoría de series de Fourier. La práctica consta de dos partes. En la parte primera estudiamos las distintas formas usuales de extender periódicamente una función cualquiera. En la segunda parte, más extensa que la primera, te voy a proponer que calcules unas veces de forma simbólica y otras de forma numérica  los coeficientes de Fourier de unas cuantas funciones y que hagas algunas animaciones gráficas para comprobar cómo los polinomios de Fourier se van aproximando a la función. Además deberás representar en cada caso los espectros de amplitudes y de fases. 
La razón por la que te propongo hacer todo esto en un solo cuaderno de Mathematica es la gran unidad e interdependencia que hay en todo lo que vamos a hacer. Por ello me ha parecido que sería más cómodo tenerlo todo junto en lugar de repartir el contenido en varios cuadernos, lo que haría más incómodo el trabajo pues tendríamos que tenerlos todos abiertos para hacer frecuentes consultas y referencias.  
Esta práctica es en realidad dos prácticas (por lo menos...). Cada parte de la práctica se entregará en las fechas que se indiquen en clase y en la página Web de la asignatura.

Parte primera. Periodización de funciones

Recuerda que una función \(TraditionalForm\`f : ℝ -> ℝ\) se llama periódica si hay un número T ≠  0tal que f(x + T) = f(x) para todo x ∈ ℝ en cuyo caso se dice que T es un periodo de f. Si f es periódica con periodo T se verifica que f(x + n T) = f(x) para todo entero n por lo que n T también es un periodo de f. Las funciones constantes son funciones periódicas que admiten como periodo cualquier número T ≠ 0. Salvo este caso, si f es una función periódica no constante y que al menos es continua en un punto (condición que cumple cualquier función "interesante") se verifica que f tiene un periodo positivo mínimo. Es a este periodo al que solemos referirnos cuando usamos la expresión "f es una función periódica con periodo T". Así, por ejemplo, decimos que la función seno es periódica con periodo 2π y la función tangente es periódica con periodo π. Precisamente, el ejemplo de la función tangente nos indica que no es preciso exigir que una función periódica esté definida en todo ℝ, es suficiente que esté definida en un conjunto A ⊂ ℝ, por lo general una unión de intervalos, tal que si x ∈ A también x + T ∈ A.

Para representar una función f por una serie de Fourier debemos elegir un intervalo acotado I en el que queremos que la serie represente a la función y debemos considerar que fuera de dicho intervalo la función f repite periódicamente los valores que toma en el mismo; esto es, consideramos a f como una función periódica con periodo igual a la longitud del intervalo I. Esto es lo que llamamos "periodizar a f en el intervalo I". Naturalmente, si la función f ya es por sí misma periódica con periodo T y elegimos como intervalo I un intervalo cuya longitud sea igual a T, por ejemplo el intervalo [-T/2, T/2] o el [0, T], la serie de Fourier obtenida representará a la función en todo ℝ.
No olvides que puedes cambiar los valores de una función en un conjunto finito de puntos sin que ello afecte para nada a sus coeficientes de Fourier. Igualmente, si una función está definida en un intervalo excepto en un núnero finito de puntos del mismo puedes definir como quieras dicha función en esos puntos sin que ello afecte para nada a sus coeficientes de Fourier. La única condición que se precisa para calcular los coeficientes de Fourier de una función en un intervalo es que la función sea integrable en dicho intervalo y esa condición no se ve alterada si modificas la función en un conjunto finito de puntos. 
Si periodizas una función en un intervalo y la función tiene valores diferentes en los extremos del mismo, aparecen discontinuidades de salto las cuales, al hacer la representación gráfica, Mathematica une por un segmento vertical. En la inicialización de este cuaderno he modificado el comando Plot[ ] para que una por una línea discontinua las discontinuidades de salto de una función pero tienes que decirle tú a Mathematica dónde están dichas discontinuidades. La sintaxis del comando modificado es "Plot[funcion[x],{x,a,d1,d2,d3,...,b},opciones]" y representa la gráfica de "funcion[x]" en el intervalo [a,b] con discontinuidades de salto en los puntos d1,d2,d3,... . El nuevo comando solamente modifica a Plot[ ] cuando la lista {a,d1,d2,d3,...,b} tiene longitud mayor que 2, es decir, el comando Plot[ ] sigue funcionando igual que siempre si no pones ningún punto de discontinuidad. Aquí tienes un ejemplo (debes ejecutar la celda de inicialización).

f[
          x_] = (x - 3)/Abs[x - 3] + (x - 4)/Abs[x - 4] + (x - 7)/
            Abs[x - 7];  
Plot[f[x], {x, 0, 3, 4, 7, 9}, AspectRatio -> Automatic, 
        Ticks -> {{3, 4, 7}, Automatic}, 
        PlotStyle -> {Hue[ .65], Thickness[ .007]}];

Ejercicio 1

Plot[Mod[x, 1], {x, -3, -2, -1, 0, 1, 2, 3}, 
        AspectRatio -> Automatic, Ticks -> {Range[-3, 3], {1}}, 
        PlotStyle -> {Hue[0], Thickness[ .005]}, 
        ImageSize -> {435.25, 72.5}];

Plot[{Sin[π x], Sin[π Mod[x, 1]]}, {x, -3, 3}, 
        PlotStyle -> {Hue[0], {Thickness[ .007], Hue[ .65]}}, 
        AspectRatio -> Automatic];

Ejercicio 2

Ejercicio 3

Parte segunda: Coeficientes de Fourier. Aproximación por polinomios de              		                                                          		
		  Fourier.

Recuerda que la serie de Fourier de una función T-periódica integrable f es la serie de funciones

(Subscript[a, 0])/
        2 + (∑_(k ⩾ 1))((Subscript[a, k]) cos(
                2  π k t/T)  +  
            (Subscript[b, k]) sen(2 π k t/T))

donde los (Subscript[a, k])se llaman coeficientes coseno y los (Subscript[b, k]) coefiecientes seno de f y vienen dados por:

(Subscript[a, k])=(2/T)(∫_(-T/2)^(T/2))f(t) cos(
                    2  π k t/
                        T) dt     (k ≥ 
                      0),       (Subscript[b, k])=(2/T)(∫_(-T/2)^(T/2))f(t)sen(2  π k t/T)dt      (k ≥ 1)

Equivalentemente, podemos escribir dicha serie en la forma

\(TraditionalForm\`∑\+\(k ∈ ℤ\)c\_k\ e\^\(2  i\ π\ k\ t/T\)\)

donde (Subscript[c, k]) son los coeficientes de Fourier de f y vienen dados por:

(Subscript[c, k])=(1/T)(∫_(-T/2)^(T/2))\(\(f(
                    t)\) \(e\^\(\(-2\) i\ π\ k\ t/
                          T\)\) \(dt\)\(\ \ \ \ \ \)\((k ∈ ℤ)\)\(\ \ \ \)\)

La función

(Subscript[S, N])(t)= (∑_(k = -N)^N)(Subscript[c, k]) (e^(2  i π k t/T))=(Subscript[a, 0])/2+(∑_(k = 1)^N)((Subscript[a, k]) cos(2  π k t/T)  +  
                      (Subscript[b, k]) sen(2 π k t/T))

se llama el polinomio de Fourier de orden N de f.

Te propongo ahora, entre otras cosas, que calcules de forma simbólica los coeficientes de Fourier complejos así como los coeficientes coseno y seno de algunas funciones sencillas. Al hacer estos cálculos es frecuente que aparezcan algunas de las expresiones Sin[n π], Cos[n π], Exp[i π n], Exp[-i π n], (((-1))^n), donde n es un número entero. Conviene sustituir estas expresiones por sus correspondientes valores los cuales dependen de la paridad de n. En la celda siguiente he definido un comando que hace estas simplificaciones. La sintaxis del comando es "simplifica[expresion,variable]" y lo que hace es atribuir localmente a "variable" el carácter de número entero y aplicar a "expresion" (se supone que "expresion" depende de "variable") un conjunto de reglas de sustitución apropiadas para reemplazar las expresiones antes señaladas por sus valores.

Con Mathematica 5 podemos usar el comando "Simplify[expresion,condiciones]"; por ejemplo podemos decirle a Mathematica 5 que simplifique una expresión en la que interviene una variable k y le informamos de que dicha variable toma valores enteros positivos con Simplify[expresion[k], k ∈ Integers && k > 0]. Además en Mathematica 5 están implementadas las reglas de sustitución que usa el comando que se define en la siguiente celda; por eso, si trabajas con Mathematica 5 o posterior no es preciso que evalúes dicha celda (aunque si lo haces tampoco causará ningún problema).

simplifica[expr__, n_]:=Module[{par, impar, entero, reglas, st},impar[i_] := OddQ[i];par[i_] := EvenQ[i];  (* 
                cambiamos el nombre para modificar la función *) 
impar[i_  j_] := impar[i]  &&  impar[j];  
par[i_  j_] := par[i]  || par[j];
  
impar[
                i_  +  
                  j_] := (impar[i]  &&  par[j] ) ||  (par[i]  &&  
                    impar[j] );  
par[
                i_  +  
                  j_] := (par[i]  &&  par[j] ) ||  (impar[i]  &&  
                    impar[j] );
  entero[x_] := IntegerQ[x];entero[x_ y_] := entero[x] && entero[y];entero[x_  +  y_] := entero[x] && entero[y];
  
reglas = {Sin[ q_?entero  π] :> 0, 
                Cos[q_?entero  π] :> (((-1))^q), 
                Exp[I π q_?
                        entero] :> (((-1))^q), 
                Exp[-I π q_?
                        entero] :> (((-1))^q), 
                0^q_?entero :> 
                  0, ((-1)^m_?entero)^
                    q_?entero :> (-1)^(m q), (-1)^(q_?impar  m_) :>  (-1)^m, (-1)^
                    q_?par :> 
                  1, (-1)^
                    q_?impar :> -1, (-1)^
                    m_ :> (((-1))^m)};
st /:  entero[st] = True;  
Simplify[
                Simplify[expr /. {n -> st}] //. reglas] /. 
              st -> n  ]

En el caso de que hayas evaluado la celda anterior, para comprobar si el comando funciona correctamente ejecuta la siguiente celda.

b[n_] := ((4 Sin[((n π)/2)]^3 Cos[(((n + 1) π)/2)]^n)/((((-1))^(3  n)) +  n π))
simplifica[b[2  n], n]
simplifica[b[2  n - 1], n]

Si estás trabajando con Mathematica 5 ejecuta la celda siguiente y compara los resultados obtenidos.

Simplify[b[2  n], n ∈ Integers]  
Simplify[b[2  n - 1], n ∈ Integers]  
  
simplifica[
      b[n], n]
Simplify[b[n], n ∈ Integers]  
Clear[
      b]

Te aconsejo que seas precavido al usar el comando "simplifica[ ]", pues no estoy totalmente seguro de que sea infalible. Por ello debes comprobar  cualquier resultado que no te convenza.

En los ejercicios siguientes debes calcular, unas veces de forma simbólica y otras de forma numérica, los coeficientes de Fourier de unas cuantas funciones y hacer algunas animaciones gráficas para comprobar cómo los polinomios de Fourier se van aproximando a la función. Además deberás representar en cada caso los espectros de amplitudes y de fases.

Para representar los espectros de amplitudes y de fases debes recordar la relación entre los coefiecientes de Fourier y los coeficientes seno y coseno que viene dada por:

(Subscript[c, 0]) = ((Subscript[a, 0])/2)\(\(,\)\(\ \)\(c\_n = \(a\_n\  - \ i\ b\_n\)\/2\)\(,\)\(\ \)\(c\_\(-n\) = \(a\_n\  + \ i\ b\_n\)\/2\)\), n = 1,2,...

Recuerda que el espectro de una señal cuya serie de Fourier compleja es \(TraditionalForm\`∑\+\(k ∈ ℤ\)c\_k\ e\^\(2  i\ π\ k\ t/T\)\) se define como el conjunto de puntos {(n/T, (Subscript[c, n])) : 
          n ∈ 𝕫}. El espectro de amplitudes es el conjunto  {(n/T, Abs[(Subscript[c, n])]) : 
          n ∈ 𝕫} y el espectro de fases es el conjunto {(n/T, Arg[(Subscript[c, n])]) : 
          n ∈ 𝕫} dichos conjuntos suelen representarse por segmentos de línea que unen los puntos (n/T, 0) con los puntos del espectro correspondiente. Deben quedar, más o menos, como los ves aquí debajo.

Sugerencias. Usa en lo que sigue las notaciones a[k], b[k],c[k] para los coeficientes respectivos. Aveces no es posible sustituir k por 0 (y puede que para algunos otros valores) para obtener los valores a[0],b[0] o c[0]. Cuando esto ocurra Mathematica dará un mensaje de error diciendo que ha encontrado la expresión 1/0. En tales casos lo mejor es definir directamente esos valores (a[0]=A[f,0,{a,b}], a[1]=A[f,1,{a,b}],b[1]=B[f,1,{a,b}]...). No olvides que b[0]=0 por lo que siempre puedes definirlo así. Para calcular los coeficientes de una señal deberás definir funciones a[k_], b[k_] o c[k_]. Puede ser conveniente usar la asignación directa en las definiciones, pues de esta forma tenemos calculados todos los coeficientes de golpe con una sola integración. Si usas la asignación diferida, cada vez que llames a los coeficientes deberá realizarse la integración correspondiente. Naturalmente, esto se refiere a los cálculos simbólicos, en los cálculos numéricos debes usar siempre la asignación diferida. Cuando pases de un ejercicio al siguiente no olvides limpiar todo lo que pueda molestar; al menos un "Clear[a,b,c]" parece necesario.

Ejercicio 4

Ejercicio 5

Ejercicio 6

Ejercicio 7

Ejercicio 8

Ejercicio 9