Sistemas de ecuaciones lineales

Introducción

Esta sección da herramientas para resolver los problemas del tipo $A x = b$.

Discretización

Para los sistemas que se modelan con la ecuación de Laplace, si hacemos algunos manipuleos con Taylor ($f(x + h) + f(x - h)$) llegamos a que se pueden discretizar así:

(1)
\begin{align} u''(x,y) = \frac{u(x-h, y) + u(x+h, y) + u(x,y-h) + u(x,y+h) - 4u(x,y)}{h^2} \end{align}

donde $h$ es el espaciamiento entre los puntos discretizados.

Para una sola variable esto se convierte en

(2)
\begin{align} y''(x) = \frac{y(x-h) +y(x+h)+y(x)+y(x) - 4y(x)}{h^2} = \frac{y(x-h)+y(x+h)-2y(x)}{h^2} \end{align}

más algún término de error.

Pasando a los nodos en el caso de una variable queda

(3)
\begin{align} y''(x_i) = \frac{y_{i-1} +y_{i+1} - 2y_i}{h^2} \end{align}

entonces de acá extraemos un sistema de ecuaciones (necesitamos condiciones de borde).

Método de eliminación gaussiana (GEM)

Este método transforma un sistema lineal en otro equivalente (con las mismas soluciones) para resolverlo más fácilmente. También se lo llama Gauss-Jordan.

Se hace:

  • Multiplicando una fila por una constante.
  • Sumando una ecuación con el múltiplo de otra.
  • Intercambiando ecuaciones.

De esta manera se debe obtener una matriz triangular superior o inferior, logrando despejar fácilmente una solución para luego reemplazarla e ir obteniendo las demás.

A tener en cuenta: si el determinante de la matriz es distinto de cero, el sistema tiene solución única.

Pivoteo parcial

La onda con el pivoteo es hacer Gauss-Jordan pero eligiendo el mejor valor para hacer cuentas, de manera tal que los errores de las cuentas de la maquinita sean los menores. Así, se implementa un criterio para reordenar filas de manera tal que se queden los valores mayores al momento de dividir. El método consiste en elegir en la iteración i-ésima, el elemento de la columna que sea el de mayor valor absoluto entre todos los que están por encima de la diagonal, el menor valor de $p\geq i$ tal que

(4)
\begin{align} |a_{pi}| = \max\limits_{i\leq k\leq n}|a_{ki}| \end{align}

para poder hacer el intercambio Fila $i$ por Fila $p$.

Ejemplo

La matriz:

(5)
\begin{align} \left( \begin{array}{rrr|r} 5 & 4 & 2 & 1 \\ 4 & 1 & 3 & 0 \\ 9 & 7 & 1 & 2 \end{array} \right) \end{align}

quedaría:

(6)
\begin{align} \left( \begin{array}{rrr|r} 9 & 7 & 1 & 2 \\ 5 & 4 & 2 & 1 \\ 4 & 1 & 3 & 0 \end{array} \right) \end{align}

para el primer paso en que se logren los ceros inferiores en la primera columna. Luego se procedería a hacer lo mismo para las siguientes dos filas.

Pivoteo parcial escalado

Cuando hay algún valor muy zarpado en las ecuaciones en comparación con el resto, el pivoteo solito no alcanza. Ahí hay que usar escalado. El escalado implica usar un factor de escala para cada fila, que se define

(7)
\begin{align} s_k = \max\limits_{1\leq j\leq n}|a_{kj}| \end{align}

Ahora, antes de eliminar la variable $x_i$ el intercambio de filas Filai <-> Filap se hace tomando el primer entero $p \ne i$ tal que

(8)
\begin{align} \frac{|a_{pi}|}{s_p}=\max\limits_{i\leq k\leq n}\frac{|a_{ki}|}{s_k} \end{align}

El pivoteo parcial agrega 3/2(n²-n) comparaciones y (n²+n)/2 -1 divisiones, por lo que para lograr una precisión mayor hace falta más poder de cómputo.

Ejemplo:

(9)
\begin{align} 2.11 x_1 &- 4.21 x_2 + 0.921 x_3 & =2.01\\ 4.01 x_1 &+10.2 x_2 - 1.12 x_3 & = -3.09\\ 1.09 x_1 &+ 0.987 x_2 + 0.832 x_3 & = 4.21 \end{align}

Este sistema que es de $n = 3$ tiene como matriz a

(10)
\begin{align} \left( \begin{array}{rrr|r} 2.11 & -4.21 & 0.921 & 2.01 \\ 4.01 & 10.2 & -1.12 & -3.09 \\ 1.09 & 0.987 & 0.831 & 4.21 \end{array} \right) \end{align}

Luego buscamos los $s_k$

