5: Álgebra lineal
( \newcommand{\kernel}{\mathrm{null}\,}\)
Consideremos el sistema de ecuacionesn lineales en incógnitas, dado por
a11x1+a12x2+⋯+a1nxn=b1,a21x1+a22x2+⋯+a2nxn=b2,⋮⋮an1x1+an2x2+⋯+annxn=bn.
Podemos escribir este sistema como la ecuación matricial
Ax=b,
con
A=(a11a12⋯a1na21a22⋯a2n⋮⋮⋱⋮an1an2⋯ann),x=(x1x2⋮xn),b=(b1b2⋮bn)
En este capítulo se detalla la solución numérica de (5.1).
Eliminación Gaussiana
El algoritmo numérico estándar para resolver un sistema de ecuaciones lineales se llama Eliminación Gaussiana. Podemos ilustrar este algoritmo con el ejemplo.
Considerar el sistema de ecuaciones
−3x1+2x2−x3=−1,6x1−6x2+7x3=−7,3x1−4x2+4x3=−6.
Para realizar la eliminación gaussiana, formamos una Matriz Aumentada combinando la matriz A con el vector de columna→b:
(−32−1−16−67−73−44−6).
Luego se realiza la reducción de filas en esta matriz. Las operaciones permitidas son (1) multiplicar cualquier fila por una constante, (2) agregar múltiplo de una fila a otra fila, (3) intercambiar el orden de cualquier fila. El objetivo es convertir la matriz original en una matriz triangular superior. Comenzamos con la primera fila de la matriz y trabajamos nuestro camino hacia abajo de la siguiente manera. Primero multiplicamos la primera fila por 2 y la agregamos a la segunda fila, y agregamos la primera fila a la tercera fila:
(−32−1−10−25−90−23−7).
Después vamos a la segunda fila. Multiplicamos esta fila por−1 y la agregamos a la tercera fila:
(−32−1−10−25−900−22).
Las ecuaciones resultantes se pueden determinar a partir de la matriz y están dadas por
−3x1+2x2−x3=−1−2x2+5x3=−9−2x3=2.
Estas ecuaciones se pueden resolver mediante sustitución hacia atrás, partiendo de la última ecuación y trabajando hacia atrás. Tenemos
−2x3=2→x3=−1−2x2=−9−5x3=−4→x2=2−3x1=−1−2x2+x3=−6→x1=2
Por lo tanto,
(x1x2x3)=(22−1)
Descomposición de LU
El proceso de Eliminación Gaussiana también da como resultado la factorización de la matriz A para
A=LU,
dondeL es una matriz triangular inferior yU es una matriz triangular superior. Usando la misma matriz A que en la última sección, mostramos cómo se realiza esta factorización. Tenemos
(−32−16−673−44)→(−32−10−250−23)=M1 A
donde
M1 A=(100210101)(−32−16−673−44)=(−32−10−250−23).
Tenga en cuenta que la matrizM1 realiza la eliminación de filas en la primera columna. Dos veces la primera fila se agrega a la segunda fila y una veces la primera fila se agrega a la tercera fila. Las entradas de la columna deM1 provienen de2=−(6/−3) y1=−(3/−3) según se requiera para la eliminación de fila. El número−3 se llama el pivote.
El siguiente paso es
(−32−10−250−23)→(−32−10−2500−2)=M2(M1 A),
donde
M2(M1 A)=(1000100−11)(−32−10−250−23)=(−32−10−2500−2).
Aquí,M2 multiplica la segunda fila por−1=−(−2/−2) y la agrega a la tercera fila. El pivote es−2.
Ahora tenemos
M2M1 A=U
o
A=M−11M−12U.
Las matrices inversas son fáciles de encontrar. La matrizM1 multiplica la primera fila por 2 y la agrega a la segunda fila, y multiplica la primera fila por 1 y la agrega a la tercera fila. Para invertir estas operaciones, necesitamos multiplicar la primera fila por−2 y agregarla a la segunda fila, y multiplicar la primera fila por−1 y agregarla a la tercera fila. Para verificar, con
M1M−11=I,
tenemos
(100210101)(100−210−101)=(100010001).
Del mismo modo,
M−12=(100010011)
Por lo tanto,
L=M−11M−12
está dado por
L=(100−210−101)(100010011)=(100−210−111)
que es triangular inferior. Los elementos fuera de la diagonal deM−11 y simplementeM−12 se combinan para formar L. Nuestra descomposición de LU es por lo tanto
(−32−16−673−44)=(100−210−111)(−32−10−2500−2).
CAPÍTULO 5. ALGEBRA LINEAL Otra característica agradable de la descomposición de LU es que se puede hacer sobrescribiendo A, por lo tanto, ahorrando memoria si la matrizA es muy grande.
La descomposición de LU es útil cuando uno necesita resolverAx=b parax cuandoA es fijo y hay muchos diferentesb, primero determinaL yU usa la eliminación gaussiana. Entonces uno escribe
(LU)x=L(Ux)=b.
Dejamos
y=Ux,
y primero resolver
Ly =b
paray por sustitución directa. Luego resolvemos
Ux=y
parax por sustitución hacia atrás. Si contamos operaciones, podemos demostrar que resolver (LU)x=b es sustancialmente más rápido una vezL yU están en la mano que resolviendoAx=b directamente por eliminación gaussiana.
Ahora ilustramos la solución deLUx=b usar nuestro ejemplo anterior, donde
L=(100−210−111),U=(−32−10−2500−2),b=(−1−7−6)
Cony=Ux, primero resolvemosLy=b, es decir
(100−210−111)(y1y2y3)=(−1−7−6).
Uso de sustitución directa
y1=−1,y2=−7+2y1=−9,y3=−6+y1−y2=2.
Ahora resolvemosUx=y, es decir
(−32−10−2500−2)(x1x2x3)=(−1−92).
Usando la sustitución hacia atrás,
−2x3=2→x3=−1,−2x2=−9−5x3=−4→x2=2,−3x1=−1−2x2+x3=−6→x1=2,
y una vez más hemos determinado
(x1x2x3)=(22−1).
Pivote parcial
Al realizar la eliminación gaussiana, el elemento diagonal que se utiliza durante el procedimiento de eliminación se denomina pivote. Para obtener el múltiplo correcto, se utiliza el pivote como divisor a los elementos debajo del pivote. La eliminación gaussiana en esta forma fallará si el pivote es cero. En esta situación, se debe realizar un intercambio de filas.
Incluso si el pivote no es idéntico a cero, un valor pequeño puede resultar en grandes errores de redondeo. Para matrices muy grandes, uno puede perder fácilmente toda la precisión en la solución. Para evitar estos errores de redondeo derivados de pequeños pivotes, se realizan intercambios de filas, y esta técnica se denomina pivotamiento parcial (el pivotamiento parcial contrasta con el pivotamiento completo, donde se intercambian tanto filas como columnas). Ilustraremos con el ejemplo la descomposición de LU mediante pivotamiento parcial.
Considerar
A=(−22−16−673−84).
Intercambiamos filas para colocar el elemento más grande (en valor absoluto) en el pivote, oa11, posición. Es decir,
A→(6−67−22−13−84)=P12 A
donde
P12=(010100001)
es una matriz de permutación que cuando se multiplica a la izquierda intercambia la primera y segunda filas de una matriz. Tenga en cuenta queP−112=P12. El paso de eliminación es entonces
P12 A→(6−67004/30−51/2)=M1P12 A
donde
M1=(1001/310−1/201).
El paso final requiere un intercambio de filas más:
M1P12 A→(6−670−51/2004/3)=P23M1P12 A=U.
Dado que las matrices de permutación dadas porP son sus propias inversas, podemos escribir nuestro resultado como
(P23M1P23)P23P12 A=U.
MultiplicaciónM de a la izquierda porP intercambios de filas mientras que la multiplicación a la derecha por columnas deP intercambios. Es decir,
P23(1001/310−1/201)P23=(100−1/2011/310)P23=(100−1/2101/301).
El resultado neto enM1 es un intercambio de los elementos no diagonales1/3 y−1/2.
Entonces podemos multiplicar por la inversa de(P23M1P23) para obtener
P23P12 A=(P23M1P23)−1U,
que escribimos como
PA = LU.
En lugar de L, MATLAB escribirá esto como
A=(P−1 L)U.
Por conveniencia, solo denotaremos(P−1 L) porL, pero porL aquí nos referimos a una matriz triangular inferior permutada.
Programación MATLAB
En MATLAB, para resolverAx=b parax usar la eliminación gaussiana, un tipo
donde∖ resuelve elx uso del algoritmo más eficiente disponible, dependiendo de la forma de A. Si A es unan×n matriz general, entonces primero se encuentra la descomposición de LU de A usando pivotamiento parcial, y luegox se determina a partir de sustitución permutada hacia adelante y hacia atrás. SiA es triangular superior o inferior, entonces se usa directamente la sustitución hacia adelante o hacia atrás (o su versión permutada).
Si hay muchos lados diferentes de la mano derecha, uno primero encontraría directamente la descomposición LU de A usando una llamada a una función, y luego resolvería usando∖. Es decir, uno iteraría para diferentesb las siguientes expresiones:
donde la segunda y tercera líneas se pueden acortar a
donde se requiere el paréntesis. En conferencia, voy a demostrar estas soluciones en MATLAB usando la matrizA=[−3,2,−1;6,−6,7;3,−4,4] y el lado derechob=[−1;−7;−6], que es el ejemplo que hicimos a mano. Aunque no detallamos el algoritmo aquí, MATLAB también puede resolver el problema del valor propio del álgebra lineal. Aquí, el problema matemático a resolver viene dado por
Ax=λx,
dondeA es una matriz cuadrada,λ son los valores propios, y losx′ s asociados son los vectores propios. La subrutina MATLAB que resuelve este problema es eig.m. Para encontrar solo los valores propios deA, un tipo
Para encontrar tanto los valores propios como los vectores propios, un tipo
Puede encontrar más información en las páginas de ayuda de MATLAB. Una de las buenas características de la programación en MATLAB es que no se comete un gran pecado si uno usa una función incorporada sin gastar el tiempo requerido para comprender completamente el algoritmo subyacente.