Antes de empezar con el tema de la práctica de hoy, quiero hacer notar un detalle que solo he hecho notar teoricamente, pero del cual no he ofrecido evidencia empirica, empecemos nombrando algunos métodos numéricos que podemos utilizar, por ejemplo, para obtener las raices de un polinomio, podemos utilizar la función roots:
Antes que nada, definimos el polinomio del cual queremos obtener las raices reales:
In [ ]:
f = lambda x: x**2 - 2*x + 1
e importamos la función roots de la libreria numpy:
In [ ]:
from numpy import roots
Si no estamos seguros de como utilizar la función, podemos obtener la documentación rapidamente al postponer un signo de interrogación, ?, al nombre de la función y ejecutarel código; este aparecera en una ventana en la parte inferior y se puede cerrar en cuanto hayas terminado de leerla o reducir el espacio que ocupa en tu navegador:
In [ ]:
roots?
De la documentación podemos ver que la forma de utilizarla es dandole un arreglo de numeros, con los coeficientes del polinomio:
In [ ]:
roots([1, -2, 1])
y de la misma manera, podemos asegurarnos de que estas raices son adecuadas, al sustituir este valor en la función y comprobar que el valor que nos devuelve es $0$:
In [ ]:
f(1)
Ahora quiero hacer notar que podemos revisar cuanto tiempo le toma a la computadora hacer esta operación, al utilizar la función %%timeit al inicio de la celda de código, el motor de Jupyter correra la función $1000$ veces y nos reportará el tiempo promedio que le tomó hacer esta operación:
In [ ]:
%%timeit
rs = roots([1, -2, 1])
Sin embargo esta función no puede darnos la solución para ecuaciones no lineales, para esto podemos utilizar fsolve, primero definamos una ecuación no lineal:
In [ ]:
from numpy import sin
In [ ]:
g = lambda x: sin(x) + 1
importamos la función fsolve de la libreria de optimización de scipy:
In [ ]:
from scipy.optimize import fsolve
y la utilizamos dandole la función a resolver, y la aproximación o valor inicial:
In [ ]:
fsolve(g, 0)
Cabe notar que a esta función podemos darle mas de una aproximación inicial y me dará el resultado de el algoritmo para cada una de las aproximaciones lineales:
In [ ]:
fsolve(g, [0, 5, 10, 15])
Este metodo tambien funcionará para ecuaciones polinómicas:
In [ ]:
fsolve(f, 0)
Sin embargo va a ser considerablemente mas lento que las funciones dedicadas para ello:
In [ ]:
%%timeit
rs = fsolve(f, 0)
La función roots, utiliza un algoritmo diferente basado en el binomio de Newton, por lo que solo sirve para las ecuación polinomiales, y la función fsolve utiliza métodos numéricos parecidos a los que hemos estado estudiando, pero de una sofisticación mayor, por lo que la moraleja de esta historia es:
Las funciones programadas en las librerías van a ser mas rápidas siempre, pero hay que saber cuando utilizar cada una y cuando pueden fallar.
Dentro de la libreria de calculo numérico numpy, existe una definición que utilizaremos para utilizar matrices como elementos de nuestro ambiente de programación, primero tenemos que importar la definición de matrix:
In [ ]:
from numpy import matrix
para utilizarla, tengo que darle una lista que contenga una lista por cada fila de la matriz a definir, empecemos definiendo la matriz para el siguiente sistema de ecuaciones lineales:
$$ \begin{align} x_1 + 4x_2 + 10x_3 &= 5 \\ 2x_1 + 3x_2 + 5x_3 &= 2 \\ x_1 + 5x_2 + 2x_3 &= 10 \end{align} $$para el cual, podemos definir $A$ y $b$, de tal manera que:
$$ Ax = b $$en donde:
$$ x= \begin{pmatrix} x_1 \\ x_2 \\ x_3 \end{pmatrix} \quad A= \begin{pmatrix} 1 & 4 & 10 \\ 2 & 3 & 5 \\ 1 & 5 & 2 \end{pmatrix} \quad b = \begin{pmatrix} 5 \\ 2 \\ 10 \end{pmatrix} $$
In [ ]:
A = matrix([[1, 4, 10],
[2, 3, 5],
[1, 5, 2]])
b = matrix([[5],
[2],
[10]])
Por lo que la solución mas simple, en terminos de código, es calcular la matriz inversa y multiplicarla por el vector $b$
In [ ]:
A.I
In [ ]:
x = A.I*b
x
Sin embargo, ya que hablamos de la rapidez de los métodos numéricos, chequemos la velocidad de esta solución:
In [ ]:
%%timeit
x = A.I*b
En este caso, tenemos un método alternativo el cual utiliza algoritmos iterativos para el cálculo de esta solución; primero tenemos que importar la función solve de la librería de algebra lineal de numpy:
In [ ]:
from numpy.linalg import solve
y se utiliza dandole la matriz $A$ y el vector $b$:
In [ ]:
solve(A, b)
al revisar el tiempo que tomará esta solución, podemos ver una mejoria significativa:
In [ ]:
%%timeit
xs = solve(A, b)
Por ultimo, veamos el método de Cramer, podemos implementarlo facilmente, importando la función para calcular el determinante de una matriz:
In [ ]:
from numpy.linalg import det
In [ ]:
det(A)
de aqui sabemos que el determinante general tiene un valor no nulo, tambien necesitaremos sustituir el valor del vector $b$ en columnas especificas de la matriz $A$, haremos esto obteniendo la columna especifica y apilandolas horizontalmente, de ahi el nombre de la función que necesitamos hstack (horizontal stack):
In [ ]:
from numpy import hstack
La primer matriz a obtener es una en la cual la primer columna sea $b$:
In [ ]:
A1 = hstack((b, A[:,1], A[:,2]))
A1
la segunda en donde la segunda columna sea $b$:
In [ ]:
A2 = hstack((A[:,0], b, A[:,2]))
A2
y asi:
In [ ]:
A3 = hstack((A[:,0], A[:,1], b))
A3
con lo que solo tenemos que calcular el determinante de cada una de estas matrices entre el determinante general:
In [ ]:
%%timeit
x1 = det(A1)/det(A)
x2 = det(A2)/det(A)
x3 = det(A3)/det(A)
In [ ]:
x1, x2, x3
In [ ]:
det([[1,2,1], [4,5,10], [8, -1, 2]])
y como puedes ver, obtenemos los mismos valores que en los otros dos métodos:
solve? Implementa el código para checar estos datos.