# Sistemas de ecuaciones lineales

...pero antes
------------


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:


```python
f = lambda x: x**2 - 2*x + 1
```

e importamos la función ```roots``` de la libreria ```numpy```:


```python
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:


```python
roots?
```

De la documentación podemos ver que la forma de utilizarla es dandole un arreglo de numeros, con los coeficientes del polinomio: 


```python
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$:


```python
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:


```python
%%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:


```python
from numpy import sin
```


```python
g = lambda x: sin(x) + 1
```

importamos la función ```fsolve``` de la libreria de optimización de ```scipy```:


```python
from scipy.optimize import fsolve
```

y la utilizamos dandole la función a resolver, y la aproximación o valor inicial:


```python
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:


```python
fsolve(g, [0, 5, 10, 15])
```

Este metodo tambien funcionará para ecuaciones polinómicas:


```python
fsolve(f, 0)
```

Sin embargo va a ser considerablemente mas lento que las funciones dedicadas para ello:


```python
%%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.

Matrices
-------

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```:


```python
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}
$$


```python
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$


```python
A.I
```


```python
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:


```python
%%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```:


```python
from numpy.linalg import solve
```

y se utiliza dandole la matriz $A$ y el vector $b$:


```python
solve(A, b)
```

al revisar el tiempo que tomará esta solución, podemos ver una mejoria significativa:


```python
%%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:


```python
from numpy.linalg import det
```


```python
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):


```python
from numpy import hstack
```

La primer matriz a obtener es una en la cual la primer columna sea $b$:


```python
A1 = hstack((b, A[:,1], A[:,2]))
A1
```

la segunda en donde la segunda columna sea $b$:


```python
A2 = hstack((A[:,0], b, A[:,2]))
A2
```

y asi:


```python
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:


```python
%%timeit
x1 = det(A1)/det(A)
x2 = det(A2)/det(A)
x3 = det(A3)/det(A)
```


```python
x1, x2, x3
```


```python
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:

Problemas
=====

1. Implementa una función, que dadas la matrix $A$ y el vector $b$, obtenga la solución de un sistema de ecuaciones lineal (utiliza el método de Cramer).
2. Obten la solución de el siguiente sistema de ecuaciones lineal:
$$
\begin{align}
x_1 + 2x_2 + 1x_3 &= 2 \\
4x_1 + 5x_2 + 10x_3 &= -1 \\
8x_1 - x_2 + 2x_3 &= 10
\end{align}
$$
3. La solución obtenida con este método, ¿Es más rápida que el método de la inversa? ¿Es más rápida que la función ```solve```? Implementa el código para checar estos datos.
4. (Opcional) ¿Que problemas puede tener el método de Cramer? ¿Que código es necesario agregar para checar este problema?

a05e05ab720086da4cc0aa8e3bdc65fecb728423
ipynb
Jupyter Notebook
Practicas/P3/Practica 3 - Sistemas de ecuaciones lineales.ipynb
robblack007/clase-metodos-numericos
f16a11b4bba6b1927ee47677f000e8fa53f8b39b
MIT
Practicas/P3/Practica 3 - Sistemas de ecuaciones lineales.ipynb
robblack007/clase-metodos-numericos
f16a11b4bba6b1927ee47677f000e8fa53f8b39b
MIT
Practicas/P3/Practica 3 - Sistemas de ecuaciones lineales.ipynb
robblack007/clase-metodos-numericos
f16a11b4bba6b1927ee47677f000e8fa53f8b39b
MIT
Qwen/Qwen-72B
1. YES
2. YES