Aproximación

Contenidos

Para que inserte las figuras en la notebook, en lugar de abrir cada figura en una ventana nueva

In [1]:
%matplotlib inline

Con from __future__ import division hacemos que el resultado de la división con / dé siempre como resultado un float. Si no lo incluyéramos la división con / entre dos enteros tendría como resultado otro entero.

Importamos los módulos matplotlib.pyplot y numpy.

In [2]:
from __future__ import division
import numpy as np
import matplotlib.pyplot as plt

Precisamos el formato en que queremos imprimir

In [3]:
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)=\cos(x)$ de la cual solo disponemos de los valores en los $5$ puntos

$$x_0=-1,x_1=-0.5,x_2=0,x_3=0.5,x_4=1$$

El polinomio será

$$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_{0})=y_{0} & \quad & a_{0}+a_{1}x_{0}+a_{2}x_{0}^{2}=y_{0}\\ 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_{3}+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} \end{array} $$

Y el sistema a resolver es:

$$ \left(\begin{array}{ccccc} 1 & x_{0} & x_{0}^{2} \\ 1 & x_{1} & x_{1}^{2} \\ 1 & x_{2} & x_{2}^{2} \\ 1 & x_{3} & x_{3}^{2} \\ 1 & x_{4} & x_{4}^{2} \end{array}\right)\left(\begin{array}{c} a_{0}\\ a_{1}\\ a_{2} \end{array}\right)=\left(\begin{array}{c} y_{0}\\ y_{1}\\ y_{2}\\ y_{3}\\ y_{4} \end{array}\right) $$

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_{0} & x_{1} & x_{2} & x_{3}& x_{4}\\ x_{0}^{2} & x_{1}^{2} & x_{2}^{2} & x_{3}^{2} & x_{4}^{2} \end{array}\right) \left(\begin{array}{ccccc} 1 & x_{0} & x_{0}^{2} \\ 1 & x_{1} & x_{1}^{2} \\ 1 & x_{2} & x_{2}^{2} \\ 1 & x_{3} & x_{3}^{2} \\ 1 & x_{4} & x_{4}^{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_{0} & x_{1} & x_{2} & x_{3}& x_{4}\\ x_{0}^{2} & x_{1}^{2} & x_{2}^{2} & x_{3}^{2} & x_{4}^{2} \end{array}\right) \left(\begin{array}{c} y_{0}\\ y_{1}\\ y_{2}\\ y_{3}\\ y_{4} \end{array}\right) $$

Y ahora tenemos un sistema determinado de $3$ ecuaciones con $3$ incógnitas:

$$ \left(\begin{array}{ccc} n+1 & \sum_{k=0}^{n}x_{k} & \sum_{k=0}^{n}x_{k}^{2}\\ \sum_{k=0}^{n}x_{k} & \sum_{k=0}^{n}x_{k}^{2} & \sum_{k=0}^{n}x_{k}^{3}\\ \sum_{k=0}^{n}x_{k}^{2} & \sum_{k=0}^{n}x_{k}^{3} & \sum_{k=0}^{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_{k=0}^{n}y_{k}\\ \sum_{k=0}^{n}x_{k}y_{k}\\ \sum_{k=0}^{n}x_{k}^{2}y_{k} \end{array}\right) \quad \quad (1) $$

Ejercicio 1

Escribir una programa que calcule el polinomio de aproximación de grado dos para los nodos $x_0=-1,x_1=-0.5,x_2=0,x_3=0.5,x_4=1$ de la función $f(x)=\cos(x)$. Para ello:

  • Construir la matriz de coeficientes y la matriz de términos independientes del sistema (1) e imprimirlos.
  • Resolver el sistema usando la orden numpy.linalg.solve e imprimir la solución.
  • Obtener los valores del polinomio cuyos coeficientes hemos obtenido usando la orden numpy.polyval en el intervalo $[-1,1]$.
  • Dibujar, en el intervalor $[-1,1]$:
    • Los nodos.
    • El polinomio de aproximación obtenido.
In [4]:
%run Ejercicio1.py

 Matriz del sistema

[[ 5.    0.    2.5 ]
 [ 0.    2.5   0.  ]
 [ 2.5   0.    2.12]]


 Término independiente

[ 3.84  0.    1.52]


 Solución del sistema

[ 0.99  0.   -0.46]

Ejercicio 2

Escribir una programa que calcule el polinomio de aproximación de grado cuatro para 10 nodos equiespaciados en el intervalo $[-1,1]$ de la función $f(x)=\cos(\arctan(x))-e^{x^2}\log(x+2)$. Dibujar, en el intervalor $[-1,1]$:

  • Los nodos.
  • El polinomio de aproximación obtenido.
In [5]:
%run Ejercicio2.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.

In [6]:
a2 = np.polyfit(x,y,2)
print(a2)

[-1.09 -1.16  0.31]

Aproximación polinómica continua

Pero 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_{-1}^{1}1 dx & \int_{-1}^{1}x dx & \int_{-1}^{1}x^2 dx\\ \int_{-1}^{1}x dx & \int_{-1}^{1}x^2 dx & \int_{-1}^{1}x^3 dx\\ \int_{-1}^{1}x^2 dx & \int_{-1}^{1}x^3 dx & \int_{-1}^{1}x^4 dx \end{array}\right)\left(\begin{array}{c} a_{0}\\ a_{1}\\ a_{2} \end{array}\right)=\left(\begin{array}{c} \int_{-1}^{1}f(x)dx\\ \int_{-1}^{1}x f(x)dx\\ \int_{-1}^{1}x^2 f(x)dx \end{array}\right)\quad \quad (2) $$

Ejercicio 3

Escribir una programa que calcule el polinomio de aproximación de grado dos para la función $f(x)=\cos(x)$ en el intervalo $[-1,1]$. Para ello:

  • Construir la matriz de coeficientes y la matriz de términos independientes del sistema (2) e imprimirlos. Utilizar la función scipy.integrate.quad.
  • Resolver el sistema usando la orden numpy.linalg.solve y escribir los coeficientes del polinomio solución de mayor a menor grado.
  • Obtener los valores del polinomio cuyos coeficientes hemos obtenido usando la orden numpy.polyval en el intervalo $[-1,1]$.
  • Dibujar, en el intervalor $[-1,1]$:
    • La función $f(x)=\cos(x)$
    • El polinomio de aproximación obtenido.
  • Calcular el error relativo $$E_r = \frac{\left\Vert \cos(x)-P(x)\right\Vert}{\left\Vert \cos(x)\right\Vert}$$ usando la función numpy.linalg.norm.
In [7]:
%run Ejercicio3.py

 Matriz del sistema

[[ 2.    0.    0.67]
 [ 0.    0.67  0.  ]
 [ 0.67  0.    0.4 ]]


 Término independiente

[ 1.68  0.    0.48]


 Coeficientes del polinomio (de mayor a menor grado):

[-0.47  0.    1.  ]

Er =  0.00388705029515

Ejercicio 4

Dada la función $f(x)=\cos(\arctan(x))-e^{x^2}\log(x+2)$, calcular la aproximación polinómica de $f$ de grado 4. Dar los coeficientes, el error relativo y dibujar la función y la aproximación polinómica, como en el ejemplo anterior, en el intervalo $[-1,1]$.

In [8]:
%run Ejercicio4.py

Coeficientes del polinomio (de mayor a menor grado):
[-0.08 -1.03 -1.   -0.38  0.3 ]
Error relativo:
0.0346826444554

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

In [9]:
from scipy.special import legendre

def L(x,n):
    Leg = legendre(n)
    y = Leg(x)
    return y

xx = np.linspace(-1,1,100)
plt.plot(xx,0*xx,'k-')
for i in range(6):
    plt.plot(xx,L(xx,i),label=r'$L_'+str(i)+'$') 
plt.legend(bbox_to_anchor=(1, 1), loc='upper left', ncol=1)  
plt.title('Polinomios de Legendre')
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) $$

