# Álgebra Lineal Numérica

El álgebra lineal es el área de las matemáticas que estudia los espacios vectoriales y las transformaciones lineales entre dichos espacios. Independientemente del espacio vectorial que trabajemos, mientras sea de dimensión finita, todos elementos del espacio se pueden representar como un **vector** (una $n$-ada de elementos del campo del espacio) y toda transofmación lineal se representa como una **matriz** (una colección de $n \times n$ elementos del campo del espacio)

## El problema $\mathbf{Ax}=\mathbf{b}$ y sus equivalencias

En matemáticas y física, muchas veces surgen problemas en los que tenemos que resolver un sistema de ecuaciones lineales de la siguiente forma:

$$
\begin{align}
a_{11} x_1 + a_{12} x_2 + \ldots + a_{1n} x_n &= b_1 \\
a_{21} x_1 + a_{22} x_2 + \ldots + a_{2n} x_n &= b_2 \\
\vdots \quad \quad \quad \vdots \quad \quad \quad \vdots \\
a_{n1} x_1 + a_{n2} x_2 + \ldots + a_{nn} x_n &= b_n \\
\end{align}
$$

Que se pueden representar, de forma matricial, como

$$
 \begin{pmatrix}
  a_{1,1} & a_{1,2} & \cdots & a_{1,n} \\
  a_{2,1} & a_{2,2} & \cdots & a_{2,n} \\
  \vdots  & \vdots  & \ddots & \vdots  \\
  a_{n,1} & a_{n,2} & \cdots & a_{n,n} 
 \end{pmatrix} 
 \begin{pmatrix}
 x_1 \\
 x_2 \\
 \vdots \\
 x_n
 \end{pmatrix} 
 =
 \begin{pmatrix}
 b_1 \\
 b_2 \\
 \vdots \\
 b_n
 \end{pmatrix} 
$$


Si definimos

$$
\mathbf{A} = 
\begin{pmatrix}
  a_{1,1} & a_{1,2} & \cdots & a_{1,n} \\
  a_{2,1} & a_{2,2} & \cdots & a_{2,n} \\
  \vdots  & \vdots  & \ddots & \vdots  \\
  a_{n,1} & a_{n,2} & \cdots & a_{n,n} 
 \end{pmatrix} \quad , \quad \mathbf{x} = \begin{pmatrix}
 x_1 \\
 x_2 \\
 \vdots \\
 x_n
 \end{pmatrix} \quad \text{y} \quad \mathbf{b}= \begin{pmatrix}
 b_1 \\
 b_2 \\
 \vdots \\
 b_n
 \end{pmatrix} $$
 
 Podemos escribir el sistema de manera compacta:
 
 $$\mathbf{Ax} = \mathbf{b}$$

### Equivalencias del problema

Muchos problemas de álgebra lineal se pueden plantear como problemas de resolver un sistema de ecuaciones lineales. Por ejemplo:

* Encontrar la inversa de una matriz, es decir, encontrar $\mathbf{B}$ tal que $\mathbf{A}\mathbf{B} = \mathbf{Id}$

* Encontrar intersecciones de rectas o planos en el espacio

* Encontrar matrices de cambio de base.

* Encontrar los coeficientes una transformación lineal aplicada en $n$ vectores linealmente independientes

* Resolver ecuaciones diferenciales (parciales y ordinarias) discretizadas o numéricas.

### Sistemas de ecuaciones lineales problemáticos

En principio, no necesariamente se da el caso donde la matriz $A$ es de $n\times n$, puede suceder que tengamos menos o más ecuaciones que incógnitas. Tratar con esos casos es **complicado**. Así, siempre vamos a realizar la suposición de que *la matriz es de $n\times n$*.

Aunque nuestra matriz $A$ sea de $n\times n$, puede darse el caso de que nuestro sistema de ecuaciones lineales **no acepte una solución única**. Nuevamente, por simplicidad, siempre asumiremos que nuestro sistema admite una solución **única**. Condiciones equivalentes para eso son:

* $\mathbf{A}$ es invertible

* $\text{det}(\mathbf{A}) \neq 0$ 

* $\text{rango}(\mathbf{A}) = n$

## ¿Cómo resolvemos un sistema lineal? 

