Revisión actual del 09:33 15 feb 2024

1 Introducción

A lo largo de la Historia los matemáticos se han encontrado con problemas que [math]\textit{a priori}[/math] son imposibles de manejar, por ejemplo, debido a su generalidad. Sin embargo, esta es sin duda un arma de doble filo. Por un lado, el planteamiento de problemas tan generales pueden parecer inabarcables debido a la infinidad de posibles casos, pero por otro, su solución nos proporciona una poderosa herramienta que nos permite simplificar multitud de problemas que nos parecían en un principio intratables. Un ejemplo claro son las series de Fourier, que nos permiten aproximar de forma sencilla y muy precisa infinitud de funciones mediante una serie infinita.

2 Definición

Sean el espacio de Hilbert [math]L^2([-\pi,\pi])[/math] y una base de Hilbert de dicho espacio [math]\{ \frac{1}{2},\cos(nx),\sin(nx)\}_{n \in \mathbb{N}}[/math]. Sea [math] f \in L^2 ([-\pi,\pi]) [/math] entonces se llama serie de Fourier a la serie convergente tal que:

Los términos [math] a_0,a_n, b_n \in \mathbb{R} [/math] reciben el nombre de coeficientes de Fourier y se calculan como sigue:

Nótese que la familia [math]\{ \frac{1}{2},\cos(nx),\sin(nx)\}_{n \in \mathbb{N}}[/math] es una base de Hilbert salvo reescalamiento; no están normalizadas. Pese a esto, la teoría sigue siendo válida pues simplemente hay que dividir los coeficientes de Fourier por sus respectivas normas.

Otra cosa importante a recalcar es que la definición anterior se ha hecho para funciones definidas sobre el compacto [math][-\pi,\pi][/math]. Sin embargo, esto se puede definir para cualquier otro intervalo de la forma [math][-T,T][/math] mediante el cambio de variable [math]y=x\frac{T}{\pi}[/math], obteniendo la base:

A partir de esto también es posible generalizarlo a intervalos de la forma [math] [a,b] [/math]. Bastaría utilizar la base del espacio [math]L^2([-\frac{b-a}{2},\frac{b-a}{2}])[/math], ya que estaría conformada por funciones periódicas de periodo [math]b-a[/math], que coincidiría con la longitud del intervalo de definición de la función que queremos aproximar. De igual forma podemos trabajar con intervalos abiertos en los dos o en alguno de los extremos, pues un punto es un conjunto de medida nula y al trabajar con integrales esto no afecta al resultado final.

3 Base trigonométrica

En los siguientes apartados se mostrarán una serie de gráficas junto con sus respectivos códigos de MatLab.

Para entender el significado de la definición anterior, veamos cómo se comportan los elementos de la base en el compacto [math][-1,1][/math]. Por lo que hemos dicho, esta base resulta:

A continuación se muestra una gráfica de los diez primeros elementos de la base de la forma [math]\cos(n\pi x)[/math].

Término de la base [math] \cos(n\pi x)[/math]

close all
clear all
N=10;
xx=-1:10^(-3):1;
yy=zeros(1,length(xx));
for n=[1:N]
    e=@(x) (cos(pi*n*x));
    for i=1:length(xx)
        yy(i)=e(xx(i));
    end
    subplot(2,5,n);
    xline(0);
    hold on
    yline(0);
    plot(xx,yy,'b')
    axis equal
    xlim([-1,1])
    ylim([-1,1])
    title("n="+n,'Interpreter','latex')
    hold off
end

En la figura anterior se aprecia cómo según aumenta el valor de [math] n [/math] el periodo de la función va disminuyendo. Esto también se puede ver fácilmente a partir de la expresión analítica:

El razonamiento con los elementos de la forma [math]\sin(n\pi x)[/math] es análogo. A continuación se muestran las gráficas de dichos elementos con [math]n=1,...,10 [/math].

Término de la base [math] \sin(n\pi x)[/math]

