Diferencia entre revisiones de «Series de Fourier CPP»

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import trapezoid

import numpy as np
import matplotlib.pyplot as plt

# ============================================================
# 0) SERIE DE FOURIER (texto que se imprime al ejecutar)
# ============================================================

def print_fourier_series_info():
    print("\n" + "="*70)
    print("ESTUDIO: Fenómeno de Gibbs en la onda cuadrada 2π-periódica")
    print("="*70)

    print("\nFunción (onda cuadrada 2π-periódica):")
    print("   f(x) =  1   si  0 < x < π")
    print("   f(x) = -1   si -π < x < 0")
    print("   y se extiende 2π-periódicamente: f(x+2π)=f(x).")

    print("\nSerie de Fourier:")
    print("   f(x) ~ (4/π) * Σ_{n=1..∞}  sin((2n-1)x)/(2n-1)")

    print("\nSuma parcial (N términos):")
    print("   S_N(x) = (4/π) * Σ_{n=1..N}  sin((2n-1)x)/(2n-1)")

    print("\nHechos clave (Gibbs):")
    print(" - Cerca de las discontinuidades aparecen oscilaciones.")
    print(" - Al aumentar N, las oscilaciones se 'concentran' cerca del salto,")
    print("   pero el sobreimpulso máximo NO desaparece (tiende a una constante).")
    print(" - Lejos del salto, S_N converge mejor (convergencia puntual en continuidad).")
    print("="*70 + "\n")


# ============================================================
# 1) Definición de la onda cuadrada 2π-periódica
# ============================================================

def square_wave_2pi(x):
    """
    Onda cuadrada 2π-periódica con salto en 0 (y ±π).
    Se define en el intervalo principal (-π, π].
    """
    xr = (x + np.pi) % (2*np.pi) - np.pi  # reduce a (-π, π]
    y = np.where(xr > 0, 1.0, -1.0)
    return y


# ============================================================
# 2) Suma parcial de Fourier
# ============================================================

def fourier_partial_sum_square(x, N):
    """
    S_N(x) = (4/π) * sum_{n=1..N} sin((2n-1)x)/(2n-1)
    Vectorizado para x array.
    """
    k = 2*np.arange(1, N+1) - 1  # 1,3,5,...,(2N-1)
    return (4/np.pi) * np.sum(np.sin(np.outer(k, x)) / k[:, None], axis=0)


# ============================================================
# 3) DIBUJOS: global + zoom Gibbs
# ============================================================

def plot_global_and_zoom():
    # Global
    M = 5000
    x = np.linspace(-np.pi, np.pi, M, endpoint=True)
    f = square_wave_2pi(x)

    plt.figure()
    plt.plot(x, f, linewidth=2, label="f(x) (onda cuadrada)")
    for N in [1, 3, 10, 50, 100]:
        S = fourier_partial_sum_square(x, N)
        plt.plot(x, S, label=f"S_{N}(x)")
    plt.title("Aproximación global de la onda cuadrada con sumas parciales")
    plt.xlabel("x")
    plt.ylabel("amplitud")
    plt.ylim(-1.6, 1.6)
    plt.grid(True, alpha=0.3)
    plt.legend()
    plt.show()

    # Zoom cerca de x=0
    zoom = 0.5
    xz = np.linspace(-zoom, zoom, 3000)
    fz = square_wave_2pi(xz)

    plt.figure()
    plt.plot(xz, fz, linewidth=2, label="f(x)")
    for N in [5, 10, 20, 50, 100]:
        Sz = fourier_partial_sum_square(xz, N)
        plt.plot(xz, Sz, label=f"S_{N}(x)")
    plt.title("Fenómeno de Gibbs (zoom cerca de la discontinuidad x=0)")
    plt.xlabel("x")
    plt.ylabel("amplitud")
    plt.ylim(-1.6, 1.6)
    plt.grid(True, alpha=0.3)
    plt.legend()
    plt.show()


# ============================================================
# 4) ANÁLISIS NUMÉRICO del Gibbs: overshoot/undershoot vs N
# ============================================================