Existen varias maneras analíticas y exactas de resolver un sistema de ecuaciones lineales (Regla de Kraemer, Eliminación Gaussiana, Sustitución, etc). Sin embargo, la gran mayoría de ellas son **tediosas**, **mecánicas** y largas de realizar manualmente. En particular, cuando tenemos más de $4$ ecuaciones, se vuelve extremadamente tedioso resolverlas.

Debido a su tediosidad y a su importancia en todas las ramas, resulta ideal utilizar la computadora para resolver dichos sistemas. 

## Preludio computacional: vectores y matrices en Julia

En Julia podemos trabajar de manera nativa con matrices y vectores. Es claro que un **arreglo de flotantes** puede interpretarse como un **vector** por ser un conjunto ordenado de elementos del mismo tipo.

Dentro de Julia, es posible operar a los arreglos entre sí de la misma manera que operamos a los vectores: con suma y producto por un escalar


```julia
vec1 = [1.0,1.0,0.0]
vec2 = [1.0,0.0,1.0]
print("suma: ")
println(vec1+vec2)
alpha = 5.0
print("producto por un escalar: ")
println(alpha*vec1)
```

    suma: [2.0, 1.0, 1.0]
    producto por un escalar: [5.0, 5.0, 0.0]


Para representar a las matrices, una de las ventajas de Julia es que podemos construir **arreglos $n$-dimensionales**. Por ahora, nos basta con construir arreglos $2$ dimensionales. La sintaxis para construir un arreglo 2 dimensional de manera explícita es la siguiente:


```julia
# ponemos los renglones como arreglos 1D separados por ESPACIOS, no por comas
# cada renglon se separa entre sí por un punto y coma
A = [[1.0 2.0 3.0] ; [4.0 5.0 6.0] ]
```




    2×3 Array{Float64,2}:
     1.0  2.0  3.0
     4.0  5.0  6.0



Podemos concatenar también vectores de manera horizontal (i.e. como si fueran columnas) para construir una matriz, utilizando la siguiente sintaxis


```julia
A2 = [vec1 vec2]
```




    3×2 Array{Float64,2}:
     1.0  1.0
     1.0  0.0
     0.0  1.0



Es posible transponer una matriz utilizando la función `transpose`


```julia
transpose(A2)
```




    2×3 LinearAlgebra.Transpose{Float64,Array{Float64,2}}:
     1.0  1.0  0.0
     1.0  0.0  1.0



Si ahora nos preguntamos cuál es el tipo de una matriz, podemos usar la función `typeof`


```julia
typeof(A)
```




    Array{Float64,2}



Notemos que ahora la función `typeof(A)` nos muestra un `2` después de `Float64`, lo que nos indica que nuestro objeto `A` es un arreglo de dimensión 2. Para los vectores normales, utilizamos arreglos de dimensión $1$. Las dimensiones del arreglo se pueden consultar con la función `size`


```julia
size(A2)
```




    (3, 2)



Podemos multiplicar una matrix por un vector si coinciden las dimensiones


```julia
# matriz de 3x3
A4 = [[1 2 3] ; [2 4 5] ; [7 8 9]]
# vector de 3 entradas
vec = [1,0,0]
# su multiplicación
println(A4*vec)
```

    [1, 2, 7]



```julia
# matriz de 4x3
A5 = [[0 1 4] ; [1 2 3] ; [2 4 5] ; [7 8 9]]
# vector de 3 entradas
vec = [1,0,0]
# su multiplicación
println(A5*vec)
```

    [0, 1, 2, 7]


Si no coinciden las columnas de la matriz con la longitud del vector, intentar la multiplicación nos arrojará un error:


```julia
# matriz de 3x2
A6 = [[1 2] ; [2 4] ; [7 8]]
# vector de 3 entradas
vec = [1,0,0]
# no está definida la multiplicación
println(A6*vec)
```

Si le pedimos a Julia que nos imprima una matriz, el formato es bastante feo


```julia
# impresión fea
println(A6)
```

    [1 2; 2 4; 7 8]


Podemos utilizar la función `display(A)` para imprimirlo como una matriz. La sintaxis es la siguiente:


```julia
# impresión bonita
display(A6)
```


    3×2 Array{Int64,2}:
     1  2
     2  4
     7  8