close all
clear all
N=10;
xx=-1:10^(-3):1;
yy=zeros(1,length(xx));
for n=[1:N]
    e=@(x) (sin(pi*n*x));
    for i=1:length(xx)
        yy(i)=e(xx(i));
    end
    subplot(2,5,n);
    xline(0);
    hold on
    yline(0);
    plot(xx,yy,'b')
    axis equal
    xlim([-1,1])
    ylim([-1,1])
    title("n="+n,'Interpreter','latex')
    hold off
end

Por último, tenemos el elemento constante de la base [math] f(x)=\frac{1}{2}[/math], cuya gráfica es la que sigue.

Término de la base [math] 1/2[/math]

close all
clear all
xx=-1:10^(-3):1;
yy=zeros(1,length(xx));
e=@(x) (1/2);
for i=1:length(xx)
    yy(i)=e(xx(i));
end
xline(0);
hold on
yline(0);
plot(xx,yy,'b','LineWidth',1.5)
axis equal
xlim([-1,1])
ylim([-1,1])
hold off

Habíamos visto en el apartado anterior que la serie de Fourier se define como una combinación lineal infinita de unas ciertas funciones. Por otro lado acabamos de ver cómo estas son funciones periódicas cuya frecuencia aumenta según lo hace [math]n[/math], junto a una constante. Es decir, que en el fondo lo que estamos haciendo es descomponer una función dada como una suma de ondas de diferentes frecuencias, como se verá más adelante.

4 Aproximación de una función continua

Veamos ahora cómo utilizar la serie de Fourier para aproximar funciones continuas.

Sean [math]g(x)[/math] una función continua en el intervalo [math][0,b][/math] tal que [math]g(0)=0[/math] y [math]f[/math] su extensión impar:

Por construcción esta función es continua en [math][-b, b] \setminus \{0\}[/math]. Por otro lado, es fácil verificar la continuidad en [math]0[/math] mediante límites laterales. Por tanto, la extensión de [math]g[/math] que hemos definido es continua en [math][-b,b][/math] y su aproximación mediante la serie de Fourier converge puntualmente en [math](-b,b)[/math]. En cuanto a los extremos sabemos que debe converger al punto medio de las imágenes de los extremos, es decir, a [math]0[/math]:

Una vez hecho esto sólo queda calcular los coeficientes de Fourier para la base correspondiente a este intervalo, que siguiendo lo explicado en la definición quedaría [math]\lbrace \frac{1}{2},\cos(\frac{n\pi}{b} x),\sin(\frac{n\pi}{b} x)\rbrace_{n \in \mathbb{N}}[/math]. Así, los coeficientes a calcular son:

Es aquí donde podemos ver la ventaja de haber definido la extensión impar de la función dada, pues los coeficientes [math]a_0[/math] y [math]\{b_n\}_{n \in \mathbb{N}}[/math] al ser integrales de funciones impares son nulos. Con lo que concluimos que para calcular la serie de Fourier basta con hallar los coeficientes [math]\{a_n\}_{n \in \mathbb{N}}[/math] y así:

Apliquemos lo anterior a la función [math]f(x)=x(1-x)[/math] en el compacto [math][0,1][/math], a continuación se muestra la definición de la función [math]f[/math] y su correspondiente gráfica.

function [y] = f(x)
    y=(x.*(1-x));
end

Gráfico de la función [math] f(x)=x\cdot (1-x), D(f)=[0,1][/math]

close all
clear all

xx=0:10^(-3):1;
yy=zeros(1,length(xx));
for i=1:length(xx)
    yy(i)=f(xx(i));
end
plot(xx,yy,'r','LineWidth',1.5)
hold on
xline(0);
yline(0);
axis equal
xlim([-1,1])
ylim([-1,1])
hold off

Veamos ahora la definición de su extensión impar y su gráfica. La función queda definida en el siguiente código.

function [y] = f(x)
    if x>=0
        y=(x.*(1-x));
    else
        y=-(-x.*(1+x));
    end
end

Extensión impar de la función [math] f(x)=x\cdot (1-x), D(f)=[0,1][/math]

clear all
xx=0:10^(-3):1;
zz=-1:10^(-3):0;
yy=zeros(1,length(xx));
ww=zeros(1,length(zz));
for i=1:length(xx)
    yy(i)=f(xx(i));