(11)
\begin{eqnarray} s_1 \text{m\'aximo de la fila 1} = \max\limits_{1\leq j \leq n} |a_{1j}| = 4.21 \\ s_2 \text{m\'aximo de la fila 2} = \max\limits_{1\leq j \leq n} |a_{2j}| = 10.2 \\ s_3 \text{m\'aximo de la fila 3} = \max\limits_{1\leq j \leq n} |a_{3j}| = 1.09 \end{eqnarray}

Luego queda buscar el pivote, para ello se utiliza la ecuación 8 que busca el $\max\limits_{1\leq k\ n}\frac{|a_{k1}|}{s_k}$

(12)
\begin{align} \frac{|a_{11}|}{s_1}= \frac{2.11}{4.21} = 0.501 \text{,}\quad \frac{|a_{21}|}{s_2}= \frac{4.01}{10.2} = 0.393 \quad\text{y}\quad \frac{|a_{31}|}{s_3}= \frac{1.09}{1.09} = 1 \end{align}

De ahí agarramos el mayor que es $\frac{|a_{31}|}{s_3}$ y se hace el intercambio Fila1 <-> Fila3. La matriz 10 pasa a quedar

(13)
\begin{align} \left( \begin{array}{rrr|r} 1.09 & 0.987 & 0.831 & 4.21 \\ 4.01 & 10.2 & -1.12 & -3.09 \\ 2.11 & -4.21 & 0.921 & 2.01 \end{array} \right) \end{align}

Los multiplicadores quedan $m_{21} = 3.68$ y $m_{31}=1.94$ y eliminando, 13 queda como

(14)
\begin{align} \left( \begin{array}{rrr|r} 1.09 & 0.987 & 0.831 & 4.21 \\ 0 & 6.57 & -4.18 & -18.6 \\ 0 & -6.12 & -0.689 & -6.16 \end{array} \right) \end{align}

Repetimos una vez más la búsqueda de los máximos, manteniendo la escala que usamos antes. Como

(15)
\begin{align} \frac{|a_{22}|}{s_2} = \frac{6.57}{10.2} = 0.644 < \frac{|a_{32}|}{s_3} = \frac{6.12}{1.09} = 5.61 \end{align}

El pivote pasa a ser el de la fila 3, por lo que se intercambian Fila2 <-> Fila3 y queda

(16)
\begin{align} \left( \begin{array}{rrr|r} 1.09 & 0.987 & 0.831 & 4.21 \\ 0 & -6.12 & -0.689 & -6.16 \\ 0 & 6.57 & -4.18 & -18.6 \end{array} \right) \end{align}

y el $m_{32}=\frac{6.57}{-6.12}=-1.07$. Eliminando queda

(17)
\begin{align} \left( \begin{array}{rrr|r} 1.09 & 0.987 & 0.831 & 4.21 \\ 0 & -6.12 & -0.689 & -6.16 \\ 0 & 0 & -4.92 & -25.2 \end{array} \right) \end{align}

Sustituyendo hacia atrás: $x = (-0.431\; 0.430\; 5.12)$.

Factorización LU

La idea es factorizar a la matriz $A$ en una triangular superior y otra triangular inferior ($L$=lower, $U$=upper). Esto permite resolver la ecuación $A x = b$ con orden N (una vez factorizado $A$). Como factorizar $A$ tiene el mismo orden que gauss, esto se suele hacer cuando $b$ varía.

Haciendo GEM sobre la matriz original obtenemos una matriz inferior a la que denominamos $U$. Después convertimos las operaciones elementales que utilizamos para obtener $U$ y las convertimos en matrices (a las cuales llamamos $E_i$). Luego invertimos las matrices $E_i$ y las multiplicamos en orden contrario (por lo tanto se supone que no reordenamos filas para hacer LU).

Por ejemplo si tenemos el sistema

(18)
\begin{equation} Ax = b \end{equation}

y para obtener $U$ tubimos tuvimos que multiplicar a $A$ por $E_1$ $E_2$ y $E_3$
nos queda

(19)
\begin{equation} E_3 E_2 E_1 Ax = E_3 E_2 E_1 b \end{equation}

donde llamamos $U$ a $E_3 E_2 E_1 A$. Hasta acá podríamos decir que resolvimos el ejercicio haciendo sustitución hacia atrás, pero el problema es que por cada $b$ hay que recalcular el lado derecho de la ecuación.
Partiendo de

(20)
\begin{equation} E_3 E_2 E_1 A = U \end{equation}

multiplicando por las matrices inversas en orden inverso obtenemos

(21)
\begin{equation} E_1' E_2' E_3' E_3 E_2 E_1 A = E_1' E_2' E_3' U \end{equation}

la matriz resultante $E_1' E_2' E_3'$ es una matriz inferior, a la cual llamamos $L$. Hay que notar que $E_i' E_i$ es la matriz identidad, entonces

(22)
\begin{equation} A = E_1' E_2' E_3' U \end{equation}

llamando $L = E_1' E_2' E_3'$ logramos factorizar A en dos matrices, una inferior $L$ y otra superior $U$

(23)
\begin{equation} A = LU \end{equation}

para resolver un problema ahora remplazamos el nuevo a en la ecuación lineal 1.

(24)
\begin{equation} LU x = b \end{equation}

llamando $z = Ux$ podemos reescribirla asi

(25)
\begin{equation} Lz = b \end{equation}

como $L$ es inferior el problema sale con sustitución forward y una vez obtenido $z$ calculamos $x$ como

(26)
\begin{equation} Ux = z \end{equation}

como $U$ es superior se resuelve con sustitución backward.

Factorización de Cholesky

Para matrices definidas positivas, existe una matriz $L$ tal que $A=LL^T$. Los elementos de $L$ se calculan así:

  • La diagonal: $a_ii=\sum_{k=1}^i l_{ik}^2$
  • Lo que no es diagonal: $a_ij=\sum{k=1}^j l_{ik}l_{ij}$

y después de ahí hay que despejar los $l_{ij}$.

Links Interesantes

http://www.unalmed.edu.co/~ifasmar/cap3.pdf

Unless otherwise stated, the content of this page is licensed under Creative Commons Attribution-ShareAlike 3.0 License