Para acceder a los elementos de una matriz, utilizamos la sintaxis `A[i,j]`


```julia
display(A)
println(A[1,2])
println(A[2,3])
```


    2×3 Array{Float64,2}:
     1.0  2.0  3.0
     4.0  5.0  6.0


    2.0
    6.0


Podemos acceder a renglones o columnas completos utilizando `:` en lugar de `i` o `j`, respectivamente


```julia
# primera renglón
println(A[1,:])
# tercer columna
println(A[:,3])
```

    [1.0, 2.0, 3.0]
    [3.0, 6.0]


Notemos que en julia podemos utilizar una rango o arreglo de enteros como índice de otro arreglo para obtener el subarreglo con los elementos de dichos índices. Por ejemplo:


```julia
arr1 = [2,4,8,16,32,128]
println(arr1[[2,5]])
println(arr1[1:3])
```

    [4, 32]
    [2, 4, 8]


Lo mismo sucede con las matrices: podemos obtener una submatriz de esa forma


```julia
mat1 = [[1 2 3 4] ; [5 6 7 8]]
display(mat1)
display(mat1[1:2,3:4])
display(mat1[1:2,2:4])
```


    2×4 Array{Int64,2}:
     1  2  3  4
     5  6  7  8



    2×2 Array{Int64,2}:
     3  4
     7  8



    2×3 Array{Int64,2}:
     2  3  4
     6  7  8


## Caso sencillo: $A$ es triangular superior

Regresando a querer resolver el problema $\mathbf{Ax} = \mathbf{b}$

Antes que atacar el problema de manera general, consideremos el caso sencillo en el que la matriz $\mathbf{A}$ 
 es una matriz triangular superior, es decir, de la forma.
 
 $$
 \mathbf{A} = 
\begin{pmatrix}
  a_{1,1} & a_{1,2} & \cdots & a_{1,n} \\
  0 & a_{2,2} & \cdots & a_{2,n} \\
  \vdots  & \vdots  & \ddots & \vdots  \\
   0 & \ldots  & a_{n-1,n-1} & a_{n-1,n} \\
 0 & \ldots  & 0 & a_{n,n} \\
 \end{pmatrix}
 $$
 
 Matemáticamente, eso quiere decir que para $i>j$, $A_{ij} = 0$

Cuando la matriz es así, es muy sencillo resolver el sistema de ecuaciones.

$$
\begin{pmatrix}
  a_{1,1} & a_{1,2} & \cdots & a_{1,n} \\
  0 & a_{2,2} & \cdots & a_{2,n} \\
  \vdots  & \vdots  & \ddots & \vdots  \\
   0 & \ldots  & a_{n-1,n-1} & a_{n-1,n} \\
 0 & \ldots  & 0 & a_{n,n} \\
 \end{pmatrix}
 \begin{pmatrix}
 x_1 \\
 x_2 \\
 \vdots \\
 x_{n-1} \\
 x_n
 \end{pmatrix} 
 =
 \begin{pmatrix}
 b_1 \\
 b_2 \\
 \vdots \\
 b_{n-1} \\
 b_n
 \end{pmatrix}
$$ 

### Ejercicio 1

(i) Suponiendo que $\mathbf{A}$ es triangular superior y que no tiene ceros en la diagonal, encuentra la expresión (analítica) para $x_n$ en términos de las entradas de $A$ y las entradas del vector solución $\mathbf{b}$. 

(ii) Ya que tengas la expresión para $x_n$, ahora encuentra una expresión para $x_{n-1}$ (que también debe de depender del valor de $x_n$) y para $x_{n-2}$ (que depende también del valor de $x_n$ y $x_{n-1}$ )

(iii) Generaliza las expresiones anteriores para obtener una expresión general para cualquier $x_i$ en términos de las entradas de $A$, las entradas del vector solución $\mathbf{b}$ y de otros valores $x_k$ con $k>i$ 

### Ejercicio 2

Implementa una función `solTriSuperior(A,b)` que tome como argumento un arreglo 2D `A`, que representa a una matriz **triangular superior, sin ceros en la diagonal**, de $n\times n$ , un vector solución `b` de longitud $n$ y que regresa un arreglo `xs` con la solución del sistema $\mathbf{Ax} = \mathbf{b}$ calculada usando la fórmula del inciso (iii) del ejercicio anterior.