end
for i=1:length(xx)
    ww(i)=f(zz(i));
end
plot(xx,yy,'r','LineWidth',1.5)
hold on
plot(zz,ww,'y','LineWidth',1.5)
xline(0);
yline(0);
axis equal
xlim([-1,1])
ylim([-1,1])
hold off

Para ver cómo converge la serie de Fourier se han calculado los [math]n[/math] primeros términos para [math]n=1,5,10 [/math] y se han dibujado en la misma gráfica junto a la función a aproximar. A lo largo del artículo se denotará como [math]f_n[/math] a la suma de los [math]n[/math] primeros términos de la serie de Fourier.

Gráficas de las funciones [math]f(x)[/math] y [math]f_{n}(x)[/math] para [math]n=1,5,10[/math]

close all
clear all
xx=-1:10^(-3):1;
yy=zeros(1,length(xx));
for i=1:length(xx)
    yy(i)=f(xx(i));
end
l=1;
puntos=[1,5,10];
for n=puntos
    subplot(1,length(puntos),l)
    plot(xx,yy,'r','LineWidth',1.5)
    hold on
    axis equal
    a=zeros(1,n);
    fn=zeros(1,length(xx));
    for k=1:n
        a(k)=trapz(xx,yy.*sin(k.*pi.*xx));
        fn=fn+a(k).*sin(k.*pi.*xx);
    end
    plot(xx,fn,'b')
    xline(0);
    yline(0);
    axis equal
    xlim([-1,1])
    ylim([-1,1])
    legend({'$f(x)$','$f_{n}(x)$'},'Interpreter','latex','Location','southeast')
    title('n= '+string(n),'Interpreter','latex')
    hold off
    l=l+1;
end

En las gráficas anteriores no se aprecia como mejora la aproximación según aumenta el valor de [math]n[/math], sin embargo, esto se aprecia muy bien en las que se muestran a continuación. En la primera, se observa cómo disminuye rápidamente el error en la norma de [math]L^2(-1,1)[/math] según aumentamos el valor de [math]n[/math]. La segunda gráfica representa la misma idea pero empleando la norma uniforme.

Gráficos de errores

close all
clear all
xx=-1:10^(-3):1;
yy=zeros(1,length(xx));
for i=1:length(xx)
    yy(i)=f(xx(i));
end
puntos=[1:50];
errl2=zeros(1,length(puntos));
erruniforme=zeros(1,length(puntos));
for n=puntos
    a=zeros(1,n);
    fn=zeros(1,length(xx));
    for k=1:n
        a(k)=trapz(xx,yy.*sin(k.*pi.*xx));
        fn=fn+a(k).*sin(k.*pi.*xx);
    end
    errl2(n)=trapz(xx,(yy-fn).^2).^(1/2);
    erruniforme(n)=max(abs(yy-fn));
end
subplot(2,1,1)
plot(puntos,errl2,'LineWidth',1.5)
subplot(2,1,2)
plot(puntos,erruniforme,'LineWidth',1.5)

Se aprecia un comportamiento similar a la función [math] f(x)=\frac{1}{x^2}[/math].

5 Aproximación de una función discontinua

Una vez hemos visto el ejemplo de una aproximación de una función continua tras extenderla de forma impar, veamos cómo hacerlo si la extensión es par.

Sean [math]g(x)[/math] una función definida en el intervalo [math][0,b][/math] y [math]f[/math] su extensión par:

Al igual que con la extensión impar, una vez hecho esto sólo queda calcular los coeficientes de Fourier para la base correspondiente a este intervalo, que siguiendo lo explicado en la definición quedaría [math]\lbrace \frac{1}{2},\cos(\frac{n\pi}{b} x),\sin(\frac{n\pi}{b} x)\rbrace_{n \in \mathbb{N}}[/math]. Por tanto, los coeficientes a calcular son:

Es aquí donde podemos ver que la idea detrás de ambas extensiones es en el fondo la misma, forzar a que una familia de coeficientes sean todos nulos; en este caso [math]\{a_n\}_{n \in \mathbb{N}}[/math]. Al igual que antes esto se debe a que son integrales de funciones impares. Con lo que concluimos que para calcular la serie de Fourier basta con hallar los coeficientes [math]a_0[/math] y [math]\{b_n\}_{n \in \mathbb{N}}[/math]. Así:

Apliquemos lo anterior a la función [math]g(x)=1_{x \leq \frac{1}{2}}(x)[/math] definida en el compacto [math][0,1][/math], cuya definición y gráfica se muestra a continuación.

function [y] = g(x)
    if x>1/2
        y=0;
    else
        y=1;
    end
end

Gráfico de la función [math] g(x)=1_{x\leq 1/2}(x), D(g)=[0,1][/math]

close all
clear all
xx=0:10^(-3):1;
yy=zeros(1,length(xx));
for i=1:length(xx)
    yy(i)=g(xx(i));
end
plot(xx,yy,'r','LineWidth',1.5)
hold on
xline(0);
yline(0);
axis equal
xlim([-1,1])
ylim([0,1.5])
hold off

Por otro lado, dibujamos la gráfica de la extensión par de [math]g(x)[/math]:

function [y] = g(x)
    if abs(x)>1/2
        y=0;
    else
        y=1;
    end
end

Extensión par de la función [math] g(x)=1_{x\leq 1/2}(x), D(g)=[0,1][/math]

close all
clear all
xx=0:10^(-3):1;
zz=-1:10^(-3):0;
yy=zeros(1,length(xx));
ww=zeros(1,length(zz));
for i=1:length(xx)
    yy(i)=g(xx(i));
end
for i=1:length(xx)
    ww(i)=g(zz(i));
end
plot(xx,yy,'r','LineWidth',1.5)
hold on
plot(zz,ww,'y','LineWidth',1.5)
xline(0);
yline(0);
axis equal
xlim([-1,1])
ylim([0,1.5])
hold off

Es fácil ver que esta función es continua en todo su dominio salvo en [math]x=\frac{1}{2}[/math] y [math]x=-\frac{1}{2}[/math], por lo que la serie convergerá puntualmente en todo el compacto salvo en esos puntos; en los que valdrá [math]\frac{1}{2}[/math]. Esto se aprecia muy bien en la siguiente figura, en la que se muestra la serie de Fourier con [math]n=1,5,10[/math].

Gráficas de las funciones [math]g(x)[/math] y [math]f_{n}(x)[/math] para [math]n=1,5,10[/math]

close all
clear all
xx=-1:10^(-3):1;
yy=zeros(1,length(xx));
for i=1:length(xx)
    yy(i)=g(xx(i));
end
l=1;
puntos=[1,5,10];
for n=puntos
    subplot(1,length(puntos),l)
    plot(xx,yy,'r','LineWidth',1.5)
    hold on
    axis equal
    b=zeros(1,n);
    fn=1/2*ones(1,length(xx));
    for k=1:n
        b(k)=trapz(xx,yy.*cos(k.*pi.*xx));
        fn=fn+b(k).*cos(k.*pi.*xx);
    end
    plot(xx,fn,'b')
    xline(0);
    yline(0);
    axis equal
    xlim([-1,1])
    ylim([0,1.5])
    legend({'$f(x)$','$f_{n}(x)$'},'Interpreter','latex','Location','south')
    title('n= '+string(n),'Interpreter','latex')
    hold off
    l=l+1;
end

El hecho de que la serie converja puntualmente en todo punto salvo en los de discontinuidad da lugar a la aparición de oscilaciones cada vez más rápidas según aumentamos el valor de [math]n[/math]. Este fenómeno recibe el nombre de fenómeno de Gibbs. Para aliviar este efecto se emplean las sumas de Cesàro que se definen como:

A continuación se muestran dichas sumas junto con la función que buscábamos aproximar para [math]N=1,5,10[/math]. En estas gráficas se observa cómo la suma aproxima mejor la función según aumenta el número de términos, y además, cómo que desaparece el fenómeno de Gibbs en las discontinuidades de la función original.

Gráficas de las funciones [math]g(x)[/math] y [math]S_{N}(x)[/math] para distintos [math]N=1,5,10[/math]

