Aproximación¶
Importamos los módulos matplotlib.pyplot, numpy y las funciones de numpy para polinomios.
import numpy as np
import matplotlib.pyplot as plt
import numpy.polynomial.polynomial as pol
Precisamos el formato en que queremos imprimir
np.set_printoptions(precision = 2) # solo dos decimales
np.set_printoptions(suppress = True) # no usar notación exponencial
Aproximación polinómica discreta¶
Supongamos que buscamos un polinomio de grado $2$ que aproxime a la función $f(x)=\mathrm{sen}(x)$ de la cual disponemos valores para $5$ puntos
$$x_1=0,x_2=0.5,x_3=1,x_4=1.5,x_5=2$$El polinomio serÃa
$$P_2(x)=a_0+a_1x+a_2x^2$$Y tendremos tres incógnitas $a_0,a_1,a_2$,
Las ecuaciones serÃan
$$ \begin{array}{lll} P_{2}(x_{1})=y_{1} & \quad & a_{0}+a_{1}x_{1}+a_{2}x_{1}^{2}=y_{1}\\ P_{2}(x_{2})=y_{2} & \quad & a_{0}+a_{1}x_{2}+a_{2}x_{2}^{2}=y_{2}\\ P_{2}(x_{3})=y_{3} & \quad & a_{0}+a_{1}x_{3}+a_{2}x_{3}^{2}=y_{3}\\ P_{2}(x_{4})=y_{4} & \quad & a_{0}+a_{1}x_{4}+a_{2}x_{4}^{2}=y_{4}\\ P_{2}(x_{5})=y_{5} & \quad & a_{0}+a_{1}x_{5}+a_{2}x_{5}^{2}=y_{5} \end{array} $$Y el sistema a resolver es:
$$ \left(\begin{array}{ccccc} 1 & x_{1} & x_{1}^{2} \\ 1 & x_{2} & x_{2}^{2} \\ 1 & x_{3} & x_{3}^{2} \\ 1 & x_{4} & x_{4}^{2} \\ 1 & x_{5} & x_{5}^{2} \end{array}\right)\left(\begin{array}{c} a_{0}\\ a_{1}\\ a_{2} \end{array}\right)=\left(\begin{array}{c} y_{1}\\ y_{2}\\ y_{3}\\ y_{4}\\ y_{5} \end{array}\right)\quad \quad (1) $$es decir
$$ Vp=y\quad \quad (1) $$Tenemos más ecuaciones que incógnitas y el sistema no tiene solución. Pero si multiplicamos los dos miembros por la traspuesta de la matriz de coeficientes:
$$ \left(\begin{array}{ccccc} 1 & 1 & 1 & 1 & 1 \\ x_{1} & x_{2} & x_{3}& x_{4} &x_{5}\\ x_{1}^{2} & x_{2}^{2} & x_{3}^{2} & x_{4}^{2} & x_{5}^{2} \end{array}\right) \left(\begin{array}{ccccc} 1 & x_{1} & x_{1}^{2} \\ 1 & x_{2} & x_{2}^{2} \\ 1 & x_{3} & x_{3}^{2} \\ 1 & x_{4} & x_{4}^{2} \\ 1 & x_{5} & x_{5}^{2} \end{array}\right) \left(\begin{array}{c} a_{0}\\ a_{1}\\ a_{2} \end{array}\right)= \left(\begin{array}{ccccc} 1 & 1 & 1 & 1 & 1 \\ x_{1} & x_{2} & x_{3}& x_{4} &x_{5}\\ x_{1}^{2} & x_{2}^{2} & x_{3}^{2} & x_{4}^{2} & x_{5}^{2} \end{array}\right) \left(\begin{array}{c} y_{1}\\ y_{2}\\ y_{3}\\ y_{4}\\ y_{5} \end{array}\right) $$Y ahora tenemos un sistema determinado de $3$ ecuaciones con $3$ incógnitas:
$$ \left(\begin{array}{ccc} n & \sum\limits_{k=1}^{n}x_{k} & \sum\limits_{k=1}^{n}x_{k}^{2}\\ \sum\limits_{k=1}^{n}x_{k} & \sum\limits_{k=1}^{n}x_{k}^{2} & \sum\limits_{k=1}^{n}x_{k}^{3}\\ \sum_\limits{k=1}^{n}x_{k}^{2} & \sum\limits_{k=1}^{n}x_{k}^{3} & \sum\limits_{k=1}^{n}x_{k}^{4} \end{array}\right) \left(\begin{array}{c} a_{0}\\ a_{1}\\ a_{2} \end{array}\right)=\left(\begin{array}{c} \sum\limits_{k=1}^{n}y_{k}\\ \sum\limits_{k=1}^{n}x_{k}y_{k}\\ \sum\limits_{k=1}^{n}x_{k}^{2}y_{k} \end{array}\right) \quad \quad (2) $$es decir tenemos un sistema
$$Cp=d \quad \quad (2)$$que del que tenemos como datos $C$ y $d$ y queremos calcularemos los coeficientes del polinomio de aproximación $p.$
Escribir una función aprox1(f,g,a,b,n) que calcule y dibuje el polinomio de aproximación de grado g para los n puntos equiespaciados de la función f en el intervalo [a,b] (tomando los extremos del intervalo como nodos). La función no devolverá nada: hará los cálculos necesarios para calcular el polinomio de aproximación y dibujará el polinomio de aproximación y los puntos a aproximar. Para ello seguir los pasos siguientes
- Construir los nodos
xcon la funciónnp.linspacey obtenery = f(x) - Construir la matriz de coeficientes $C$ y la matriz de términos independientes $d$ del sistema (2).
- Construir la matriz $V$ del sistema (1) y luego obtener
C = np.dot(V.T,V)yd = np.dot(V.T,y), dondeV.Tes la matriz tranpuesta de $V$ o - Construir directamente C usando el la orden python
np.sumque suma los elementos de un vector y teniendo en cuenta que, por ejemplo,x**2es un vector cuyos elementos son los dexelevados al cuadrado.
- Construir la matriz $V$ del sistema (1) y luego obtener
- Resolver el sistema $Cp=d$ usando la orden
p = np.linalg.solve(C,d). - Obtener 50 puntos en el intervalo $[a,b]$ y los correspondientes valores del polinomio
pusando la ordenpol.polyval. - Dibujar, en el intervalor $[a,b]$:
- Los nodos.
- El polinomio de aproximación obtenido.
Utilizar esta función para obtener y dibujar:
- El polinomio de aproximación de grado $g = 2$ para la función
utilizando $n = 5$ puntos equiespaciados en el intervalo $[a,b]=[0,2].$
- El polinomio de aproximación de grado $g = 4$ para la función
utilizando $n = 10$ puntos equiespaciados en el intervalo $[a,b]=[-2,0].$
%run Ejercicio1.py
Podemos obtener los coeficientes del polinomio de aproximación con la orden polyfit utilizando como argumento de entrada el grado del polinomio de aproximación.
x = np.linspace(0,2,5)
y = np.sin(x)
p = pol.polyfit(x,y,2)
print(p)
Utilizando el dataset Auto MPG de UC Irvine Machine Learning Repository vamos a estimar la potencia (horsepower) necesaria para un coche de un peso dado (weight) y las millas por galón (mpg) que serÃa capaz de recorrer. Los modelos de coche que se manejan en esta base de datos son del año 1970 al 1982. La base de datos también está en kaggle.
Utilizando la librerÃa pandas, leemos el fichero de datos cars.csv y lo almacenamos como un dataframe, que es el equivalente a una hoja de Excel: los datos están organizados en una matriz con una muestra por fila y una variable por columna. Los datos no tienen por qué ser homogéneos. El dataframe puede contener variables de distintos tipos (categóricos, ordinales, numéricos continuos y discretos,...)
import pandas as pd
df = pd.read_csv('http://www.unioviedo.es/compnum/laboratorios_py/new/cars.csv',sep=',')
Otra forma de leer el fichero serÃa descargárselo al directorio de trabajo desde cars y leerlo con
df = pd.read_csv('cars.csv',sep=',') # lee datos separados por comas
Vamos a echar un vistazo a los datos.
df
Para extraer la columna horsepower
x = df['horsepower']
Utilizando el fichero cars que contine datos de Auto MPG de UC Irvine Machine Learning Repository, estimar la potencia (horsepower) necesaria para un coche de un peso dado (weight) y las millas por galón (mpg) que serÃa capaz de recorrer en ciudad. Para ello
- Extraer las variables weight y horsepower.
- Utilizando
pol.polyfitcalcular el polinomio de grado 1 que aproxima estos valores y tomando como variable $x$ weight y como variable $y$ horsepower. - Utilizando
pol.polyvalestimar la potencia (horsepower) de un coche que pesa $3000$ libras. - Dibujar los puntos, el polinomio de aproximación y el punto sobre la recta para un peso de $3000$ libras.
- Extraer la variable mpg y tomando como variable $x$ horsepower y como variable $y$ mpg calcular el polinomio de grado 2 que aproxima estos puntos.
- Con la potencia estimada en el apartado para el coche de $3000$ libras, estimar las millas por galón que podrÃa recorrer este coche en ciudad utilizando el polinomio de aproximación
- Dibujar los datos, el polinomio de aproximación respectivo y el punto sobre el polinomio de grado 3 con el que estimamos las millas por galón en ciudad.
%run Ejercicio2.py
Aproximación polinómica continua¶
Si en lugar de conocer solo algunos puntos de la función conocemos la función completa en un intervalo podemos hacer el ajuste contÃnuo. En este caso, los sumatorios se convierten en integrales:
$$ \left(\begin{array}{ccc} \int_{a}^{b}1 dx & \int_{a}^{b}x dx & \int_{a}^{b}x^2 dx\\ \int_{a}^{b}x dx & \int_{a}^{b}x^2 dx & \int_{a}^{b}x^3 dx\\ \int_{a}^{b}x^2 dx & \int_{a}^{b}x^3 dx & \int_{a}^{b}x^4 dx \end{array}\right)\left(\begin{array}{c} a_{0}\\ a_{1}\\ a_{2} \end{array}\right)=\left(\begin{array}{c} \int_{a}^{b}f(x)dx\\ \int_{a}^{b}x f(x)dx\\ \int_{a}^{b}x^2 f(x)dx \end{array}\right)\quad \quad (3) $$Escribir una función aprox2(f,g,a,b) que calcule y dibuje el polinomio de aproximación de grado g para la función f en el intervalo [a,b]. Para ello seguir los pasos siguientes
- Construir la matriz de coeficientes y la matriz de términos independientes del sistema (3). Utilizar la función
scipy.integrate.quad. Por ejemplo, para calcular el valor de la integral $$I = \int_{a}^{b}f(x)dx$$dondefrom scipy.integrate import quad Ia = quad(f,a,b)[0]
fes una función lambda. Para obtener las funciones integrando, puedes definir una función lambda a partir de otra. Por ejemplog = lambda x: x * f(x)
- Resolver el sistema usando la orden
np.linalg.solvey escribir los coeficientes del polinomio solución de mayor a menor grado. - Obtener los valores del polinomio cuyos coeficientes hemos obtenido usando la orden
pol.polyvalen el intervalo $[a,b]$. - Dibujar, en el intervalor $[a,b]$:
- La función.
- El polinomio de aproximación obtenido.
Utilizar esta función para obtener y dibujar:
- El polinomio de aproximación de grado $g = 2$ para la función
en el intervalo $[a,b]=[0,2].$
- El polinomio de aproximación de grado $g = 4$ para la función
en el intervalo $[a,b]=[-2,0].$
%run Ejercicio3.py
Ejercicios propuestos¶
Aproximación con polinomios ortogonales¶
Si suponemos que hemos estado usando el producto escalar de funciones
$$ \left\langle f(x),g(x)\right\rangle=\int_{-1}^ 1 f(x)g(x)dx$$y la base de polinomios $\{P_0,P_1,P_2\}=\{1,x,x^2\}$ podemos reescribir el sistema como
$$ \left(\begin{array}{ccc} \left\langle P_0,P_0\right\rangle & \left\langle P_0,P_1\right\rangle & \left\langle P_0,P_2\right\rangle\\ \left\langle P_1,P_0\right\rangle & \left\langle P_1,P_1\right\rangle & \left\langle P_1,P_2\right\rangle\\ \left\langle P_2,P_0\right\rangle & \left\langle P_2,P_1\right\rangle & \left\langle P_2,P_2\right\rangle \end{array}\right)\left(\begin{array}{c} a_{0}\\ a_{1}\\ a_{2} \end{array}\right)=\left(\begin{array}{c} \left\langle P_0,f(x)\right\rangle\\ \left\langle P_1,f(x)\right\rangle\\ \left\langle P_2,f(x)\right\rangle \end{array}\right) $$Podemos mejorar el método anterior si conseguimos que la matriz de coeficientes sea diagonal. Lo conseguirÃamos si tuviéramos una base de polinomios ortogonales $\{L_0,L_1,L_2\}$. Entonces se tendrÃa
$$ \left(\begin{array}{ccc} \left\langle L_0,L_0\right\rangle & 0 & 0\\ 0& \left\langle L_1,L_1\right\rangle & 0\\ 0 & 0 & \left\langle L_2,L_2\right\rangle \end{array}\right)\left(\begin{array}{c} a_{0}\\ a_{1}\\ a_{2} \end{array}\right)=\left(\begin{array}{c} \left\langle L_0,f(x)\right\rangle\\ \left\langle L_1,f(x)\right\rangle\\ \left\langle L_2,f(x)\right\rangle \end{array}\right) $$Y ahora tenemos que
$$ a_0=\frac{\left\langle L_0,f(x)\right\rangle}{\left\langle L_0,L_0\right\rangle} \quad a_1=\frac{\left\langle L_1,f(x)\right\rangle}{\left\langle L_1,L_1\right\rangle} \quad a_2=\frac{\left\langle L_2,f(x)\right\rangle}{\left\langle L_2,L_2\right\rangle}\quad \quad (3) $$Estos polinomios existen y son los Polinomios de Legendre. Podemos obtener los polinomios de Legendre del módulo scipy.special
from scipy.special import eval_legendre
xp = np.linspace(-1,1,100)
plt.figure(figsize=(8,6))
plt.plot(xp,0*xp,'k-')
for i in range(6):
plt.plot(xp,eval_legendre(i,xp),label=r'$L_'+str(i)+'$')
plt.legend(bbox_to_anchor=(1, 1), fontsize=16)
plt.title('Polinomios de Legendre',fontsize=16)
plt.show()
Si resolvemos el problema anterior usando los polinomios ortogonales, obtenemos el mismo resultado. Y el polinomio de aproximación será
$$ P_2(x)=a_0\,L_0(x)+a_1\,L_1(x)+a_2\,L_2(x) \quad \quad (4) $$Pero estos polinomiosa de Legendre están definidos para $[-1,1].$
Pero podemos obtenerlos para cualquier intervalo $[a,b]$ realizando una transformación lineal.
Queremos hacer un cambio de variable que nos lleve del intervalo $[-1,1]$ al intervalo $[a,b].$
Si hacemos un cambio de variable lineal
$$x = m\,t +n,$$y si $x = -1$, entonces $t = a$ y
$$-1 = a\,m+n.$$Y si $x = 1$, tenemos que $t = b$ y
$$1 = b\,m+n.$$Si multiplicamos la primera ecuación por $-1$ y sumamos estas dos ecuaciones
$$2 = (b-a)\,m \qquad m = \dfrac{2}{b-a}$$Ahora sustituimos este valor en la primera ecuación
$$-1 = a\,\dfrac{2}{b-a} +n$$Y despejando $n$
$$n =-\dfrac{a+b}{b-a}$$Por lo tanto, el cambio de variable es
$$\fbox{$x=\dfrac{2t-(a+b)}{b-a}$}$$Y si queremos usar los polinomios de Legendre en un intervalo que no sea $[-1,1]$
a = -1.; b = 3.
xp = np.linspace(a,b,100)
plt.figure(figsize=(8,6))
plt.plot(xp,0*xp,'k-')
for i in range(6):
plt.plot(xp,eval_legendre(i,(2*xp-(a+b))/(b-a)),label=r'$L_'+str(i)+'$')
plt.legend(bbox_to_anchor=(1, 1),fontsize=16)
plt.title('Polinomios de Legendre',fontsize=16)
plt.show()
Escribir una función aprox3(f,g,a,b) que calcula el polinomio de aproximación de grado g para la función f en el intervalo [a,b]. Para ello:
- Calcular los coeficientes (3). Utilizar la función
scipy.integrate.quady los polinomios de Legendre descipy.specialtransformados para el intervalo $[a,b]$ - Obtener los valores del polinomio (4) en el intervalo $[a,b]$.
- Dibujar, en el intervalor $[0,2]$:
- La función $f(x)=\cos(x)$
- El polinomio de aproximación de grado $2$ obtenido usando la fórmula (4).
Calcular también la aproximación polinómica de grado 4 de $f_2(x)=\cos(\arctan(x))-\log(x+5)$. Dibujar la función y la aproximación polinómica como en el ejemplo anterior.
%run Ejercicio4.py
Aproximación con funciones trigonométricas: Fourier¶
Si $f$ es una función periódica de periodo $T$ la aproximación mediante funciones trigonométricas resulta especialmente adeduada. UtilizarÃamos el producto escalar
$$ \left\langle f(x),g(x)\right\rangle=\int_{\lambda}^{\lambda+T} f(x)g(x)dx$$donde el intervalo de integración es un intervalo que abarca un periodo completo.
La base de funciones trigonométricas
$$ \left\{ 1,\cos\left(\frac{2\pi x}{T}\right),\mathrm{sen} \left(\frac{2\pi x}{T}\right),\cos\left(2\frac{2\pi x}{T}\right),\mathrm{sen} \left(2\frac{2\pi x}{T}\right),\ldots,\cos\left(n\frac{2\pi x}{T}\right),\mathrm{sen} \left(n\frac{2\pi x}{T}\right)\right\} $$que es una base ortogonal para el producto escalar dado. Entonces, podemos escribir el sistema $Ax=b$ donde
$$ A= \left(\begin{array}{cccccc} \int_{\lambda}^{\lambda+T} 1^2 dx & 0 & 0 & \cdots & 0 & 0\\ 0& \int_{\lambda}^{\lambda+T} \cos^2\left(\frac{2\pi x}{T}\right) dx & 0& \cdots & 0 & 0\\ 0 & 0 & \int_{\lambda}^{\lambda+T} \mathrm{sen}^2 \left(\frac{2\pi x}{T}\right) dx & \cdots & 0 & 0\\ \cdots & \cdots & \cdots & \cdots & \cdots & \cdots\\ 0 & 0 & 0 & \cdots & \int_{\lambda}^{\lambda+T} \cos^2 \left(n\frac{2\pi x}{T}\right) dx & 0\\ 0 & 0 & 0 & \cdots & 0 & \int_{\lambda}^{\lambda+T} \mathrm{sen}^2 \left(n\frac{2\pi x}{T}\right) dx\\ \end{array}\right)\left(\begin{array}{c} \frac{a_{0}}{2}\\ a_{1}\\ b_{1}\\ \cdots\\ a_{n}\\ b_{n} \end{array}\right)= $$$$ =\left(\begin{array}{c} \int_{\lambda}^{\lambda+T} f(x) dx\\ \int_{\lambda}^{\lambda+T} f(x)\cos \left(\frac{2\pi x}{T}\right) dx\\ \int_{\lambda}^{\lambda+T} f(x)\mathrm{sen} \left(\frac{2\pi x}{T}\right) dx\\ \cdots\\ \int_{\lambda}^{\lambda+T} f(x)\cos \left(n\frac{2\pi x}{T}\right) dx\\ \int_{\lambda}^{\lambda+T} f(x)\mathrm{sen} \left(n\frac{2\pi x}{T}\right) dx \end{array}\right) $$Y la función que aproximarÃa a $f$ serÃa de la forma:
$$ F(x)=\frac{a_0}{2}+a_1\cos\left(\frac{2\pi x}{T}\right)+b_1\mathrm{sen} \left(\frac{2\pi x}{T}\right)+a_2\cos\left(2\frac{2\pi x}{T}\right)+b_2\mathrm{sen} \left(2\frac{2\pi x}{T}\right)+\cdots+a_n\cos\left(n\frac{2\pi x}{T}\right)+b_n\mathrm{sen} \left(n\frac{2\pi x}{T}\right) $$Y ahora tenemos que
$$ \frac{a_0}{2}=\frac{\int_{\lambda}^{\lambda+T} f(x) dx}{\int_{\lambda}^{\lambda+T} 1^2 dx} \quad a_1=\frac{\int_{\lambda}^{\lambda+T} f(x)\cos \left(\frac{2\pi x}{T}\right) dx}{\int_{\lambda}^{\lambda+T} \cos^2\left(\frac{2\pi x}{T}\right)dx} \quad b_1=\frac{\int_{\lambda}^{\lambda+T} f(x)\mathrm{sen} \left(\frac{2\pi x}{T}\right) dx}{\int_{\lambda}^{\lambda+T} \mathrm{sen}^2\left(\frac{2\pi x}{T}\right)dx} ,\ldots, \\ a_n=\frac{\int_{\lambda}^{\lambda+T} f(x)\cos \left(n\frac{2\pi x}{T}\right) dx}{\int_{\lambda}^{\lambda+T} \cos^2\left(n\frac{2\pi x}{T}\right)dx} \quad b_n=\frac{\int_{\lambda}^{\lambda+T} f(x)\mathrm{sen} \left(n\frac{2\pi x}{T}\right) dx}{\int_{\lambda}^{\lambda+T} \mathrm{sen}^2\left(n\frac{2\pi x}{T}\right)dx} $$Como se tiene que
$$ \int_{\lambda}^{\lambda+T} 1^2 dx= T, \quad \int_{\lambda}^{\lambda+T} \cos^2\left(k\frac{2\pi x}{T}\right)dx = \frac{T}{2}, \quad \int_{\lambda}^{\lambda+T} \mathrm{sen}^2\left(k\frac{2\pi x}{T}\right)dx = \frac{T}{2} $$para cualquier $k$ entero positivo, los coeficientes de las funciones de la base serán
$$ a_0=\frac{2}{T}\int_{\lambda}^{\lambda+T} f(x) dx, \quad a_1=\frac{2}{T}\int_{\lambda}^{\lambda+T} f(x) \cos\left(\frac{2\pi x}{T}\right) dx, \quad b_1=\frac{2}{T}\int_{\lambda}^{\lambda+T} f(x)\mathrm{sen} \left(\frac{2\pi x}{T}\right) dx, \ldots, a_n=\frac{2}{T}\int_{\lambda}^{\lambda+T} f(x)\cos \left(n\frac{2\pi x}{T}\right) dx, \quad b_n=\frac{2}{T}\int_{\lambda}^{\lambda+T} f(x)\mathrm{sen} \left(n\frac{2\pi x}{T}\right) dx. $$Vamos a dibujar la serie de Fourier desde el término $1$ hasta el término $5$ para la función $f(x)=x$ en el intervalo $[0,3]$ de periodo $T=3$.
%run Ejercicio5.py
