Revisión del 19:38 13 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 periodicas 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 afecta de ninguna forma 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 cierta funciones. Por otro lado acabamos de ver cómo éstas 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[/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[/math] y así:

Apliquemos lo anterior a la función [math]f(x)=x\cdot (1-x)[/math] en el compacto [math][0,1][/math], cuya gráfica y su correspondiente extensión impar se muestran a continuación.

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

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.

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.

TENEMOS QUE ESTIMAR A QUÉ FUNCIÓN SE PARECE

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)

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]f(x)=1_{x \leq \frac{1}{2}}(x)[/math] definida en el compacto [math][0,1][/math], cuya gráfica y su correspondiente extensión par se muestran a continuación.

Gráfico de la función [math] f(x)=1_{x\leq 1/2}(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)=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

Extensión par de la función [math] f(x)=1_{x\leq 1/2}(x), D(f)=[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]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)=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 no sólo cómo según aumenta el número de términos, la suma aproxima mejor la función, sino que desaparece el fenómeno de Gibbs en las discontinuidades de la función original.

Gráficas de las funciones [math]f(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')

6 Cambio de intervalo

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

7 Base trigonométrica compleja

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