Ver el código de «Series de Fourier(PDM)»

{{ 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>