Ejercicio 5

Escribir una programa que calcule el polinomio de aproximación de grado dos para la función $f(x)=\cos(x)$ en el intervalo $[-1,1]$. Para ello:

  • Calcular los coeficientes (3). Utilizar la función integrate.quad de scipy y los polinomios de Legendre de scipy.special.
  • Obtener los valores del polinomio (4) en el intervalo $[-1,1]$.
  • Dibujar, en el intervalor $[-1,1]$:
    • La función $f(x)=\cos(x)$
    • El polinomio de aproximación obtenido.
  • Calcular el error relativo $$E_r = \frac{\left\Vert \cos(x)-P(x)\right\Vert}{\left\Vert \cos(x)\right\Vert}$$ usando la función numpy.linalg.norm.
In [10]:
%run Ejercicio5.py

Er =  0.00388705029515

Ejercicio 6

Dada la función $f(x)=\cos(\arctan(x))-e^{x^2}\log(x+2)$, calcular la aproximación polinómica de $f$ de grado 4. Dar el error relativo y dibujar la función y la aproximación polinómica como en el ejemplo anterior.

In [11]:
%run Ejercicio6.py

Er =  0.0346826444554

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$.

Ejercicio 7

Dibujar la serie de Fourier de orden 5 para la función $f(x)=x$ en el intervalo $[0,3]$ de periodo $T=3$.

In [12]:
%run Ejercicio7.py

Ejercicio 8

Dibujar la serie de Fourier de orden 6 asociada a la función $f(x)=(x-\pi)^2$ en el intervalo $[0,2\pi]$ de periodo $T=2\pi$.

In [13]:
%run Ejercicio8.py