close all
clear all
xx=-1:10^(-3):1;
yy=zeros(1,length(xx));
for i=1:length(xx)
    yy(i)=g(xx(i));
end
l=1;
puntos=[1,5,10];
for n=puntos
    fN=zeros(n,length(xx));
    subplot(1,length(puntos),l)
    plot(xx,yy,'r','LineWidth',1.5)
    hold on
    axis equal
    b=zeros(1,n);
    fn=1/2*ones(1,length(xx));
    for k=1:n
        b(k)=trapz(xx,yy.*cos(k.*pi.*xx));
        fn=fn+b(k).*cos(k.*pi.*xx);
        fN(k,:)=fn;
    end
    SN=zeros(n,length(xx));
    for i=[1:n]
        SN(i,:)=ones(1,i)*fN(1:i,:)/(i+1);
    end
    plot(xx,SN(n,:),'b')
    xline(0);
    yline(0);
    axis equal
    xlim([-1,1])
    ylim([0,1.5])
    legend({'$f(x)$','$S_{N}(x)$'},'Interpreter','latex','Location','south')
    title('N= '+string(n),'Interpreter','latex')
    hold off
    l=l+1;
end

Al igual que el apartado anterior, concluimos calculando el error cometido en la aproximación. En este caso tomaremos como aproximación la obtenida mediante las sumas de Cesàro.

Gráficos de errores

close all
clear all
xx=-1:10^(-3):1;
yy=zeros(1,length(xx));
for i=1:length(xx)
    yy(i)=g(xx(i));
end
plot(xx,yy)
hold on
puntos=[1:200];
errl2=zeros(1,length(puntos));
erruniforme=zeros(1,length(puntos));
fN=zeros(length(puntos),length(xx));
for n=puntos
    b=zeros(1,n);
    fn=1/2*ones(1,length(xx)); 
    for k=1:n
        b(k)=trapz(xx,yy.*cos(k.*pi.*xx));
        fn=fn+b(k).*cos(k.*pi.*xx);
    end
    errl2(n)=trapz(xx,(yy-fn).^2).^(1/2);
    erruniforme(n)=max(abs(yy-fn));
    fN(n,:)=fn;
end
SN=zeros(length(puntos),length(xx));
errl2SN=zeros(1,length(puntos));
erruniformeSN=zeros(1,length(puntos));
for i=puntos
    SN(i,:)=ones(1,i)*fN(1:i,:)/(i+1);
    errl2SN(i)=trapz(xx,(yy-SN(i,:)).^2).^(1/2);
    erruniformeSN(i)=max(abs(yy-SN(i,:)));
end
subplot(2,1,1)
plot(puntos,errl2,'LineWidth',1.5)
xlabel('$n$','Interpreter','latex')
subplot(2,1,2)
plot(puntos,erruniforme,'LineWidth',1.5)
xlabel('$n$','Interpreter','latex')

Llama la atención que la disminución del error en la norma del supremo es lineal a diferencia del anterior ejemplo. En la norma de [math]L^2([-1,1])[/math] se comporta de manera similar, pero más suave como consecuencia de la eliminación de las oscilaciones.

6 Cambio de intervalo

Hasta ahora se han visto tan solo series de Fourier asociadas a funciones definidas en intervalos simétricos con respecto al [math] 0 [/math] o reducibles a serlo, mediante extensiones pares o impares. De cara a formular desarrollos de Fourier en un intervalo genérico, se suele usar una estrategia que consiste en extender nuestra función por periodicidad, y tomar el tramo que contenga al [math] 0 [/math] para el desarrollo. Una vez hecho esto, se extiende la serie por periodicidad de la misma manera para quedarnos con el tramo que contiene al intervalo de definición de la función inicial.

Pongamos como ejemplo la siguiente función

Podemos extenderla por periodicidad con periodo [math]3-1=2[/math] de acuerdo a lo mencionado en la definición. Así, [math]f(x)=f(x+2k),\thinspace k\in \mathbb{Z}[/math].

Una vez visto esto, basta utilizar la base de [math]L^{2}\left( -1,1 \right)[/math]

