|
|
| Línea 1: |
Línea 1: |
| − | {{ TrabajoED | Series de Fourier(DPM)| [[:Categoría:EDP|EDP]]|[[:Categoría:EDP25/26|2025-26]] | Diego García Raposo
| |
| | | | |
| − | Paula Dopico Muñoz
| |
| − |
| |
| − | Manuel Herreros Zarco}}
| |
| − |
| |
| − | ==Series de Fourier DPM==
| |
| − |
| |
| − | ===Poster===
| |
| − | [[Archivo:FinalPosterDefinitivo.jpeg|center|800px]]]]
| |
| − |
| |
| − | ===Códigos===
| |
| − | <pre>
| |
| − | import numpy as np
| |
| − | import matplotlib.pyplot as plt
| |
| − |
| |
| − |
| |
| − | T = 1.0 ### Longitud del intervalo
| |
| − | n_max = 3 # La frecuencua hasta la que pinta senos y cosenos para representar la base lineal
| |
| − | M = 3000 ### Número de puntos en X
| |
| − |
| |
| − | x = np.linspace(-T, T, M) ### Genera un vector de M puntos entre -T y T
| |
| − |
| |
| − | plt.figure(figsize=(10,5)) ### Tamaño de la gráfica, pa q sea grande o chiquita pero no afecta a los puntos
| |
| − |
| |
| − |
| |
| − | phi0 = (1/np.sqrt(2*T)) * np.ones_like(x) ### Genera un vector de 3000 ptos con cada uno la función cte de la basepara luego representarlo
| |
| − | plt.plot(x, phi0, linewidth=3) ### Pinta la linea generada por la función phi0
| |
| − |
| |
| − |
| |
| − | for n in range(1, n_max + 1): ### Pinta las funciones seno y coseno asociadas a cada frecuencia en el plot
| |
| − | phi_c = (1/np.sqrt(T)) * np.cos(n*np.pi*x/T) ### funcion coseno
| |
| − | phi_s = (1/np.sqrt(T)) * np.sin(n*np.pi*x/T) ### funcion seno
| |
| − |
| |
| − | plt.plot(x, phi_c)
| |
| − | plt.plot(x, phi_s)
| |
| − |
| |
| − | plt.title("Base trigonométrica ortonormal en [-T, T]") ### SIstema de representación de python
| |
| − | plt.xlabel("x")
| |
| − | plt.grid(True)
| |
| − | plt.tight_layout()
| |
| − | plt.show()
| |
| − |
| |
| − | <pre>
| |
| − | #nuevo código
| |
| − | import math
| |
| − | import numpy as np
| |
| − | import matplotlib.pyplot as plt
| |
| − |
| |
| − |
| |
| − | # --------- función de ejemplo ----------
| |
| − | def f(x):
| |
| − | return x # cambia aquí si quieres
| |
| − |
| |
| − |
| |
| − | # --------- coeficientes numéricos ----------
| |
| − | def fourier_coeffs_numeric(L, N, M=20000):
| |
| − | xs = np.linspace(0.0, L, M)
| |
| − | fx = f(xs)
| |
| − |
| |
| − | a0 = (2.0 / L) * np.trapz(fx, xs)
| |
| − |
| |
| − | a = np.zeros(N + 1)
| |
| − | b = np.zeros(N + 1)
| |
| − |
| |
| − | for n in range(1, N + 1):
| |
| − | c = np.cos(2.0 * math.pi * n * xs / L)
| |
| − | s = np.sin(2.0 * math.pi * n * xs / L)
| |
| − |
| |
| − | a[n] = (2.0 / L) * np.trapz(fx * c, xs)
| |
| − | b[n] = (2.0 / L) * np.trapz(fx * s, xs)
| |
| − |
| |
| − | return a0, a, b
| |
| − |
| |
| − |
| |
| − | def S_N(x, L, N, a0, a, b):
| |
| − | val = a0 / 2.0
| |
| − | for n in range(1, N + 1):
| |
| − | val += (
| |
| − | a[n] * np.cos(2.0 * math.pi * n * x / L)
| |
| − | + b[n] * np.sin(2.0 * math.pi * n * x / L)
| |
| − | )
| |
| − | return val
| |
| − |
| |
| − |
| |
| − | # --------- visualizar error en muchos puntos ----------
| |
| − | def visualize_error(Ls, N, num_points=500):
| |
| − |
| |
| − | xs = np.linspace(0.001, min(Ls), num_points) # intervalo común
| |
| − |
| |
| − | plt.figure()
| |
| − |
| |
| − | for L in Ls:
| |
| − |
| |
| − | a0, a, b = fourier_coeffs_numeric(L, N)
| |
| − |
| |
| − | approx = S_N(xs, L, N, a0, a, b)
| |
| − | exact = f(xs)
| |
| − |
| |
| − | error = np.abs(approx - exact)
| |
| − |
| |
| − | plt.plot(xs, error, label=f"L={L}")
| |
| − |
| |
| − | plt.yscale("log")
| |
| − | plt.xlabel("x")
| |
| − | plt.ylabel("|S_N(x) - f(x)|")
| |
| − | plt.title(f"Error puntual para N={N}")
| |
| − | plt.legend()
| |
| − | plt.grid(True, which="both")
| |
| − | plt.show()
| |
| − |
| |
| − |
| |
| − | if __name__ == "__main__":
| |
| − |
| |
| − | Ls = [1.0, 4.0]
| |
| − | N = 20
| |
| − |
| |
| − | visualize_error(Ls, N, num_points=800)
| |
| − |
| |
| − |
| |
| − |
| |
| − | #nuevo código
| |
| − | import numpy as np
| |
| − | import matplotlib.pyplot as plt
| |
| − |
| |
| − | # =====================================================
| |
| − | # 1) INTERVALO Y FUNCIÓN
| |
| − | # =====================================================
| |
| − | a, b = 0.0, 1.0 ### Inicio y final del intervalo
| |
| − | L = b - a ### Longitud intervalo
| |
| − |
| |
| − | def f(x):
| |
| − | return np.where((x >= 0.0) & (x < 0.5), 1.0, 0.0) ### Funcion indicatriz
| |
| − |
| |
| − | # =====================================================
| |
| − | # 2) MUESTREO
| |
| − | # =====================================================
| |
| − | M = 8192 ### Numero de puntos en x
| |
| − | x = np.linspace(a, b, M, endpoint=False) ### Genera un vector de M puntos entre a y b, pero no incluye el punto b
| |
| − | fx = f(x) ### Aplica la función a cada elementos del vector del linspace x tanto es un vector de lomgitud M
| |
| − |
| |
| − | # =====================================================
| |
| − | # 3) COEFICIENTES REALES DIRECTOS (SIN FFT)
| |
| − | # =====================================================
| |
| − | def coeffs_real_direct(fx, N, a=0.0, L=1.0):
| |
| − | M = fx.size
| |
| − | j = np.arange(M)
| |
| − | xj = a + L*j/M # mismos puntos que linspace(endpoint=False)
| |
| − |
| |
| − | a0 = fx.mean() # (1/M) sum f_j
| |
| − |
| |
| − | k = np.arange(1, N+1) ### np array, es decir una lista de las frecuencias ordenadas desde 0 hasta N, siendo N la mayor de la suma parcial deseada
| |
| − | theta = 2*np.pi * k[:, None] * (xj[None, :] - a) / L ### Crea una matriz donde el elemento Aij es 2*pi*k_i*x_j
| |
| − |
| |
| − | ak = (2/M) * (fx[None, :] * np.cos(theta)).sum(axis=1) ### Calcula un vector donde cada elemento es el coeficiente asociado a la frecuencua ki
| |
| − | bk = (2/M) * (fx[None, :] * np.sin(theta)).sum(axis=1) ### Lo mismo pero el seno
| |
| − |
| |
| − | return a0, ak, bk
| |
| − |
| |
| − | # =====================================================
| |
| − | # 4) RECONSTRUCCIÓN S_N EN BASE REAL
| |
| − | # =====================================================
| |
| − | def partial_sum_real_direct(N, x, fx, a=0.0, L=1.0):
| |
| − | a0, ak, bk = coeffs_real_direct(fx, N, a=a, L=L)
| |
| − |
| |
| − | k = np.arange(1, N+1) ### Vector de frecuencias
| |
| − | theta = 2*np.pi * k[:, None] * (x[None, :] - a) / L ### Genera una matriz elemento Aij es igual que antes pero con el cambio de intervalo para poder aplicar fourier
| |
| − |
| |
| − | sN = a0 \
| |
| − | + (ak[:, None] * np.cos(theta)).sum(axis=0) \ ### multiplica cada coeficiente con su coseno respectivo, en forma de vector por cada punto del linspace
| |
| − | + (bk[:, None] * np.sin(theta)).sum(axis=0) ### lo mismo para el seno
| |
| − |
| |
| − | return sN ### devuelve la serie de fourier
| |
| − |
| |
| − | # =====================================================
| |
| − | # 5) VISUALIZACIÓN
| |
| − | # =====================================================
| |
| − | Ns = [1, 3, 5, 10, 40, 80] ### Representación python
| |
| − |
| |
| − | cols = 3
| |
| − | rows = (len(Ns) + cols - 1) // cols
| |
| − |
| |
| − | fig, axes = plt.subplots(rows, cols, figsize=(12, 3*rows))
| |
| − | axes = axes.flatten()
| |
| − |
| |
| − | for i, N in enumerate(Ns):
| |
| − | sN = partial_sum_real_direct(N, x, fx, a=a, L=L)
| |
| − |
| |
| − | axes[i].plot(x, fx, color='black', linewidth=2, label="f(x)")
| |
| − | axes[i].plot(x, sN, color='red', linewidth=1.6, label="S_N(x) (sin FFT)")
| |
| − | axes[i].set_title(f"N = {N}")
| |
| − | axes[i].grid(True)
| |
| − |
| |
| − | if i == 0:
| |
| − | axes[i].legend(loc="best")
| |
| − |
| |
| − | # Quitar paneles vacíos
| |
| − | for j in range(i+1, len(axes)):
| |
| − | fig.delaxes(axes[j])
| |
| − |
| |
| − | plt.tight_layout()
| |
| − | plt.show()
| |
| − |
| |
| − | #nuevo código
| |
| − |
| |
| − | import numpy as np
| |
| − | import matplotlib.pyplot as plt
| |
| − |
| |
| − | # =====================================================
| |
| − | # CONFIG
| |
| − | # =====================================================
| |
| − | a, b = 0.0, 1.0
| |
| − | L = b - a
| |
| − |
| |
| − | M_fft = 8192
| |
| − | h_L2 = 1e-5
| |
| − | Ns = [1, 2, 4, 8, 16, 32, 64, 128, 256, 512]
| |
| − |
| |
| − | m_max = 3 # vamos a construir f0, f1, f2, f3 a mano
| |
| − |
| |
| − | # =====================================================
| |
| − | # 1) f0 y sus primitivas exactas f1,f2,f3
| |
| − | # =====================================================
| |
| − | def f0(x):
| |
| − | return np.where((x >= 0.0) & (x < 0.5), 1.0, 0.0)
| |
| − |
| |
| − | def f1(x):
| |
| − | # ∫0^x f0(t) dt
| |
| − | return np.where(x < 0.5, x, 0.5)
| |
| − |
| |
| − | def f2(x):
| |
| − | # ∫0^x f1(t) dt
| |
| − | # x<0.5: x^2/2
| |
| − | # x>=0.5: (1/2)x - 1/8
| |
| − | return np.where(x < 0.5, 0.5*x**2, 0.5*x - 0.125)
| |
| − |
| |
| − | def f3(x):
| |
| − | # ∫0^x f2(t) dt
| |
| − | # x<0.5: x^3/6
| |
| − | # x>=0.5: (1/4)x^2 - (1/8)x + 1/48
| |
| − | return np.where(x < 0.5, (x**3)/6.0, 0.25*x**2 - 0.125*x + (1.0/48.0))
| |
| − |
| |
| − | f_funcs = [f0, f1, f2, f3] # lista de funciones exactas
| |
| − |
| |
| − | # =====================================================
| |
| − | # 2) MALLAS Y VALORES "TRUE" EN L2
| |
| − | # =====================================================
| |
| − | x_L2 = np.arange(a, b, h_L2)
| |
| − | f_list = [f_funcs[m](x_L2) for m in range(m_max + 1)] # f0..f3 evaluadas en x_L2
| |
| − |
| |
| − | # =====================================================
| |
| − | # 3) Fourier por FFT sobre una malla uniforme para coeficientes
| |
| − | # (ojo: extensión periódica implícita de lo que haya en [0,1))
| |
| − | # =====================================================
| |
| − | x_fft = np.linspace(a, b, M_fft, endpoint=False)
| |
| − |
| |
| − | def fourier_coeffs_up_to_from_samples(fx_samples):
| |
| − | """
| |
| − | Devuelve coeffs_up_to(N) -> (a0, ak, bk) para:
| |
| − | f(x) ~ a0 + sum_{k=1}^N [ak cos(2πk(x-a)/L) + bk sin(2πk(x-a)/L)].
| |
| − | """
| |
| − | C = np.fft.fft(fx_samples) / M_fft
| |
| − | a0 = C[0].real
| |
| − |
| |
| − | def coeffs_up_to(N):
| |
| − | ck = C[1:N+1]
| |
| − | ak = 2.0 * ck.real
| |
| − | bk = -2.0 * ck.imag
| |
| − | return a0, ak, bk
| |
| − |
| |
| − | return coeffs_up_to
| |
| − |
| |
| − | def partial_sum_real(coeffs_up_to, N, xgrid):
| |
| − | a0, ak, bk = coeffs_up_to(N)
| |
| − | s = a0 * np.ones_like(xgrid, dtype=float)
| |
| − | for k in range(1, N+1):
| |
| − | ang = 2*np.pi*k*(xgrid - a)/L
| |
| − | s += ak[k-1]*np.cos(ang) + bk[k-1]*np.sin(ang)
| |
| − | return s
| |
| − |
| |
| − | def L2_error(f_true_vals, s_vals, xgrid):
| |
| − | e = f_true_vals - s_vals
| |
| − | return np.sqrt(np.trapz(e**2, xgrid))
| |
| − |
| |
| − | # Precomputamos coeffs_up_to para cada m (muestreando en x_fft)
| |
| − | coeffs_list = []
| |
| − | for m in range(m_max + 1):
| |
| − | f_m_on_fft = f_funcs[m](x_fft) # como es explícita, no hace falta interp
| |
| − | coeffs_list.append(fourier_coeffs_up_to_from_samples(f_m_on_fft))
| |
| − |
| |
| − | # =====================================================
| |
| − | # 4) EXPERIMENTO: error L2 vs N para cada m
| |
| − | # =====================================================
| |
| − | plt.figure(figsize=(7,5))
| |
| − |
| |
| − | for m in range(m_max + 1):
| |
| − | f_true = f_list[m]
| |
| − | coeffs_up_to = coeffs_list[m]
| |
| − | E = []
| |
| − | for N in Ns:
| |
| − | sN = partial_sum_real(coeffs_up_to, N, x_L2)
| |
| − | E.append(L2_error(f_true, sN, x_L2))
| |
| − | plt.loglog(Ns, E, marker='o', label=f"m={m} (f_{m} exacta)")
| |
| − |
| |
| − | plt.grid(True, which="both")
| |
| − | plt.xlabel("N")
| |
| − | plt.ylabel(r"$\|f_m - S_N\|_{L^2(0,1)}$")
| |
| − | plt.title("Convergencia L2 vs N usando primitivas EXACTAS (sin integrar numéricamente)")
| |
| − | plt.legend()
| |
| − | plt.tight_layout()
| |
| − | plt.show()
| |
| − |
| |
| − | # =====================================================
| |
| − | # 5) Visualización de f_m exactas
| |
| − | # =====================================================
| |
| − | fig, axes = plt.subplots(m_max + 1, 1, figsize=(9, 2.4*(m_max+1)))
| |
| − | if m_max == 0:
| |
| − | axes = [axes]
| |
| − |
| |
| − | for m in range(m_max + 1):
| |
| − | axes[m].plot(x_L2, f_list[m], linewidth=1.8)
| |
| − | axes[m].grid(True)
| |
| − | axes[m].set_title(f"f_{m}(x) exacta (integración aplicada {m} veces)")
| |
| − |
| |
| − | plt.tight_layout()
| |
| − | plt.show()
| |
| − | </pre>
| |