def gibbs_numerical_study(Nmax=300, window=0.3):
    """
    Mide numéricamente:
      - overshoot a la derecha del salto: max(S_N) - 1
      - undershoot a la izquierda: (-1) - min(S_N) con signo ajustado
    Normaliza por el salto (que vale 2).
    """
    jump = 2.0  # salto de -1 a +1

    x_right = np.linspace(0, window, 6000, endpoint=True)
    x_left  = np.linspace(-window, 0, 6000, endpoint=True)

    overshoot_vals = np.zeros(Nmax)
    undershoot_vals = np.zeros(Nmax)

    for N in range(1, Nmax+1):
        Sr = fourier_partial_sum_square(x_right, N)
        Sl = fourier_partial_sum_square(x_left, N)

        overshoot_vals[N-1] = np.max(Sr) - 1.0
        undershoot_vals[N-1] = -1.0 - np.min(Sl)

    overshoot_rel = overshoot_vals / jump
    undershoot_rel = undershoot_vals / jump

    # Gráfica overshoot/undershoot
    N_axis = np.arange(1, Nmax+1)
    plt.figure()
    plt.plot(N_axis, overshoot_rel, label="Sobreoscilación / salto")
    plt.plot(N_axis, undershoot_rel, label="Suboscilación / salto")
    plt.title("Medida numérica del Gibbs (normalizado por el salto)")
    plt.xlabel("N (número de términos)")
    plt.ylabel("amplitud relativa")
    plt.grid(True, alpha=0.3)
    plt.legend()
    plt.show()

    # Resumen en consola
    print("\n" + "-"*70)
    print(f"Resumen numérico local (buscado en ventana ±{window} alrededor del salto x=0)")
    print("El salto es 2 (de -1 a +1). Se muestra el sobreimpulso relativo over/jump.")
    print("-"*70)

    Ns_show = [5, 10, 20, 50, 100, 200, Nmax]
    print("N\tmax(S_N) der.\tmin(S_N) izq.\tOvershoot\tUndershoot\tOver/jump")
    for N in Ns_show:
        Sr = fourier_partial_sum_square(x_right, N)
        Sl = fourier_partial_sum_square(x_left, N)
        max_right = np.max(Sr)
        min_left  = np.min(Sl)
        over = max_right - 1.0
        under = -1.0 - min_left
        print(f"{N}\t{max_right:+.6f}\t\t{min_left:+.6f}\t\t{over:.6f}\t\t{under:.6f}\t\t{over/jump:.6f}")

    # Reporte final (valores límite aproximados numéricamente)
    print("\nObservación:")
    print(" - Al crecer N, el sobreimpulso relativo (over/jump) se estabiliza cerca de un valor constante.")
    print(" - Aunque el ancho de la zona oscilatoria se hace más pequeño, la altura máxima no tiende a 0.")
    print("-"*70 + "\n")


# ============================================================
# 5) CONVERGENCIA LEJOS DE LA DISCONTINUIDAD
# ============================================================

def convergence_away_from_jumps(Nmax=300, exclude=0.15):
    """
    Mide el error supremo ||S_N - f||_∞ evitando:
      - vecindad de x=0
      - vecindad de x=±π
    para mostrar que en zonas continuas el error baja con N.
    """
    M = 6000
    x = np.linspace(-np.pi, np.pi, M, endpoint=True)
    f = square_wave_2pi(x)

    mask = (np.abs(x) > exclude) & (np.abs(np.abs(x) - np.pi) > exclude)

    errors_sup = np.zeros(Nmax)
    for N in range(1, Nmax+1):
        S = fourier_partial_sum_square(x, N)
        errors_sup[N-1] = np.max(np.abs(S[mask] - f[mask]))

    N_axis = np.arange(1, Nmax+1)
    plt.figure()
    plt.plot(N_axis, errors_sup)
    plt.title("Error supremo lejos de discontinuidades (zona continua)")
    plt.xlabel("N (número de términos)")
    plt.ylabel(r"$\|S_N - f\|_\infty$ (fuera de saltos)")
    plt.grid(True, alpha=0.3)
    plt.show()

    # Print breve
    print("Error supremo lejos de saltos:")
    print(f" - N=10:  {errors_sup[9]:.6f}")
    print(f" - N=50:  {errors_sup[49]:.6f}")
    print(f" - N=100: {errors_sup[99]:.6f}")
    print(f" - N={Nmax}: {errors_sup[Nmax-1]:.6f}")
    print()


# ============================================================
# MAIN: que salga TODO al ejecutar
# ============================================================