y extender la serie por periodicidad hasta que [math][1,3]=Dom(f)[/math].

Gráficas de las funciones [math]f(x)[/math] y [math]f_{n}(x)[/math] para distintos [math]n[/math]

close all
clear all
xx=-1:10^(-3):1;
haux=@(x)(x.*exp(-x));
h=@(x)((x+2).*exp(-x-2));
yy=zeros(1,length(xx));
for i=1:length(xx)
    yy(i)=h(xx(i));
end
puntos=[5,10,20];
i=1;
for n=puntos
    a=zeros(1,n);
    b=zeros(1,n);
    fn=trapz(xx,yy)*ones(1,length(xx))/(2);
    for k=1:n
        a(k)=trapz(xx,yy.*sin(k.*pi.*xx));
        b(k)=trapz(xx,yy.*cos(k.*pi.*xx));
        fn=fn+a(k).*sin(k.*pi.*xx)+b(k).*cos(k.*pi.*xx);
    end
    subplot(1,3,i)
    fplot(haux,'r','LineWidth',1.5)
    hold on
    plot(xx+2*ones(1,length(xx)),fn,'b')
    xline(0);
    yline(0);
    axis equal
    xlim([1,3])
    ylim([-1,1])
    legend({'$f(x)$','$f_{n}(x)$'},'Interpreter','latex','Location','southeast')
    title('n= '+string(n),'Interpreter','latex')
    hold off
    hold off
    i=i+1;
end

Nótese que como estamos aproximando una función que da valores distintos a los extremos del intervalo de definición, podemos observar de nuevo el fenómeno de Gibbs.

7 Base trigonométrica compleja

Otra formulación de la base trigonométrica para formar series de Fourier es la exponencial compleja. En el espacio [math] L^{2}(-\pi,\pi) [/math], tendríamos la siguiente base de funciones:

Los coeficientes de la serie por tanto para una función [math]f\in L^{2}(-\pi,\pi)[/math] se obtendrían mediante la fórmula

y con esto se tendría

Pongamos como ejemplo realizar el desarrollo de la siguiente función:

Por lo visto anteriormente, basta realizar la aproximación de la función en la base de [math] L^{2}(-\frac{1}{2},\frac{1}{2}) [/math]. Realizando el cambio de variable pertinente, se obtiene que la base es [math]\lbrace e^{i2\pi nx}\rbrace_{n\in\mathbb{Z}}[/math] y se obtendrían los coeficientes con las integrales mencionadas anteriormente. Nótese que ahora las funciones de la base no solo son ortogonales, también son ortonormales.

Gráficas de las funciones [math]f(x)[/math] y [math]f_{n}(x)[/math] para distintos [math]n[/math]

close all
clear all
xx=-1/2:10^(-3):1/2;
haux=@(x) (4*x.*(1/2-x).^2);
h=@(x)(4*(x+1/2).*(1/2-x-1/2).^2);
yy=zeros(1,length(xx));
for n=1:length(xx)
    yy(n)=h(xx(n));
end
l=1;
puntos=[5,10,20];
for n=puntos
    a=zeros(1,n);
    b=zeros(1,n);
    fn=trapz(xx,yy).*ones(1,length(xx));
    for k=1:n
        a(k)=trapz(xx,yy.*conj(exp(1i.*2*k.*pi.*xx)));
        b(k)=trapz(xx,yy.*conj(exp(1i.*2*(-k).*pi.*xx)));

        fn=fn+a(k).*exp(1i.*2*k.*pi.*xx)+b(k).*exp(1i.*2*(-k).*pi.*xx);
    end
    subplot(1,3,l)
    plot(xx+ones(1,length(xx))/2,haux(xx+ones(1,length(xx))/2),'r','LineWidth',1.5)
    hold on
    plot(xx+ones(1,length(xx))/2,fn,'b')
    xline(0);
    yline(0);
    axis equal
    xlim([0,1])
    ylim([-0.25,1.25])
    legend({'$f(x)$','$f_{n}(x)$'},'Interpreter','latex','Location','northwest')
    title('n= '+string(n),'Interpreter','latex')
    hold off
    l=l+1;
end

7.1 Referencias