Series de Fourier. Yan, Otelo, Miguel
| Trabajo realizado por estudiantes | |
|---|---|
| Título | Series de Fourier. Grupo 6-A |
| Asignatura | EDP |
| Curso | 2023-24 |
| Autores | Miguel Cazorla Pedraza
Otelo Gallego Ayala Yan Wang |
| Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura | |
A lo largo de este artículo trataremos de entender la aproximación de funciones por series trigonométricas, más concretamente, por Series de Fourier. Para ello aproximaremos varias funciones con diferentes propiedades, para así llegar a comprender mejor cómo funcionan estas aproximaciones en casos diferentes.
Contenido
1 Series de Fourier
No tendría sentido intentar entender la aproximación de funciones por series trigonométricas, sin ni siquiera saber que son estás series. Así, comenzamos con una breve explicación teórica de las Series de Fourier.
Tomamos el espacio de Hilbert [math] L^2([-\pi,\pi]) [/math], con producto escalar asociado [math] ( f,g ) _{L^2} = \int_{-\pi}^{\pi} f \cdot g \hspace{5px} dx [/math], norma asociada a dicho producto escalar [math] || f ||_{L^2} = \sqrt{\langle f,g \rangle _{L^2}} [/math] y la base ortogonal [math] \{\frac{1}{2},cos(nx),sen(nx)\}_{n\in \mathbb{N}} [/math] compuesta por funciones [math] 2\pi [/math] periódicas (omitiremos la comprobación de que la base es ortogonal y de que el espacio es de Hilbert). Por ser un sistema ortogonal, podemos aproximar cualquier función [math] f(x) \in L^2([-\pi,\pi]) [/math] mediante la serie
donde los coeficientes [math] a_0, a_n [/math] y [math] b_n [/math] reciben el nombre de Coeficientes de Fourier y vienen dados (por ser la base antes mencionada una base ortognal) por las expresiones
Para finalizar con esta introducción teórica, veamos rápidamente como cambian los elementos de nuestra base trigonométrica dependiendo del intervalo de definición de las funciones de [math] L^2 [/math]. Si el intervalo es la forma [math] [-T,T] [/math], nótese que mediante el cambio de variable [math] y=\frac{xT}{\pi} [/math] en la base anterior obtenemos [math] \{\frac{1}{2},cos(\frac{ny\pi}{T}),sen(\frac{ny\pi}{T})\}_{n\in \mathbb{N}} [/math], siendo esta una base ortogonal de [math] L^2([-T,T]) [/math].
Esta idea se puede generalizar a intervalos de la forma [math] [a,b] [/math] para cualesquiera [math] a [/math] y [math] b [/math] reales. Bastaría tomar la base trigonométrica del espacio [math] L^2([-\frac{b-a}{2},\frac{b-a}{2}]) [/math] , puesto que este intervalo tiene la misma longitud que [math] [a,b] [/math], y los elementos de dicha base son funciones periódicas cuyo periodo es precisamente la longitud [math] b-a [/math].
2 La base trigonométrica
Antes de comenzar con la aproximación de funciones, lo primero que haremos será dibujar los diez primeros términos de nuestra base trigonométrica, [math] \{\frac{1}{2},cos(nx\pi),sen(nx\pi)\}_{n\in \mathbb{N}} [/math] , con el fin de hacernos una idea de la forma geométrica de los elementos de nuestra base. Para dibujarlos utilizamos el siguiente código de matlab:
x = linspace(-1, 1, 1000); % Generar valores de x en el intervalo [-1, 1]
n = 1:10; % Números de términos
% Graficar los términos en una sola imagen con subgráficas
figure;
for i = 1:length(n)
subplot(5, 2, i); % Crear una subgráfica
plot(x, 1/2*ones(size(x)), 'r--', x, cos(n(i)*pi*x), 'g', x, sin(n(i)*pi*x), 'b'); % Graficar los términos
title(['Términos para n = ', num2str(n(i))]); % Título de la subgráfica
xlabel('x'); % Etiqueta del eje x
ylabel(['y_', num2str(n(i))]); % Etiqueta del eje y
legend({'1/2', ['cos(', num2str(n(i)), '\pi x)'], ['sin(', num2str(n(i)), '\pi x)']}); % Leyenda
end
Estas representaciones muestran algo que ya podíamos intuir por la forma de los elementos de nuestra base; a medida que el valor de [math] n [/math] aumenta, el periodo de las funciones [math] cos(nx\pi) [/math] y [math] sen(nx\pi) [/math] disminuye.
3 Aproximación de función continua
Una vez representados algunos de los términos de nuestra base, usaremos estos para aproximar mediante una serie la función [math] f(x)=x(1-x) [/math] en el intervalo [math] [0,1] [/math]. Para facilitarnos la tarea, primero extenderemos de forma impar la función [math] f(x) [/math] al intervalo [math] [-1,1] [/math]. Nótese que los elementos de la forma [math]\frac{1}{2} [/math] y [math] cos(nx\pi) [/math] son pares, y por la definición antes explicada de los coeficientes de Fourier se tiene que en el dominio [math] [-1,1] [/math] dichos coeficientes serán nulos, por lo que la serie estará formada únicamente por términos de la forma [math] b_{n}sen(nx\pi) [/math].
% Definir la función
f = @(x) x.*(1-x);
% Definir el intervalo
x = linspace(0, 1, 100); % Generar 100 puntos entre 0 y 1
x2 = linspace(-1, 0, 100);
% Evaluar la función en los puntos del intervalo
y = f(x);
y1=-f(x);
% Graficar la función
hold on
plot(x2, y1, 'b', 'LineWidth', 2);
plot(x,y,'b','LineWidth', 2)
xlabel('x');
ylabel('Extensión impar def(x) ');
title('Gráfico de f(x) = x(1-x) extendida');
grid on;
Una vez aclarado el por qué de la extensión impar, utilizamos el método del trapecio para calcular numéricamente los coeficientes de Fourier y representamos gráficamente tanto la función [math] f(x) [/math] como las aproximaciones [math] f_{1}(x), f_{5}(x) [/math] y [math] f_{10}(x) [/math]. Además, calculamos el error en las normas [math] L^2 [/math] y uniforme en función del número [math] n [/math] de términos de la serie, y graficamos también los errores obtenidos.
Todo esto queda escrito en el siguiente código de Matlab:
f1 = @(x)(1-x).*x;
coeficientes = zeros(1, 10); %vector de coefs
N=1000; %regla del trapecio
a=0; b=1;
h=(b-a)/(N-1);
u=a:h:b;
w=ones(N,1);
w(1)=1/2; w(N)=1/2;
for k = 1:10 %integrales con la regla y prod escalar para hallar los coefs
f=(f1(u).*sin(k.*pi.*u))';
coeficientes(1,k) = h*w'*f;
end
coeficientes = 2*coeficientes; %coeficientes finales
y = zeros(10,1000); %calculamos las imágenes de las series y las guardamos en la matriz y, que contiene las de los diez primeros términos
y(1,:) = coeficientes(1,1).*sin(pi*u);
for k = 2:10
y(k,:) = y(k-1,:) + coeficientes(1,k).*sin(k*pi*u);
end
figure(1) %rep gráfica de la aproximación
hold on
plot(u,f1(u), "black")
plot(u,y(1,:), "g")
plot(u,y(5,:), "blue")
plot(u,y(10,:), 'red')
hold off
title("Aproximación de f por series de Fourier")
legend("f(x) =(1-x)x", "n =1", "n=5", "n=10")
m = 10; %cálculo de los errores
errores_L2 = zeros(1,m); errores_uniforme = zeros(1,m);
for n=1:m
integrando = (abs(f1(u)-y(n,:))).^2;
g = integrando';
int = h*w'*g;
errores_L2(1,n) = sqrt(int);
errores_uniforme(1,n) = max(abs(f1(u)-y(n,:)));
end
figure(2) %rep gráfica de los errores
hold on
v = 1:1:m;
plot(v,errores_L2)
plot(v,errores_uniforme)
hold off
title("Gráfica de errores")
legend("error en norma L^2", "error en norma uniforme")
En la primera imagen podemos ver como ya a partir de [math] n=5 [/math] la gráfica de [math] f(x) [/math] queda prácticamente solapada por la aproximación [math] f_n(x) [/math]. Esto se confirma en el segundo gráfico, donde vemos que los errores para [math] n=3 [/math] son del orden de [math] 10^{-10} [/math], esto es, que las funciones [math] f_n(x) [/math] aproximan muy correctamente la función [math] f(x) [/math].
4 Aproximación de [math] 1_{[0,\frac{1}{2}]} [/math], una función discontinua
En el apartado apartado anterior pudimos ver como de bien se aproximaba mediante una serie geométrica la función [math] g(x)= x(1-x) [/math]. Ahora probaremos con una función discontinua, con el fin de ver si estas aproximaciones son igual de efectivas para este tipo de funciones. Para ello, trataremos de aproximar la función [math] 1_{[0,\frac{1}{2}]} [/math] definida sobre el intervalo [math] [0,1][/math].
Lo primero que haremos será, de forma parecida a lo hecho en el apartado anterior, tomar la extensión par de la función [math] 1_{[0,\frac{1}{2}]} [/math] en el intervalo [math] [-1,1][/math], para facilitarnos el cálculo de la serie trignométrica. Nótese que en este caso los coefecientes de Fourier que serán nulos serán aquellos correspondientes a las funciones [math] \{sen(nx\pi)\}_{n\in \mathbb{N}} [/math]. Una vez obtenida la serie, nos fijaremos en su imagen en el [math] [0,1][/math] para tener la aproximación buscada.
En este caso, podemos ver como en el gráfico que muestra las aproximaciones ocurre algo que no tenía lugar en el caso de la función continua; al llegar a la discontinuidad, las aproximaciones empiezan a oscilar con periodos cada vez más pequeños, hasta conseguir una pendiente acentuada a lo largo de la discontinuidad. Esto es debido a que la serie trigonométrica está compuesta por funciones regulares, y para aproximar bien puntos discontinuos tienen que hacer estas oscilaciones.
Una forma de amainar este suceso, es mediante las conocidas como sumas de Césaro. El siguiente código, genera tanto las aproximaciones originales como las nuevas aproximaciones con sumas de Césaro, y calcula los errores obtenidos con las sumas de Césaro y con las originales.
close all
m = 100;
coeficientes = zeros(1, m);
N=1000; % Números de puntos
a=0; b=1/2; % Extremos del intervalo
h=(b-a)/(N);
u=a:h:b; % Coordenadas de la partición
w=ones(N+1,1); % Vector de pesos
w(1)=1/2; w(N+1)=1/2;
for k = 1:m
f=(cos(k.*pi.*u))'; % Función base
coeficientes(1,k) = h*w'*f; % Integral hecha por la fórmula del trapecio
end
coeficientes = 2*coeficientes; % Coeficientes del coseno
u = a:1/N:1;
y = zeros(m,N+1);
y(1,:) = (1/2)+ coeficientes(1,1).*cos(pi*u); % f_1
for k = 2:m
y(k,:) = y(k-1,:) + coeficientes(1,k).*cos(k*pi*u); % f_n
end
w = zeros(1,N+1); % este vector lo usaremos para los errores
for p=1:N/2+1
w(1,p) = 1;
end
w1 = ones(1,500); %estos vectores para la rep gráfica
u1 = linspace(0,1/2,500);
w2 = zeros(1,500);
u2 = linspace(1/2,1,500);
subplot(2,2,1) %rep gráfica de la aproximación por series
hold on
plot(u1,w1,'black')
plot(u2,w2,'black')
plot(u,y(1,:), "g")
plot(u,y(5,:), "blue")
plot(u,y(10,:), 'yellow')
plot(u,y(m,:), 'red')
hold off
title("Aproximación de f por series de Fourier")
legend("f", "" , "n =1", "n=5", "n=10", "n="+ num2str(m))
xlabel("x")
ylabel("y")
errores_L2 = zeros(1,m); errores_uniforme = zeros(1,m); %calculamos los erroes
for n=1:m
integrando = (abs(w-y(n,:))).^2;
g = integrando';
int = h*w*g;
errores_L2(1,n) = sqrt(int);
errores_uniforme(1,n) = max(abs(w-y(n,:)));
end
subplot(2,2,2) %rep gráfica de los errores
hold on
v = 1:1:m;
plot(v,errores_L2)
plot(v,errores_uniforme)
hold off
title("Gráfica de errores")
legend("Error en norma L^2", "Error en norma uniforme")
xlabel("N")
ylabel("Error")
M = 100; % sumas de Cesàro
matrizS = zeros(M, N+1);
matrizSF = zeros(M, N+1);
matrizS(1,:) = y(1,:);
matrizSF(1,:) = (1/2).*matrizS(1,:);
for k = 2:M
matrizS(k,:) = matrizS(k-1,:) + y(k,:);
matrizSF(k,:) = (1/(k+1)).*matrizS(k,:);
end
subplot(2,2,3) %rep gráfica
hold on
plot(u1,w1,'black')
plot(u2,w2,'black')
Indices = [1 5 10 M];
for c=1:length(Indices)
plot(u,matrizSF(Indices(c),:))
end
hold off
title(" Sumas de Cesàro para mitigar el fenómeno de Gibbs")
legend("f","" , "n=1", "n=5", "n=10", "n=" + num2str(m))
xlabel("x")
ylabel("y")
erroresC_L2 = zeros(1,m); erroresC_uniforme = zeros(1,m); %calculo de errores de las sumas
for n=1:m
integrandoC = (abs(w-matrizSF(n,:))).^2;
g = integrandoC';
intC = h*w*g;
erroresC_L2(1,n) = sqrt(intC);
erroresC_uniforme(1,n) = max(abs(w-matrizSF(n,:)));
end
subplot(2,2,4) %rep grafica de los errores
hold on
v = 1:1:m;
plot(v,erroresC_L2)
plot(v,erroresC_uniforme)
hold off
title("Gráfica de errores")
legend("Error en norma L^2", "Error en norma uniforme")
xlabel("N")
ylabel("Error")
5 Cambio de intervalo y aproximación de [math] f(x) = xe^{−x} [/math]
En este apartado pretendemos aproximar la función [math] f(x) = xe^{−x} [/math] definida en el intervalo [math] [1,3] [/math]. Aplicando lo explicado en la introducción teórica a inicios de este artículo, la base trigonométrica de [math] L^2([1,3]) [/math] por periodicidad y por la longitud de los intervalos será la misma que la de [math] L^2([-1,1]) [/math], esto es, [math] \{\frac{1}{2},cos(nx\pi),sen(nx\pi)\}_{n\in \mathbb{N}} [/math].
A continuación, de la misma forma que se hizo en apartados anteriores, usamos dicha base para aproximar la función [math] f(x) = xe^{−x} [/math]. Esto queda representado mediante el siguiente código de Matlab:
f =@(x) x.*exp(-x);
m = 20;
coeficientes1 = zeros(1, m); coeficientes2 = zeros(1, m); %creamos los vectores vector de coefs para cada término de la base no cte
N=1000; %aplicamos la regla del trapecio para las integrales
a=1; b=3;
h=(b-a)/(N);
u=a:h:b;
w=ones(N+1,1);
w(1)=1/2; w(N+1)=1/2;
coef0 = h*w'*((f(u))'); %calculamos el a_0
for k = 1:m %calculamos el resto de coefs con el prod escalar y los guardamos
f1=(f(u).*sin(k.*pi.*u))';
f2=(f(u).*cos(k.*pi.*u))'; %function
coeficientes1(1,k) = h*w'*f1;
coeficientes2(1,k) = h*w'*f2;
end
y = zeros(m,N+1); %calculamos las imagenes de las series y las guardamos en la matriz y
y(1,:) = coef0/2+ coeficientes1(1,1).*sin(pi*u) + coeficientes2(1,1).*cos(pi*u);
for k = 2:m
y(k,:) = y(k-1,:) + coeficientes1(1,k).*sin(k*pi*u) + coeficientes2(1,k).*cos(k*pi*u);
end
figure(1) %rep gráfica
hold on
plot(u,f(u),'black')
plot(u,y(5,:), "blue")
plot(u,y(10,:), 'red')
plot(u,y(m,:), 'yellow')
xlabel('x');
ylabel('f(x), f_{5}(x), f_{10}(x), f_{20}(x)');
title('Aproximación de f(x) = xe−x con n= 5, 10 y 20');
hold off
legend("f", "f_{5}", "f_{10}", "f_{20}")
En este gráfico podemos ver que a medida que aumenta la [math] n [/math] la aproximación de la serie trigonométrica se acerca más a la función [math] f(x)[/math] y que las series tienden en los extremos al punto medio de las imágenes de [math] f(x)[/math] en estos. Además, al igual que en el apartado anterior aparece el fenómeno de Gibbs debido a la traslación hecha para crear la serie trigonométrica en el intervalo [math][-1, 1][/math].
6 Base trigonométrica compleja
En este apartado nos piden trabajar con la base trigonométrica compleja [math]\{e^{inx}\}_{n \in \mathbb{Z}}[/math] en [math][-\pi, \pi][/math]. En concreto, buscamos aproximar la función [math]f(x) = 4x(1/2-x)^2[/math] en el intervalo [math][0,1][/math]. Por lo que la base que usaremos será [math]\{e^{2\pi inx}\}_{n \in \mathbb{Z}}[/math]
Debemos tener en cuenta, como indica el enunciado, que en el espacio [math]L^2(-\pi, \pi)[/math] de funciones el producto escalar es:
[math] ( f,g ) _{L^2} = \int_{-\pi}^{\pi} f(x) \overline{g(x)} \hspace{5px} dx [/math]
Por tanto, en el espacio [math]L^2(0,1)[/math] obtendremos los coeficientes de la serie de la siguiente manera:
[math] a_n = \int_{0}^{1} f(x) \overline{e^{2\pi inx}} \hspace{5px} dx = \int_{0}^{1} f(x) e^{-2\pi inx} \hspace{5px} dx [/math]
A continuación se muestra el código empleado para aproximar la función por series con los primeros 5, 10 y 20 términos:
f =@(x) 4*x.*(1/2-x).^2;
m = 20;
coeficientes1 = zeros(1, 2*m+1); %vector de coefs
N=1000; %regla del trapecio para la integral
a=0; b=1;
h=(b-a)/(N);
u=a:h:b;
w=ones(N+1,1);
w(1)=1/2; w(N+1)=1/2;
for k = 1:m
fun1=(f(u).*(exp(-2*1i*pi*k.*u)))'; %coefs para n<0
coeficientes1(1,2*k) = h*w'*fun1;
fun2=(f(u).*(exp(2*1i*pi*k.*u)))'; %coefs para n>0
coeficientes1(1,2*k+1) = h*w'*fun2;
end
fun3=(f(u))'; %coeficiente para n=0
coeficientes1(1,1) = h*w'*fun3;
y = zeros(m+1,N+1); %calculamos las sumas y las guardamos en la matriz y
y(1,:) = coeficientes1(1,1);
for k = 1:m
y(k+1,:) = y(k,:) + coeficientes1(1,2*k).*(exp(-2*1i*k*pi.*u)) + coeficientes1(1,2*k+1).*(exp(2*1i*k*pi.*u));
end
figure(1) %rep gráfica
hold on
plot(u,f(u),'black')
plot(u,y(5,:), "blue")
plot(u,y(10,:), 'red')
plot(u,y(m,:), 'yellow')
hold off
legend("f", "n=5", "n=10", "n ="+ num2str(m))
Observaciones de la gráfica: Al igual que en el apartado anterior, podemos ver que a medida que aumenta el número de términos mejora la aproximación. Además las series en los extremos cumplen la condición de aproximarse al punto medio de la imagen de f en los extremos y aparece el fenómeno de Gibbs.