## Caso general: volver a $\mathbf{A}$ una matriz triangular superior

Generalmente, $\mathbf{A}$ no va a ser una matriz triangular superior. Sin embargo, mediante el procedimiento de [**eliminación gaussiana**](https://en.wikipedia.org/wiki/Gaussian_elimination) es posible llevarla a una de ellas. 

### Eliminación Gaussiana

La eliminación Gaussiana consiste en trabajar con la **matriz aumentada**, denotada $\mathbf{C}$, cuya expresión es

$$
\mathbf{C} = 
\begin{pmatrix}
  a_{1,1} & a_{1,2} & \cdots & a_{1,n} & b_1 \\
  a_{2,1} & a_{2,2} & \cdots & a_{2,n} & b_2 \\
  \vdots  &   & \ddots &   & \vdots \\
  a_{n,1} & a_{n,2} & \cdots & a_{n,n} & b_n 
 \end{pmatrix} 
 $$

Mediante operaciones elementales de renglones (intercambiar renglones, multiplicarlos por un escalar o sumarlos), buscamos reducir la matriz $\mathbf{C}$ a una forma **escalonada**, denotada $\mathbf{C}^*$, en la que todos los elementos por debajo de la diagonal sean $0$:

$$
\mathbf{C}^* = 
\begin{pmatrix}
  a^*_{1,1} & a^*_{1,2} & \cdots & a^*_{1,n} & b^*_1 \\
  0 & a^*_{2,2} & \cdots & a^*_{2,n} & b^*_2 \\
  \vdots  &   & \ddots &   & \vdots \\
  0 & 0 & \cdots & a^*_{n,n} & b^*_n 
 \end{pmatrix} 
 $$
 
 Esa matriz aumenta representa al sistema de ecuaciones:
 
 $$
\begin{pmatrix}
  a^*_{1,1} & a^*_{1,2} & \cdots & a^*_{1,n} \\
  0 & a^*_{2,2} & \cdots & a^*_{2,n} \\
  \vdots  &   & \ddots &   \vdots \\
  0 & 0 & \cdots & a^*_{n,n} 
 \end{pmatrix} 
 \begin{pmatrix}
x_1 \\
x_2 \\
\vdots \\
x_n 
 \end{pmatrix}=
 \begin{pmatrix}
b^*_1 \\
b^*_2 \\
\vdots \\
b^*_n 
 \end{pmatrix}
 $$
 

El sistema de ecuaciones representado por $\mathbf{C}^*$ es equivalente al representado por $\mathbf{C}$, por lo que sus soluciones son las mismas. Sin embargo, el sistema asociado a la matriz escalonada es triangular superior, por lo que es mucho más fácil de resolver, como vimos en los ejercicios anteriores.

### Ejercicio 3:

(i) Supón que $a_{1,1} \neq 0$. ¿Qué operación elemental podrías aplicarle al renglón $n$ para hacer que el elemento $a_{n,1}$ se vuelva 0?

**Sugerencia** La operación es convertir al renglón $n$ en una combinación lineal del renglón $1$ y el renglón $n$

(ii) ¿Puedes aplicar el procedimiento un procedimiento similar al del inciso anterior pero ahora para hacer $a_{n-1,1}=0$? ¿y para cualquier otro $a_{k,1}=0$ con $k \neq 1$?

(iii) Vamos ahora a la segunda columna. Supón que ya realizaste todas las operaciones elementales necesarias para que $a_{k,1} =0$ para $k \neq 1$. Nuevamente suponiendo que $a_{2,2}\neq 0$, ¿cómo puedes hacer algo parecido a las operaciones elementales de los incisos anteriores para volver $a_{k,2} =0$ para $k \neq 2$?

(iv) Generaliza todo lo visto anteriormente para encontrar las operaciones elementales necesarias (y el orden de ellas) para convertir a la matriz aumentada $\mathbf{C}$ en la matriz escalonada $\mathbf{C}^*$.

### Ejercicio 4

Implementa una función `elimGaussBasica(A,b)` que toma como argumento un arreglo 2D `A` de $n\times n$ **suponiendo que no tiene ceros en la diagonal** y un arreglo 1D `b` de longitud $n$. La función debe de construir la matriz aumentada $\mathbf{C}$ obtenida a partir de `A` y `b` y luego realizar las operaciones elementales obtenidas en el ejercicio anterior para generar la matriz escalonada $\mathbf{C}^*$. La función debe de regresar la matriz escalonada obtenida de $n\times (n+1)$.

**Sugerencia:** la función `hcat` te permite concatenar vectores y matrices horizontalmente de la siguiente forma:


```julia
A = [[1 2 3] ; [4 5 6]]
B = [8,9]
display(A)
display(B)
display(hcat(A,B))
```


    2×3 Array{Int64,2}:
     1  2  3
     4  5  6



    2-element Array{Int64,1}:
     8
     9



    2×4 Array{Int64,2}:
     1  2  3  8
     4  5  6  9


## ¿Qué sucede si hay elementos $a_{ii}=0$?

Es fundamental que los elementos de la diagonal no sean $0$ para que nuestro método funcione. Sin embargo, es posible que se de el caso en el que sean $0$. Si $a_{k,k}=0$, podemos buscar un renglón $l$ en el que $a_{l,k}\neq 0 $ y luego **intercambiar** el renglón $k$ por el $l$.

Debido a que estamos suponiendo que la solución al sistema es única, no puede nunca darse el caso de $a_{l,k}=0$ para toda $l$ pues eso implicaría que no existe solución única para el sistema

### Ejercicio 5

Implementa una función `checarDiagonal(A,b)` que toma como argumento un arreglo 2D `A` de $n\times n$ **suponiendo que no tiene ceros en la diagonal** y un arreglo 1D `b` de longitud $n$. La función debe de construir la matriz aumentada $\mathbf{C}$ obtenida a partir de `A` y `b` y luego revisar sus elementos diagonales para ver que no sean $0$. Si encuentra alguno que sea $0$, debe de intercambiar el renglón por otro en el que sea distinto de cero. La función debe de regresar una matriz aumentada  en la que no haya ceros en las diagonales.

### Ejercicio 6

Implementa una función `eliminaciónGaussiana(A,b)` que toma como argumento un arreglo 2D `A` de $n\times n$ **suponiendo que no tiene ceros en la diagonal** y un arreglo 1D `b` de longitud $n$. La función primero debe de usar `checarDiagonal(A,b)` para obtener una matriz aumentada sin ceros en las diagonales y luego debe de aplicar las operaciones elementales del ejercicio 3 para obtener una matriz escalonada $\mathbf{C}^*$. La función debe de regresar dicha matriz escalonada de $n\times (n+1)$.

### Ejercicio 7

Implementa una función `ecLineales(A,b)` que toma como argumento un arreglo 2D `A` de $n\times n$  y un arreglo 1D `b` de longitud $n$. Tu función debe de utilizar las funciones `eliminaciónGaussiana(A,b)` y `solTriSuperior(A,b)` para resolver el sistema de ecuaciones lineales y regresar un arreglo 1D de longitud $n$ correspondiente a las soluciones $x_i$.

### Ejercicio 8

Prueba tu función `ecLineales(A,b)` con el siguiente sistema

$$
\begin{align}
    x+y+z &= 1 \\
    3x-2y+w &= -4 \\
    y - w &= 2 \\
    x-2y+4z-5w &= -6
\end{align}
$$

Cuya solución está dada por $x= -0.0555$, $y=1.8333$, $z=-0.7777$ y $-w=0.1666$

9941f8a5a2310b67c5d91c4fac1a722258c4d27e
ipynb
Jupyter Notebook
files/fiscomp_2020-4/material/clase11.ipynb
sayeg84/sayeg84.github.io
18f2e36dd7252603fad8f7093dc5aa00fc721be4
MIT
files/fiscomp_2020-4/material/clase11.ipynb
sayeg84/sayeg84.github.io
18f2e36dd7252603fad8f7093dc5aa00fc721be4
MIT
files/fiscomp_2020-4/material/clase11.ipynb
sayeg84/sayeg84.github.io
18f2e36dd7252603fad8f7093dc5aa00fc721be4
MIT
Qwen/Qwen-72B
1. YES
2. YES