-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathUseful features
128 lines (103 loc) · 2.68 KB
/
Useful features
1
#### APPROXIMATION ####import mathimport numpy as npimport matplotlib.pyplot as plt%matplotlib inlinedef f(x): return np.sin(x) + 0.5 * xx = np.linspace(-2 * np.pi, 2 * np.pi, 50)plt.grid(True)plt.xlabel('x')plt.ylabel('f(x)')plt.plot(x, f(x), 'b')#### REGRESSION ####reg = np.polyfit(x, f(x), deg=1)ry = np.polyval(reg, x)plt.plot(x, f(x), 'b', label='f(x)')plt.plot(x, ry, 'r.', label='regression')plt.legend(loc=0)plt.grid(True)plt.xlabel('x')plt.ylabel('f(x)')reg = np.polyfit(x, f(x), deg=5)ry = np.polyval(reg, x)plt.plot(x, f(x), 'b', label='f(x)')plt.plot(x, ry, 'r.', label='regression')plt.legend(loc=0)plt.grid(True)plt.xlabel('x')plt.ylabel('f(x)')reg = np.polyfit(x, f(x), 7)ry = np.polyval(reg, x)plt.plot(x, f(x), 'b', label='f(x)')plt.plot(x, ry, 'r.', label='regression')plt.legend(loc=0)plt.grid(True)plt.xlabel('x')plt.ylabel('f(x)')np.allclose(f(x), ry) #checks resultsnp.sum((f(x) - ry) ** 2) / len(x) #MSEmatrix = np.zeros((3 + 1, len(x)))matrix[3, :] = x ** 3matrix[2, :] = x ** 2matrix[1, :] = xmatrix[0, :] = 1reg = np.linalg.lstsq(matrix.T, f(x))[0]regry = np.dot(reg, matrix)plt.plot(x, f(x), 'b', label='f(x)')plt.plot(x, ry, 'r.', label='regression')plt.legend(loc=0)plt.grid(True)plt.xlabel('x')plt.ylabel('f(x)')matrix[3, :] = np.sin(x)reg = np.linalg.lstsq(matrix.T, f(x))[0]ry = np.dot(reg, matrix)plt.plot(x, f(x), 'b', label='f(x)')plt.plot(x, ry, 'r.', label='regression')plt.legend(loc=0)plt.grid(True)plt.xlabel('x')plt.ylabel('f(x)')np.allclose(f(x), ry)np.sum((f(x) - ry) ** 2) / len(x)#### INTERPOLATION ####import scipy.interpolate as spix = np.linspace(-2 * np.pi, 2 * np.pi, 25)def f(x): return np.sin(x) + 0.5 * xipo = spi.splrep(x, f(x), k=1)iy = spi.splev(x, ipo)plt.plot(x, f(x), 'b', label='f(x)')plt.plot(x, iy, 'r.', label='interpolation')plt.legend(loc=0)plt.grid(True)plt.xlabel('x')plt.ylabel('f(x)')#### INTEGRATION ####import scipy.integrate as scidef f(x): return np.sin(x) + 0.5 * xa = 0.5 # left integral limitb = 9.5 # right integral limitx = np.linspace(0, 10)y = f(x)from matplotlib.patches import Polygonfig, ax = plt.subplots(figsize=(7, 5))plt.plot(x, y, 'b', linewidth=2)plt.ylim(ymin=0)# area under the function# between lower and upper limitIx = np.linspace(a, b)Iy = f(Ix)verts = [(a, 0)] + list(zip(Ix, Iy)) + [(b, 0)]poly = Polygon(verts, facecolor='0.7', edgecolor='0.5')ax.add_patch(poly)# labelsplt.text(0.75 * (a + b), 1.5, r'$\int_a^b f(x)dx$', horizontalalignment='center', fontsize=20)plt.figtext(0.9, 0.075, '$x$')plt.figtext(0.075, 0.9, '$f(x)$')ax.set_xticks((a, b))ax.set_xticklabels(('$a$', '$b$'))ax.set_yticks([f(a), f(b)])