Series de Fourier (Grupo 2 1/2)

Trabajo realizado por estudiantes
Título	Series de Fourier.
Asignatura	EDP
Curso	2023-24
Autores	Alfredo de Lorenzo, Hugo Sanz, Manuel Fdez.
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

1 Introducción

Las Series de Fourier constituyen un conjunto fundamental de herramientas matemáticas empleadas para la aproximación de funciones periódicas mediante la suma infinita de funciones seno y coseno. El término se origina en el trabajo del matemático francés Jean-Baptiste Joseph Fourier, quien desarrolló la teoría durante su investigación sobre la ecuación del calor. Se le atribuye como el pionero en el estudio sistemático de dichas series, y presentó sus resultados iniciales en 1807 y 1811.

El objetivo de las Series de Fourier radica en la capacidad de descomponer cualquier función periódica con un periodo [math]T[/math] en una serie infinita de funciones trigonométricas. En particular, si [math]f(x)[/math] es una función periódica con periodo [math]T[/math], puede expresarse como:

[math]f(x) = \frac{a_0}{2} + \sum_{n=1}^{\infty} \left( a_n \cos \left( \frac{2\pi nx}{T} \right) + b_n \sin \left( \frac{2\pi nx}{T} \right) \right)[/math]

donde [math]a_0[/math], [math]a_n[/math], y [math]b_n[/math] se denominan coeficientes de Fourier de la serie de Fourier de la función [math]f(x)[/math]. Estos coeficientes se calculan mediante integrales definidas sobre el periodo [math]T[/math] de la función.

2 Base trigonométrica

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

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

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

close all
clear
% Definir el rango de x
x_values = linspace(-1, 1, 1000);

% Definir las funciones
num_functions = 4;
seno = cell(1, num_functions);
coseno = cell(1, num_functions);
legends = cell(1, num_functions);
f = @(x) 1/2;

for n = 1:num_functions
        % Funciones seno 
        seno{n} = sin(n *pi* x_values);
        legendss{n} = ['f(x)=sin(', num2str(n), 'pix)'];
    % Funciones coseno
        coseno{n} = cos(n *pi* x_values);
        legendsc{n} = ['f(x)=cos(', num2str(n), 'pix)'];
end

% Graficar todas las funciones en una sola gráfica con colores diferentes
figure(1)
for i = 1:num_functions 
    subplot(1,4,i)
    plot(x_values, seno{i}, 'DisplayName', legendss{i});
    legend(legendss{i});
    xlabel('x');    
    ylabel('f(x)');
end


figure(2)
for i = 1:num_functions 
    subplot(1,4,i)
    plot(x_values, coseno{i}, 'DisplayName', legendsc{i});
    legend(legendsc{i});
    xlabel('x');
    ylabel('f(x)');
end

figure(3)
plot(x_values,ones(size(x_values))*1/2,'DisplayName', '1/2')
legend('f(x)=1/2');
xlabel('x');
ylabel('f(x)');
grid on;
hold off;

3 Aproximación de una función

3.1 Funciones continuas

En primer lugar se estudiarán las funciones continuas, y como se pueden aproximar mediante la serie de Fourier. Para ello se trabajará como ejemplo con la función:

[math]f(x) = x(1-x)[/math]

sobre el intervalo [math][0,1][/math].

Para poder obtener su base trigonométrica en el intervalo [math][-1,1][/math], en primer lugar se comprueba que se pueda extender de forma impar, para ello se obtiene expresión [math]-f(-x)[/math] . Una vez hecho, se de fine la extensión, que en este caso es:

[math] g(x)=\begin{cases} x(1+x), & \text{si } x \in [0,1] \\ x(1-x), & \text{si } x \in [-1,0] \end{cases}[/math]

Una vez se ha obtenido, ya está definida la función en el intervalo [math][-1,1][/math], viendo que es una función continua ya se pueden obtener los coeficientes de Fourier.

