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

2 Base trigonométrica

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

% 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

3.2 Funciones discontinuas

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

4 Cambio de intervalo

En esta sección se trabajará como modificar un intervalo dado para obtener su base trigonométrica asociada.

Sea [math] [a,b] [/math] el intervalo dado, 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].

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 además aproximar la función [math]f(x) = xe^{-x}[/math] en dicho intervalo.

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 será [math]\{ \frac{1}{2}, cos(\pi n x), sen(\pi n x)\}[/math]

Para poder aproximar la [math]f[/math] dada por la serie de Fourier, se debe trasladarla una longitud de 2 para poder trabajar con la base obtenida. Es decir, los coeficientes de Fourier obtenidos serán de [math]g(x)=f(x+2)=(x+2)e^{-(x+2)}[/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 obtener los 20 primeros coeficientes de Fourier.

% 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
i=20; 

% Inicializar un vector para almacenar los resultados
results = zeros(i, 1);

% Calcular los coeficientes de Fourier
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));
fn=@(x) a0/2;

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

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

% Mostrar los resultados
disp('Resultados para cada n del seno:');
disp(a_k);
disp('Resultados para cada n del coseno:');
disp(b_k);

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);
plot(x,f(x),'b')
hold on
plot(x,fn(x),'r--')
title('Comparación f y fn'); % Título del gráfico
legend( 'f(x)','fn(x)')
grid on; % Activar la cuadrícula
hold off

Archivo:Coches