if __name__ == "__main__":
    print_fourier_series_info()
    plot_global_and_zoom()
    gibbs_numerical_study(Nmax=300, window=0.3)
    convergence_away_from_jumps(Nmax=300, exclude=0.15)




'''

import numpy as np
import matplotlib.pyplot as plt

#
# Parámetros (fijamos f0=440)
#
f0 = 440.0                      # Hz
omega0 = 2*np.pi*f0             # rad/s
T = 1.0/f0                      # periodo

# Mallado temporal, cambiamos periodos
n_periods = 3
t = np.linspace(0, n_periods*T, 6000)

# Coeficientes a_k de la señal (solo senos)
a = {1: 1.0, 2: 0.6, 3: 0.3, 4: 0.1}

def F(t):
    return sum(a[k]*np.sin(k*omega0*t) for k in a)

def S_N(t, N):
    return sum(a[k]*np.sin(k*omega0*t) for k in a if k <= N)

#
# 1) Gráfica: señal original y sumas parciales
# 
plt.figure()
plt.plot(t, F(t), label="F(t)")

for N in [1, 2, 3, 4]:
    plt.plot(t, S_N(t, N), label=f"S_{N}(t)")

plt.xlabel("t (s)")
plt.ylabel("amplitud")
plt.title("F(t) y sumas parciales S_N(t)")
plt.legend()
plt.grid(True)
plt.show()

# 
# 2) Gráfica: 3er armónico aislado
#
plt.figure()
plt.plot(t, 0.3*np.sin(3*omega0*t), label="F_3(t) = 0.3 sin(3 ω0 t)")
plt.plot(t, F(t), label="F(t)", alpha=0.7)
plt.xlabel("t (s)")
plt.ylabel("amplitud")
plt.title("3er armónico y comparación con la señal completa")
plt.legend()
plt.grid(True)
plt.show()

#
# 3) Estudio numérico del error vs N
#   - Error L2 (aprox. por trapecios)
#   - Error sup (norma infinita en la malla)
# 
def L2_error(N):
    e = F(t) - S_N(t, N)
    # Integral sobre el intervalo visualizado [0, n_periods*T]
    return np.sqrt(np.trapz(e**2, t))

def sup_error(N):
    e = F(t) - S_N(t, N)
    return np.max(np.abs(e))

Ns = np.arange(1, 9)
err_L2 = np.array([L2_error(N) for N in Ns])
err_sup = np.array([sup_error(N) for N in Ns])

plt.figure()
plt.plot(Ns, err_L2, marker="o", label="Error L2 (numérico)")
plt.plot(Ns, err_sup, marker="o", label="Error sup (malla)")
plt.xlabel("N (número de términos)")
plt.ylabel("error")
plt.title("Convergencia numérica: error vs N")
plt.legend()
plt.grid(True)
plt.show()
'''

'''
import numpy as np
import matplotlib.pyplot as plt

# Dominio
x = np.linspace(-np.pi, np.pi, 4000)

# Onda cuadrada 2π-periódica: -1 en (-π,0), +1 en (0,π)
f = np.where(x < 0, -1.0, 1.0)
f[np.isclose(x, 0.0)] = 0.0  # en el salto el valor da igual para el dibujo

def fourier_square_partial_sum(x, N):
    """
    Suma parcial S_N de la serie de Fourier de la onda cuadrada:
    f(x) ~ (4/π) * sum_{k=1..∞} sin((2k-1)x)/(2k-1)
    """
    k = np.arange(1, N + 1)
    odd = 2 * k - 1                       # 1,3,5,7,...
    terms = np.sin(np.outer(odd, x)) / odd[:, None]
    return (4/np.pi) * terms.sum(axis=0)

N_list = [3, 10, 50]
plt.figure(figsize=(10, 4.8))
plt.plot(x, f, linewidth=2, label="Onda cuadrada f(x)")

for N in N_list:
    plt.plot(x, fourier_square_partial_sum(x, N), linewidth=1.5, label=f"Suma parcial N={N}")

plt.xlim(-np.pi, np.pi)
plt.ylim(-1.6, 1.6)
plt.xlabel("x")
plt.ylabel("amplitud")
plt.title("Fenómeno de Gibbs: aproximación de una onda cuadrada con series de Fourier")
plt.grid(True, alpha=0.25)
plt.legend()
plt.tight_layout()
plt.show()
'''