[math] \lim_{x\rightarrow 0^{-}} g(x)=0=\lim_{x\rightarrow 0^{+}} g(x) [/math]

Cabe destacar, que al estar trabajando con una función impar, como resultado todos los coeficientes asociados al coseno serán [math]0[/math], al ser pares, y de los asociados al seno uno de cada dos lo será, aquellos que sean pares por ser multiplicación de impares. Una vez calculados los términos, se verán las series de Fourier. Se han representado las series relativas a [math]n=1,5,10[/math] y se muestra como varia el error.

% Definir la función a integrar
fa_par = @(x, n) sin(pi * n * x) .* (x.*(1-x));
fb_par = @(x, n) cos(pi * n * x) .* (x.*(1-x)); 
fa_impar = @(x, n) sin(pi * n * x) .* (x.*(1+x));
fb_impar = @(x, n) cos(pi * n * x) .* (x.*(1+x));
f = @(x) (x.*(1-x));
h = @(x) (x.*(1+x));
% Tolerancia para el método del trapecio
xx=-1:10^(-3):0;
xxx=0:10^(-3):1;
%numero de terminos
i=10;
% Inicializar un vector para almacenar los resultados
b_k = zeros(i, 1);
a_k = zeros(i, 1);

% Calcular las integrales para n desde 1 hasta 20 usando el método del trapecio adaptativo
for n = 1:i
    resultado_seno = trapz(xx,fa_impar(xx, n))+trapz(xxx,fa_par(xxx, n));
    resultado_cos = trapz(xx,fb_impar(xx, n))+trapz(xxx,fb_par(xxx, n)); 

    % Almacenar el resultado
    b_k(n) = resultado_seno;
    a_k(n) = resultado_cos; 
end
%Serie de fourier para n=i
fn = @(x) 0;
for n=1:i
    fn = @(x) b_k(n) * sin(n*pi*x) + fn(x);
end

% Definir la función
x = linspace(-1, 1, 100000); % Generar 100 puntos en el rango [0, 1]
x2 = linspace(0, 1, 100000); % Generar 100 puntos en el rango [0, 1]
y = linspace(-1, 0, 100000); % Generar 100 puntos en el rango [0, 1]

% Trazar la función
figure(1)
hold on 
plot(x, fn(x), 'r--')
plot(x2, f(x2),'blue-')
plot(y, h(y),'yellow-')
title('Comparación f y fn'); % Título del gráfico
legend('fn(x)', 'f(x)', 'h(x)')
grid on; % Activar la cuadrícula
hold off

% Mostrar los resultados
disp('Resultados para coefs del seno usando el método del trapecio:');
disp(b_k);
disp('Resultados para coefs del coseno usando el método del trapecio:');
disp(a_k);

fn = @(x) 0;
for n=1:i
    fn = @(x) b_k(n) * sin(n*pi*x) + fn(x);
    abs_diff = @(x) abs(f(x) - fn(x)).^2;
    integral_result(n) = sqrt(formula_trapecio(abs_diff, 0, 1,tolerancia));
    abs_diffneg = @(x) -abs(f(x) - fn(x));
    [supremo(n), x_supremo(n)] = (fminbnd(abs_diffneg, 0, 1));
end
N=zeros(i);
for n=1:i
    N(n)=n;
end

figure(2)
hold on 
plot(N, log(integral_result),'r-')
plot(N, log(abs(x_supremo)),'b--')
xlabel('n'); % Etiqueta del eje x
ylabel('error'); % Etiqueta del eje y
title('Gráfico fn'); % Título del gráfico
grid on; % Activar la cuadrícula
hold off

Se muestra en los distintos casos, la función que es la roja, y su aproximación mediante serie de Fourier. En la gráfica del logaritmo del error, vemos el lo función roja, es el error uniforme, y la azul es el error en las normas [math] L^{2}[/math] . Se ha empleado el logaritmo, ya que en este caso, no se puede apreciar bien la forma en la que disminuye, y el logaritmo nos da una mejor idea. Se aprecia como ha medida que la n aumenta, el error disminuye, como es lógico.

