Programa en Python para calcular los coeficientes del polinomio aproximador de una función.
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()
