Autor Tema: Interpolación de Chebyshev  (Leído 30025 veces)

0 Usuarios y 1 Visitante están viendo este tema.

Carlos

  • Moderador Global
  • ****
  • Mensajes: 294
Interpolación de Chebyshev
« en: 23/Ago/2019, 13:14:30 pm »
Hoja de cálculo para calcular polinomios de interpolación de hasta 6 puntos.
Se muestran los nodos de chebyshev para conseguir con ellos una interpolación óptima con mínimo error absoluto.


« Última modificación: 22/Dic/2021, 17:00:13 pm por Carlos »

Carlos

  • Moderador Global
  • ****
  • Mensajes: 294
Re:Interpolación de Chebyshev
« Respuesta #1 en: 22/Dic/2021, 20:33:58 pm »
Programa en Python para calcular los coeficientes del polinomio aproximador de una función.


Código: (Python) [Seleccionar]
import matplotlib.pyplot as plt
import numpy as np

num_nodes = 6
interval = [-np.pi/2, np.pi/2]

def f(x):
    return np.sin(x)

def chebyshev_nodes(n, interval, closed_interval=False):
    """Return n chebyshev nodes over interval [a, b]"""
    i = np.arange(n)
    if closed_interval:
       x = np.cos((2*i)*np.pi/(2*(n-1))) # nodes over closed interval [-1, 1]
    else:
       x = np.cos((2*i+1)*np.pi/(2*n)) # nodes over open interval (-1, 1)
    a, b = interval
    return 0.5*(b - a)*x + 0.5*(b + a) # nodes over interval [a, b]

def print_coefs(coefs):
    print("\nCoefficients:")
    for i in range(len(coefs)):
        print(f"    a_{i} = {coefs[i]:.14g}")

def tics(interval, numtics):
    a, b = interval
    return np.linspace(a, b, numtics)

def plot_func(x, y, err_abs, err_rel):
    fig, ax = plt.subplots(3)
    fig.subplots_adjust(left=0.1, bottom=0.05, right=0.95, top=0.95, wspace=None, hspace=0.35)

    ax[0].plot(x, y, linewidth=1)
    ax[0].set_title('Function')
    ax[0].spines['left'].set_position('zero')
    ax[0].spines['right'].set_color('none')
    ax[0].spines['bottom'].set_position('zero')
    ax[0].spines['top'].set_color('none') 

    ax[1].plot(x, err_abs, linewidth=1)
    ax[1].set_title('Polynomial absolute error')
    ax[1].spines['left'].set_position('zero')
    ax[1].spines['right'].set_color('none')
    ax[1].spines['bottom'].set_position('zero')
    ax[1].spines['top'].set_color('none')   

    ax[2].plot(x, err_rel, linewidth=1)
    ax[2].set_title('Polynomial relative error')
    ax[2].spines['left'].set_position('zero')
    ax[2].spines['right'].set_color('none')
    ax[2].spines['bottom'].set_position('zero')
    ax[2].spines['top'].set_color('none')

    plt.show()

def test_errors(interval, poly_coefs, num_dots=1000):
    x_dots = np.linspace(interval[0], interval[1], num_dots)
    y_dots = f(x_dots)
    y_poly_dots = np.polyval(poly_coefs, x_dots)

    # Compute errors
    err_abs = y_dots - y_poly_dots
    err_abs_max = max(np.abs(err_abs))
    err_rel = np.arange(num_dots).astype(float)
    for i in range(len(x_dots)):
        if y_dots[i] == 0:
            if y_poly_dots[i] == 0:
                err_rel[i] = 0.0
            else:
                err_rel[i] = np.inf
        else:
            err_rel[i] = err_abs[i] / y_dots[i]
    err_rel_max = max(np.abs(err_rel))

    # Show errors
    print(f"\nMax absolute error = {err_abs_max:.14g}")
    print(f"Max relative error = {err_rel_max:.14g}")
    print(f"Max polynomial value = {max(y_poly_dots):.14g}")
    print(f"Min polynomial value = {min(y_poly_dots):.14g}")
    plot_func(x_dots, y_dots, err_abs, err_rel)

def main():
    x_dots = chebyshev_nodes(num_nodes, interval, closed_interval=True)
    print(f"Nodes = {x_dots}")
    y_dots = f(x_dots)
    degrees = np.arange(1, num_nodes, 2) # 2 = compute only odd coefficients
    poly_coefs = np.polynomial.polynomial.polyfit(x_dots, y_dots, degrees)
    print_coefs(poly_coefs)
    test_errors([0, interval[1]], np.flip(poly_coefs))

main()
« Última modificación: 24/Dic/2021, 14:21:15 pm por Carlos »