3.2 Funciones discontinuas

Ahora se pasará a estudiar una función discontinua, y aunque el procedimiento sea similar, lo interesante serán las conclusiones que se puedan sacar. La función que estudiaremos será:

[math] f(x)=\begin{cases} 1, & \text{si } x \in [0,1/2] \\ 0, & \text{si } x \in (1/2,1] \end{cases}[/math]

Dicha función es discontinua, y cada rama es constante, de tal forma que al hacer la extensión que en este caso será par que en este paso será trivial, nos quedará una función definida en el intervalo [math][-1,1][/math] :

[math] f(x)=\begin{cases} 0, & \text{si } x \in [-1,-1/2] \\ 1, & \text{si } x \in (1/2,1/2)\\ 0, & \text{si } x \in [1/2,1]\\ \end{cases}[/math]

Y una vez hecha la extensión par, se calcularán los coeficientes de la serie de Fourier, pero en este caso tiene los términos pares. La idea es ver las diferencias a la hora de aproximar una función continua o discontinua mediante una serie de Fourier.

clc
clear
format long
close all

g = @(x) (x > 0 & x < 0.5);
% Definir la función a integrar
funa = @(x, n) sin(pi * n * x);
funb = @(x, n) cos(pi * n * x);
func = @(x, n) cos(pi * n * x).^2;
fund = @(x, n) sin(pi * n * x).^2;
% Definir los límites de integración
a = -1; %A pesar de ser [-1,1], fuera de [-1/2,1/2] toma el valor 0
b = 1; 
xx= -1/2:10^(-3):1/2;
xxx= -1:10^(-3):1;

% Tolerancia para el método del trapecio adaptativo
tolerancia = 1e-10;
%numero de términos
i=10;

% Inicializar un vector para almacenar los resultados
ak = zeros(i, 1); %coeficientes del seno
bk = zeros(i, 1); %coeficientes del coseno


% Calcular las integrales para n desde 1 hasta i usando el método del trapecio adaptativo
for n = 1:i 
    ak(n) = trapz(xx,funa(xx,n))/trapz(xxx,fund(xxx,n));
    bk(n) = trapz(xx,funb(xx,n))/trapz(xxx,func(xxx,n));
end
a0= 1; %coeficiente de Fourier dado por 1/2
%Serie de fourier hasta n=i
fn = @(x) a0/2;
%Sumas de Cesarò
Sn = @(x) fn(x);
for n=1:i
    fn = @(x) bk(n) * cos(n*pi*x) + fn(x);
    Sn = @(x)fn(x) + Sn(x);
    Snaux=@(x) 1/(n+1) * (Sn(x));
    abs_diffn = @(x) abs(g(x) - Snaux(x)).^2;
    integral_result(n) = sqrt(integral(abs_diffn, 0, 1));
    abs_diffneg = @(x) -abs(g(x) - Snaux(x));
    [supremo(n), x_supremo(n)] = fminbnd(abs_diffneg, 0, 1);
end
Sn = @(x) 1/(i+1) * (Sn(x));

% Definir la función
x = linspace(0, 1, 100000); % Generar 100 puntos en el rango [0, 1]

% Trazar la función
hold on 
plot(x, fn(x))
plot(x,g(x)) %gráfica de la función característica
xlabel('x'); % Etiqueta del eje x
ylabel('f(x)'); % Etiqueta del eje y
title('Gráfico fn'); % Título del gráfico
grid on; % Activar la cuadrícula
hold off

% Mostrar los resultados
disp('Resultados para coefs del seno usando el método del trapecio:');
disp(ak);
disp('Resultados para coefs del coseno usando el método del trapecio:');
disp(bk);

%Vemos que se produce el llamado "fenómeno de Gibbs" por lo que definimos
%otra función usando las sumas de Cesarò

figure

