Series de Fourier (Raúl, Sofía, Jaime)
| Trabajo realizado por estudiantes | |
|---|---|
| Título | Series de Fourier |
| Asignatura | EDP |
| Curso | 2023-24 |
| Autores | Raúl Ortega
Sofía Gómez Jaime Sáenz de Miera |
| Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura | |
En este artículo vamos a estudiar qué son las series de Fourier utilizando varias funciones que ejemplifiquen todas las situaciones que nos podemos encontrar al calcular una serie de Fourier.
Contenido
1 Preliminares
Antes de definir series de Fourier, necesitamos conocer el espacio [math] L^2= \{ f: \Omega \longrightarrow \mathbb{R} : \int_{\Omega}|f(x)|^2 \, dx \lt \infty \} [/math], que es un espacio de Hilbert.
Esto nos permite definir el producto escalar como [math]\langle f,g\rangle_{L^2(\Omega)}= \int_{\Omega} f(x) \cdot g(x) \, dx [/math], con [math] f,g\in L^2(\Omega) [/math], y el módulo como [math] \| \cdot \| _ {L^2 (\Omega)} = \sqrt{\langle\cdot , \cdot \rangle_{L^2(\Omega)}} [/math].
2 Series de Fourier
El conjunto [math] \{ \frac{1}{2},\cos(nx),\sin(nx) \}_{n\in\mathbb{N}}[/math] es la base trigonométrica del espacio [math] L^2([-\pi,\pi])[/math], por tanto, toda función de [math] L^2([-\pi,\pi])[/math] puede expresarse como una combinación lineal de los elementos de esta. Además, podemos observar que las funciones de la base son [math] 2\pi[/math] periódicas, y es fácil comprobar que es una base ortogonal en [math] L^2([-\pi,\pi])[/math].
Para entender mejor la base vamos a graficar los primeros 10 términos de esta.
x=linspace(-1,1,200);
y=1/2*ones(1,200);
hold on
plot(x,y)
for n=1:4
plot(x,cos(n*pi*x))
plot(x,sin(n*pi*x))
end
plot(x,cos(5*pi*x))
title('Primeros 10 términos de la base trigonométrica');
xlabel('x');
ylabel('y');
legend('1/2', '$\cos(\pi x)$', '$\sin(\pi x)$', '$\cos(2\pi x)$', '$\sin(2\pi x)$', '$\cos(3\pi x)$', '$\sin(3\pi x)$', '$\cos(4\pi x)$', '$\sin(4\pi x)$','$\cos(5\pi x)$', 'Interpreter','latex');
Una vez hemos entendido lo que es una base trigonométrica, podemos definir la serie de Fourier de una función [math] f\in L^2([-\pi,\pi])[/math] como: [math] f \approx \frac{a_0}{2}+\sum^{N}_{n=1}(a_n\cos{(nx)}+b_n\sin{(nx)})[/math]
Como la base trigonométrica es ortogonal, podemos obtener los coeficientes de la serie de Fourier como:
[math] a_0=\frac{\langle f,\frac{1}{2}\rangle_{L^2([-\pi,\pi])}}{\| \frac{1}{2} \|^2 _ {L^2 ([-\pi,\pi])}} =\frac{1}{\pi}\int^{\pi}_{-\pi} f(x)\, dx [/math],
[math] a_n=\frac{\langle f,\cos(nx)\rangle_{L^2([-\pi,\pi])}}{\| \cos(nx) \|^2 _ {L^2 ([-\pi,\pi])}} =\frac{1}{\pi}\int^{\pi}_{-\pi} f(x)\cos(x)\, dx [/math].
[math] b_n=\frac{\langle f,\sin(nx)\rangle_{L^2([-\pi,\pi])}}{\| \sin(nx) \|^2 _ {L^2 ([-\pi,\pi])}} =\frac{1}{\pi}\int^{\pi}_{-\pi} f(x)\sin(x)\, dx [/math],
3 Extensión de funciones
Para simplificar los cálculos de la serie de Fourier podemos extender la función de forma par o impar, lo cual hace que algunos términos se cancelen. Si la función es par, los coeficientes [math] b_n[/math] serán cero, ya que al ser el seno una función impar, el producto [math] f(x)\sin(x)[/math] también lo va a ser, y la integral de una función impar sobre un intervalo simétrico se anula. Análogamente, si la función es impar los coeficientes [math] a_0,a_n,n\geq 1[/math] van a ser cero.
Para ejemplificar esto, tomaremos la función [math]f(x)=x(1-x)[/math] en el intervalo [math][0,1][/math]. La extenderemos de forma impar en el intervalo [math][-1,1][/math] y la aproximaremos usando su desarrollo en serie de Fourier.
Para hacer esto, primero vamos a comprobar que la extensión sigue siendo una función continua, esto se debe a que la original pasaba por el origen de coordenadas.
x1=0:10^(-3):1;
x2=-1:10^(-3):0;
figure(1)
hold on
plot(x1,x1.*(1-x1),'r',LineWidth=1.5)
plot(x2,x2.*(1+x2),'b',LineWidth=1.5)
title('Extension impar de $f(x)=x(1-x)$', 'interpreter','latex',FontSize=15)
xlabel('x')
ylabel('y')
Ahora calcularemos la serie de Fourier aproximando los coeficientes de la serie de forma numérica mediante el calculo de las integrales por el método del trapecio.
a=0; b=1;
u_1=a:10^(-3):b;
N=length(u_1)-1;
h=(b-a)/N;
w=ones(N+1,1);
w(1)=1/2; w(N+1)=1/2;
figure(2)
n=[1,5,10];
f_n=zeros(length(n),N+1);
for i=1:length(n)
for k=1:n(i)
g=(u_1.*(1-u_1).*sin(k*pi*u_1))';
a_k=2*h*w'*g;
f_n(i,:)=f_n(i,:)+a_k*sin(k*pi*u_1);
end
subplot(1,3,i)
hold on
plot(u_1,f_n(i,:),'r',LineWidth=1.5)
plot(u_1,u_1.*(1-u_1),'--b',LineWidth=1.5)
title('n='+string(n(i)))
legend('$f_n(x)$','$f(x)$','interpreter','latex','Location','northwest')
grid('on')
endEn la gráfica se muestra la aproximación de [math] f(x) [/math] por los primeros términos de su serie de Fourier (1,5 y 10 respectivamente). Como ya de partida la función [math] f(x) [/math] tenía una forma bastante sinusoidal en el intervalo [math] [0,1] [/math], vemos que la aproximación es bastante buena con tan sólo unos pocos términos de la serie.
Aunque a simple vista se ve que la aproximación es buena, vamos a comprobarlo utilizando los errores en la norma [math] L^2 [/math] y en la uniforme.
a=0; b=1;
u_1=a:10^(-3):b;
u_2=-1:10^(-3):0;
N=length(u_1)-1;
h=(b-a)/N;
w=ones(N+1,1);
w(1)=1/2; w(N+1)=1/2;
for n=1:10
f_n=zeros(1,N+1);
for k=1:n
g=(u_1.*(1-u_1).*sin(k*pi*u_1))';
a_k=2*h*w'*g;
f_n=f_n+a_k*sin(k*pi*u_1);
end
g1=abs(u_1.*(1-u_1)-f_n).^2;
error1(n)=(h*w'*g1')^(1/2);
error2(n)=max(abs(u_1.*(1-u_1)-f_n));
end
nn=1:1:10;
hold on
plot(nn,error1,'b',LineWidth=1.5)
plot(nn,error2,'r',LineWidth=1.5)
legend('$L^2$','Supremo','interpreter','latex')
Como podemos observar en la gráfica, cuanto mayor es la [math] n [/math], mejor aproxima la función y además, podemos ver que el error en [math] n=10 [/math] es casi 0. También, podemos observar que el error de la norma [math] L^2 [/math] es menor que el de la norma uniforme.
Ahora vamos a tomar la función [math] g(x)=1_{x\leq 1/2}(x) [/math] extendiéndola de forma par al intervalo [math] [-1,1] [/math] y aproximándola por su serie de Fourier.
x1=[0:0.01:1];
x2=[-1:0.01:0];
y1=[ones(50,1);zeros(51,1)];
y2=[zeros(50,1);ones(51,1)];
figure(1)
plot(x1,y1,LineWidth=1.5,Color='r')
ylim([-1,2])
hold on
plot(x2,y2,LineWidth=1.5,Color='b')
title('Extension par de $g(x)=1_{x\leq 1/2}(x)$','Interpreter','latex',FontSize=15)
xlabel('x')
ylabel('y')
Comprobamos que como la función era discontinua, la extensión sigue siendo discontinua, y por tanto vamos a ver que la serie de Fourier que la aproxime......
x=[0:0.01:1];
y=[ones(50,1);zeros(51,1)];
u=0:10^(-3):1;
N=length(u)-1;
w = [ones(N/2,1);zeros(N/2+1,1)];
h = 10^(-3);
nn=[1,5,10];
for j=1:length(nn)
for n=1:nn(j)
f_n=0.5*ones(1,N+1);
for k=1:n
g=(cos(k*pi*u))';
a_k=2*h*w'*g;
f_n=f_n+a_k*cos(k*pi*u);
end
end
subplot(1,3,j)
plot(x,y,'b--',LineWidth=1.5)
hold on
plot(u,f_n,'r',LineWidth=1.5)
title('n='+string(nn(j)))
legend('f(x)','$f_n(x)$','interpreter','latex')
endEn la gráfica se muestra la aproximación de [math] g(x) [/math] por los primeros términos de su serie de Fourier (1,5 y 10 respectivamente).
También se observa que cerca del punto de discontinuidad ([math] x=1/2 [/math]) los primeros términos de la serie "fallan" en aproximar bien a la función, ya que se producen muchas oscilaciones de gran amplitud. Esto es conocido como el fenómeno de Gibbs, que se produce ya que los elementos de la base (funciones trigonométricas) son funciones muy regulares, pero que están aproximando una función discontinua. Para [math] x \leq 1/2 [/math] la serie toma el valor 1 y para [math] x \geq 1/2 [/math] toma el valor 0 (ya que tiene que ser igual a [math] g(x) [/math] en casi todo punto), pero por ser periódica, en [math] x = 1/2 [/math] tomará el valor 1/2, de ahi la dificultad de que la serie truncada se comporte "bien" cerca de este punto.
Para evitar el fenómeno de Gibbs vamos a utilizar las sumas de Cesáro para aproximar [math] g(x) [/math]. Su fórmula es: [math] S_N(x)=\frac{1}{N+1} \sum_{n=0}^N f_n(x). [/math]
x=[0:0.01:1];
y=[ones(50,1);zeros(51,1)];
u=0:10^(-3):1;
N=length(u)-1;
NN=[1,5,10];
w = [ones(N/2,1);zeros(N/2+1,1)];
h = 10^(-3);
f=[];
for j=1:length(NN)
for n=1:NN(j)
f_n=0.5*ones(1,N+1);
for k=1:n
g=(cos(k*pi*u))';
a_k=2*h*w'*g;
f_n=f_n+a_k*cos(k*pi*u);
end
f=[f;f_n];
end
%sumas de Cesaro
subplot(1,3,j)
plot(x,y,'b--',LineWidth=1.5)
hold on
S=zeros(1,N+1);
for n=1:NN(j)
S=S+f(n,:);
end
S=(1./(NN(j)+1))*S;
plot(u,S,'r',LineWidth=1.5)
title('n='+string(NN(j)))
legend('f(x)','$f_n(x)$','interpreter','latex')
ylim([-0.2,1.2])
endAl observar las gráficas, podemos observar que efectivamente el fenómeno de Gibbs desaparece.
4 Cambio de intervalo
Ahora vamos a observar qué hay que hacer cuando una función en vez de estar definida en el intervalo [math] [-\pi,\pi] [/math], está definida en un intervalo [math] [a,b] [/math]. En estos casos, hacemos un cambio de variable, variando así también la base trigonométrica.
Dado el intervalo [math] [a,b] [/math] tomamos [math] T=\frac{b-a}{2} [/math], y el cambio de variable [math] y=\frac{T}{pi}x [/math]. Con este cambio, la base trigonométrica es: [math] \left\{ \frac{1}{2}, \sin\left(\frac{n\pi}{T}y\right), \cos\left(\frac{n\pi}{T}y\right) \right\}_{n\in\mathbb{N}} [/math].
Por este motivo, en los dos ejemplos anteriores, cuyo intervalo era [math] [-1,1] [/math], hemos tomado como base trigonométrica [math] \left\{ \frac{1}{2},\cos( n\pi x),\sin(n\pi x) \right\}_{n\in\mathbb{N}}[/math].
Para entender mejor cómo cambiar de base, vamos a tomar la función [math] h(x)=xe^{-x} [/math] en el intervalo [math] [1,3] [/math] y a aproximarla por su serie de Fourier utilizando la misma base que en los ejemplos anteriores. Esto se debe a que al hacer el cambio de variable para tomar la nueva base lo único que afecta es el tamaño del intervalo, debido a la periodicidad de las funciones.
a=1; b=3;
u_1=a:10^(-3):b;
N=length(u_1)-1;
h=(b-a)/N;
w=ones(N+1,1);
w(1)=1/2; w(N+1)=1/2;
g0=(u_1.*exp(-u_1))';
a_0=h*w'*g0;
nn=[5,10,20];
for j=1:length(nn)
for n=1:nn(j)
f_n=0.5*a_0*ones(1,N+1);
for k=1:n
g1=(u_1.*exp(-u_1).*sin(k*pi*u_1))';
g2=(u_1.*exp(-u_1).*cos(k*pi*u_1))';
a_k=h*w'*g1;
b_k=h*w'*g2;
f_n=f_n+a_k*sin(k*pi*u_1)+b_k*cos(k*pi*u_1);
end
end
subplot(1,3,j)
plot(u_1,f_n,'r','LineWidth',1.5)
hold on
plot(u_1,u_1.*exp(-u_1),'b--','LineWidth',1.5)
title('n='+string(nn(j)))
legend('f(x)','$f_n(x)$','interpreter','latex')
end5 Base trigonométrica compleja
Ahora vamos a tomar la base trigonométrica compleja [math] \{ e^{inx}\}_{n\in\mathbb{Z}}[/math] en [math] L^2[-\pi,\pi][/math], cuyo producto escalar se define como: [math]\langle f,g\rangle_{L^2}= \int_{\Omega} f(x) \cdot \overline {g(x)} \, dx [/math].
Esta base nos va a permitir aproximar funciones comlejas mediante series de Fourier, pero también vamos a poder aproximar funciones reales. Tomamos la función [math] f(x)=4x(1/2−x)^2 [/math] en el intervalo [math] [0,1][/math] y vamos a aproximarla con esta base compleja. Para hacer esto lo primero que tenemos que hacer es expresar la base en el intervalo [math] [0,1][/math] que es: [math] \{ e^{2\pi inx}\}_{n\in\mathbb{Z}}[/math] y una vez tenemos esta base ya podemos aproximar nuestra función mediante la serie de Fourier que es: [math] f(x)=\sum _{k=-\infty}^{\infty} c_k e^{2\pi inx}[/math]. El calculo de los coeficientes [math]c_k[/math] lo vamos a hacer utilizando el producto escalar definido para esta base compleja.
clear all
a=0; b=1;
u_1=a:10^(-3):b;
N=length(u_1)-1;
h=(b-a)/N;
w=ones(N+1,1);
w(1)=1/2; w(N+1)=1/2;
n=10;
nn=[5,10,20];
for j=1:length(nn)
f_n=zeros(1,N+1);
for k=-nn(j):nn(j)
g1=(4*u_1.*(0.5-u_1).^2.*conj(exp(-i*k*2*pi*u_1)))';
g2=(exp(i*k*2*pi*u_1).*conj(exp(i*k*2*pi*u_1)))';
c=h*w'*g1;
d=h*w'*g2;
a_k=c/d;
f_n=f_n+a_k*exp(i*k*pi*2*u_1);
end
subplot(1,3,j)
plot(u_1,f_n,'r',LineWidth=1.5)
hold on
plot(u_1,(4*u_1.*(0.5-u_1).^2),'b--',LineWidth=1.5)
title('n='+string(nn(j)))
legend('$f_n(x)$','f(x)','interpreter','latex')
end