hold on 
plot(x, Sn(x))
plot(x,g(x)) %gráfica de la función característica
xlabel('x'); % Etiqueta del eje x
ylabel('f(x)'); % Etiqueta del eje y
title('Gráfico Sn'); % Título del gráfico
grid on; % Activar la cuadrícula
hold off

N=zeros(i);
for n=1:i
    N(n)=n;
end

figure 

hold on 
plot(N, (integral_result))
plot(N,(abs(x_supremo)))
xlabel('n'); % Etiqueta del eje x
ylabel('error'); % Etiqueta del eje y
title('Gráfico fn'); % Título del gráfico
grid on; % Activar la cuadrícula
hold off

Al ver las gráficas de la función y la serie de Fourier, se observan oscilaciones en la discontinuidad. Dicho suceso es conocido como el fenómeno de Gibbs. Para mitigarlo se emplea la suma de Césaro, la cual es:

[math]S_{N}={\frac{1}{N+1}} \sum_{i=0}^{N} f_{n}[/math]

Y se puede apreciar el efecto que tiene en la función. Viendo los errores en las gráficas, se puede apreciar como el error uniforme es constantemente 0.5, y esto se debe a que la discontinuidad es de salto 1, por lo que es imposible que el supremo de la diferencia sea menor de 0.5. Otra conclusión que se puede apreciar, es como el error desciende de forma mucho más lenta en el caso de la norma [math] L^{2}[/math] al ser una función discontinua, que indica que hay mayor dificultad para adaptarse.

4 Cambio de intervalo

En esta sección se trabajará como modificar un intervalo dado para obtener su base trigonométrica asociada y su correspondiente base de Fourier. Esta es otra forma que se utiliza cuando el intervalo no es ni simétrico respecto de 0, o no se puede conseguir mediante extensiones pares o impares; tal y como se ha tratado en las secciones anteriores.

Sea [math] [a,b] [/math] el intervalo genérico, para poder obtener su base trigonométrica se debe buscar otro intervalo de forma que podamos obtener su base. Estos intervalos son de la forma [math] [-T,T] [/math], simétricos respecto a 0, y se debe verificar la siguiente condición: [math] long(-T,T)= long(a,b) [/math]. Una vez se tenga el nuevo intervalo, la función se debe extender por periodicidad para finalmente mantenernos en el intervalo original.

En dicho caso, se ha visto anteriormente como se calcula la base en dicho intervalo.

Se plantea el siguiente ejemplo para facilitar la comprensión de esta sección.

Dado el intervalo [math][1,3][/math], obtener su base trigonométrica asociada y aproximar en dicho intervalo la función

[math]f(x) = xe^{-x}[/math]

En primer lugar, para obtener la base trigonométrica se buscará el intervalo cuya longitud es 2, entonces se toma el intervalo [math][-1,1][/math].

Por tanto, según lo visto en secciones anteriores, la base trigonométrica en el espacio [math] L^{2}[-1,1][/math] será [math] \Bigl\{ \frac{1}{2}, cos(\pi n x), sen(\pi n x)\Bigl\} [/math]

Para poder aproximar la [math]f[/math] dada por la serie de Fourier en el intervalo modificado, se debe extender la función por periodicidad para volver al intervalo original. En este caso, se debe realizar una traslación de [math]2k[/math] con [math]k \in \mathbb{Z} [/math]. Es decir, los coeficientes de Fourier obtenidos serán de [math]g(x)=f(x+2)=(x+2)e^{-(x+2)}[/math] (con [math]k=1[/math]), para así, a la hora de representar [math]f_{n}[/math] en el intervalo [math][1,3][/math] se corresponda con la función [math]f(x)[/math] original.

A continuación, se muestra el código con el cual se han obtenido los coeficientes de [math]g(x)[/math] asociados al intervalo [math][-1,1][/math]. El ejemplo se muestra para tres casos, la aproximación con 5, 10 y 20 coeficientes.

% Definimos las funciones
f = @(x) (x .* exp(-x));
funa = @(x, n) sin(pi*n*x) .* ((x+2) .* exp(-(x+2)));
funb = @(x, n) cos(pi*n*x) .* ((x+2) .* exp(-(x+2)));
fun0 = @(x) ((x+2).* exp(-(x+2)));

% Tolerancia para el método del trapecio
xx=-1:10^(-3):1;

%numero de terminos
terminos=[5,10,20];
k=1
for i=terminos
    
    % Inicializar un vector para almacenar los resultados
    results = zeros(i, 1);
    
    % Calcular las integrales para n desde 1 hasta 20
    for n = 1:i
        % Utilizar la función integral para calcular la integral definida
        a_k(n) = trapz(xx,funa(xx, n)) ;
        b_k(n) = trapz(xx,funb(xx, n));
    end
    
    %veamos el termino a0
    a0 = trapz(xx,fun0(xx));
    
    % Mostrar los resultados
    disp('Resultados para cada n del seno:');
    disp(a_k);
    disp('Resultados para cada n del coseno:');
    disp(b_k);

En esta parte del código se ha obtenido la correspondiente serie de Fourier,

%Calculemos la serie de Fourier
    fn=@(x) a0/2;
    for n=1:i
        fn = @(x) a_k(n)*sin(n*pi*x)+ b_k(n)*cos(n*pi*x) + fn(x);
    end

Y finalmente obtenemos las gráficas donde se observa como se aproxima la serie de fourier a la función original [math]f(x)[/math] en [math][1,3][/math]

% Trazar la función
    x = linspace(1, 3, 1000);
    subplot(1,3,k)
    plot(x,f(x),'b')
    hold on
    plot(x,fn(x),'r--')
    title('Comparación f y fn para n='+ string(i)); % Título del gráfico
    legend( 'f(x)','fn(x)')
    grid on; % Activar la cuadrícula
    hold off
    k=k+1
end

5 Aproximación de una función con base trigonométrica compleja

Una vez trabajado con bases trigonométricas reales, se verá como se puede trabajar con una base trigonométricas compleja. La base trigonométrica compleja en un intervalo [math] [-T,T] [/math] es de la forma [math]\{ e^{\frac{\pm inx\pi}{T}} \} [/math].

La serie de Fourier compleja es de la forma,

[math]f= \sum_{i=- \infty}^{\infty} c_{k} e^{i \cdot 2nx\pi}[/math]

donde cada coeficiente complejo se obtiene de la siguiente manera,

[math] c_{k} = \frac{1}{2T} \int_{-T}^{T} f(z) e^{\frac{-i\cdot nz\pi}{T}} \, dz [/math]

En el espacio [math]L^{2}[-T,T][/math] de funciones que toman valores complejos, el producto escalar,

[math] (f,g)_{L^[2]} =\int_{-T}^{T} f(x) \overline{g(x)} \, dx [/math]

toma valores complejos. Por tanto, los coeficientes de Fourier serán número complejos, a pesar de que las funciones tomen valores reales.

Esto se particularizará para el intervalo [math][0,1][/math] y la función [math]f(x) = 4x \cdot (\frac{1}{2} - x)^{2}[/math].

Repitiendo el proceso de la sección anterior para definir el intervalo donde se puedan calcular los coeficientes (complejos) de Fourier, se debe tomar [math]T=\frac{1}{2}[/math], entonces el intervalo será [math][-\frac{1}{2},\frac{1}{2}][/math].

Entonces la base de este caso particular será:

[math] \{ e^{\pm i \cdot 2nx\pi} \} [/math]

Nuevamente como se ha modificado el intervalo dado, se debe modificar la función con la cual se obtendrán los coeficientes. En este caso la función será [math]g(x) = f(x+T) = 4(x+\frac{1}{2}) \cdot (\frac{1}{2} - (x+\frac{1}{2}))^{2}[/math]. En el siguiente código se obtendrán estos coeficientes para los primeros 5, 10 y 20 términos de la base.