Eliminación gaussiana para sistemas

de ecuaciones lineales
Fecha de recepción: Abil, 1996
E du cación Matemáti ca
Vol. JO No . I Abril 1998
Pedro Ramón Almeida Benítez y José Ramón Franco Brañas
pp. 74-88
Uni versidad de La Laguna, Dcpto de Análisis Matemático
38721. La Laguna, Tenerife, España
jfranco@ull .es

Resumen. La solución de sistemas de ecuaciones lineales es el problema central del

álgebra lineal. En efecto, muchos problemas de ingeniería precisan modelos matemáticos
cuya solución requiere métodos de algebra lineal. La fórmula llamada +regla de Cramer-V
implica determinantes y es muy laboriosa para resolver sistemas con más de cuatro o
cinco variables. En este artículo, se considera la eliminación Gaussiana y es evaluado
su coste computacional.
Las técnicas de computación para resolver grandes sistemas de ecuaciones lineales
descansan sobre métodos matriciales y la fáctdr-ización L U es uno de los más populares.
La factorización L U y sus variantes Doolittle, Crout y Cholesky son estudiadas en las
siguientes líneas.

Abstract. The solution of linear systems of equations is the central problem of linear
algebra. Indeed, many engineering lead to mathematical models whose solution requires
methods of linear algebra. The formula called Cramer's rule involving determinants is
very laborious to so/ve systems with more than four or five variables. In this paper, the
Gaussian elimination is considered and its computational costs are evaluated.
Computing techniques for solving large sets of linear equations rely on matrix methods
and LU descoposition and its variants Doolittlee, Crout and Cholesky are studied in the
next fines.

Muchos problemas de Ingeniería o Matemática Aplicada se reducen a resolver un sistema

de ecuaciones lineales. De aquí la importancia de estos y el interés de sus métodos de
resolución. La mayoría de los métodos -:-directod de resolución tienen su origen en el
algoritmo de Gauss y el objeto de estas líneas es introducir el estudio de los más
important-es: Doolittle, Crout y Cholesky.

Introducción.

Deseamos resolver el sistema den ecuaciones lineales y n incógnitas:

E¡: a¡¡X¡ + a¡2X2 + ... + ªinXn = b¡ }

E2: ª21X1 + a22x2 + ... + ª2rrn = b2

En: an¡X¡ + ª n2X2 + ... + ªnrrn = bn

e 74
1'11 E LIMINACIÓN GAUSSIANA PARA SISTEMAS DE ECUACIONES LINEALES B Pág . 75 R

que podemos escribir matricialmente

Ax b
siendo

A=
[
~ 11 .. ·
.
a 1,! ]

ª111 •. . ª1111
. ; X
_Ü
r
-~l

- b [; ,
b¡
l
donde suponemos que det (A) -:t O.

Ineficacia de la regla Cramer para resolver grandes sistemas.

Vamos a intentar resolver dicho sistema empleando el método conocido con el nom-
bre de regla de Cramer, calculando la incógnita x¡ en la forma:

a¡¡ · ·· ªt,i - l b¡ ªt ,i+l · · · ªin

X¡ =
det (A)
¿Cuál será el número de operaciones a realizar empleando éste método?. Será
preciso calcular n + 1 determinantes de orden n (*). En total:
(n + 1) · n! · (n - 1) multiplicaciones
n divisiones
Sumando:
(n + 1) · n! · (n - 1) + n = (n 2 - 1) · n! + n
Hay que advertir que hemos considerado únicamente multiplicaciones y
divisiones. Esto no es demasiado riguroso, pero dado que un ordenador suele tardar
más tiempo en realizar una multiplicación que una suma, no tomamos en consideración
éstas últimas a la hora de evaluar el número de operaciones.
Por otra parte, al resolver un problema técnico, suelen aparecer sistemas
con un elevado número de ecuaciones, varios cientos o miles de ellas en muchas ocasio-
nes. La tabla siguiente indica el número de operaciones necesarias para distintos va-
lores den:
n operac10nes
3 51
5 2885
10 358,657,210

En total, 358,657,210 operaciones, paran = 10. Si consideramos que un

ordenador personal efectúa unas 106 operaciones por segundo, para resolver el sistema

(*) Recordemos que se llama determinante de la matriz cuadrada A = (a;j) E Mnxn a la suma de los n! productos de la forma
(-l)cr(i1, i2, .. .. .in) ªu,. ª2;2 ... ª•in
en donde cada uno de éstos elementos, y sólo uno, procede de una fila cualquiera y uno, y sólo uno, de una columna cualquiera.
 Pág . 76 • E oucACióN MATEMÁTICA • Vo l. 1 O No. l · Abril 1998 • © GEi •

anterior ne cesitaría alrededo r de 6 minutos de trab ajo . Estos cálculos son muy
aproximados puesto que la velocidad de los ordenadores varía día a día y depende mucho
del tip o de ordenador.
Por otra parte, si para resolver un sistema tan pequeño (n = 1O) se nece-
sita tanto tiempo, dicho méto do es abso lutamente inserv ible para sistemas con un
" elevado" número de ecuaciones.

Nuevo intento de resolución mediante la matriz inversa.

Descartada la regla de Cramer, vamos a intentar resolverlo de la siguiente forma:

AX= b => A- 1 · A · X= A- 1 · b => X= A- 1 · b

en donde A- 1 es la inversa de la matriz A que, como sabemos, es igual a A- 1

(A*Y
det (A)'
siendo A * la matriz adjunta de la matriz A .
En este caso, el coste computacional sería:
Cálculo de A*
2
(n determinantes de orden n - 1) Cálculo det (A) Cálculo X
Mult.: n2 · (n - 1)! · (n - 2) + n! · (n - 1) + n2 =
= (n - n - 1) · n! + n 2
2
Di~: n2
En total : (n 2 - n - l ) · n! + 2n 2 .
Si hacemos una tabla para distintos valores den:

n operaciones
3 48
5 2330
10 322,963,400

resultando, para n = 1O, unos 35 millones de operaciones menos que en el método

de Cram er. Haciendo las mismas consideraciones que entonces, vemos que serían
necesarios más de 5 minutos para su resolución. Llegamos a la conclusión de que
éste método también es inaplicable a la resolución de "grandes" sistemas.
Pero, por si esto fuera poco, el mayor inconveniente de un método con
tantas operaciones lo constituye la acumulación de errores de redondeo.

El método de Gauss.

Otra forma de poder resolver el sistema podría ser la conocida con el nombre de mé-
todo de Gauss o eliminación gaussiana, que consiste ~n transformar el sistema dado
A .X= b 1
en otro equivalente

U .x = e
en el que la matriz U es triangular superior.
• E LIMINACIÓN GAUSSIANA PARA SISTEMAS DE ECUACIONES LINEALE-5
• Pág . 77 •

Tomaremos un caso práctico sencillo y desarrollaremos dicho método

detalladamente, dado su interés posterior.
Sea el sistema:

2x \ ; : ; =-= _ 1)
27
- 2, + 2y + z
Vamos a dar los siguientes pasos:
(1) Sustraemos de la segunda ecuación, la primera multiplicada por 2.
A este multiplicador lo llamaremos 121 = 2 y ha sido obtenido en la forma

l
a21la¡¡ = 4/2 = 2.
Tenemos el sistema equivalente resultante:
2x +y+ z = 1
-y - 2z = -4
-2x + 2y + z = 7
(2) Sustraemos de la tercera ecuación, la primera multiplicada por - 1.

l
El multiplicador 1 31 = a 31 !a 11 = -2/ 2 = -1. Obtenemos el sistema
equivalente:
2x + y + z = 1
-y - 2z = -4
3y + 2z = 8
(3) Sustraemos de la tercera ecuación, la segunda multiplicada por - 3, siendo
132 = a 3 / 2)!a2/ 2) = 3/-1 = - 3 (*), y resultando el sistema triangular equivalente:
2x+y+z = I)
-y - 2z = -4 ~ UX = c
-4z = -4
Los números 121 = 2, 131 = -1 y 132 = -3, son los llamados multiplica-
dores de la eliminación gaussiana. En general, 1ij representa la cantidad por la que
multiplicamos la ecuación} para sustraerla de la ecuación i y así producir un cero en la
entrada (i,j) .
A los números a 11 = a 1/ 2) = 2, a 2/ 2) = - 1 y a 33 (3) = -4, que se utilizan
para eliminar entradas que están en su misma columna, los llamaremos pivotes.
Hemos seguido los pasos anteriores de un modo sistematizado, con el fin
de que conduzcan a un algoritmo fácilmente programable.
Con esto hemos terminado la primera fase del proceso de eliminación gau-
ssiana, la llamada eliminación hacia adelante, que nos ha permitido obtener un sis-
tema reducido equivalente al inicial. Para culminar la resolución del sistema·nos queda
una segunda fase, llamada sustitución hacia atrás o regresiva, que consiste en hallar
el valor de la incógnita de la última ecuación, sustituirla en la ecuación anterior calculan-
do el valor de la otra incógnita, y así sucesivamente.
En el sistema de nuestro ejemplo, de la tercera ecuación resulta z = I y, lle-
vado éste valor a la segunda obtenemos y = 2. Finalmente, de la primera se obtiene
X = - ].
Vamos ahora a calcular el coste computacional de éste método:
a) Para la eliminación hacia adelante, esto es, para transformar el sistema
A· x = b en el U· x = c, tenemos que efectuar:

(*) Hemos llamado a 2Pl y a 32(2l a los coeficientes resultantes del paso (2), para distinguirlos de los a22 y a32 originales.

: 1
1
 Pág . 7 8 • E DUCACIÓN MATEMÁTICA • Vol. l O No. l · Ab ril 1998 • © GEi B

Transformac ión de A Segundos miembros

(multiplicamos n elementos
de la primera fila por 11 - 1
multiplicadores para anular
la primera columna)
Productos: n (n - 1) + (17 - 1) (n - 2) + ... + 3 · 2 + 2 · 1 + (n - 1) + (17 - 2) + 2 + 1 =

=
11

L, i (i - 1) +
n-1

1> = I .,
11

l -
11

L, i +
11 - I
L, i 17 (17 + 1)( 217 + 1)
i == 1 i= I 1=2 i=1 i= I 6
- 1 - 17 + 1 = 17 3/3 + n 2!2 Sn/6
11 - I
17(11 - 1)
Di visio nes: (11 - 1) + (n - 2) + ... + 2 + 1 L, i necesarias
i=I 2
para hallar los multiplicadores.
b) Para la sustitución hacia atrás tenemos que efectuar:

Producto s : (17 - 1) + (17 - 2) + ... + 2 + 1 n(l7 l) al sustituir

2
los valores de las incógnitas en las ecuaciones anteriores .

Divis iones: n.

En tota l, y teniendo en cuenta que las magnitudes de 17 y 17 2 son despreciables frente a

las de n3, puesto que nuestro objetivo es resolver sistemas para valores de 17 "grandes":
17 2 517 n(17 - l ) 173
+ 172 - -417
n3
Productos: + - + = 0(17 3/3)
3 2 6 2 3 3
17(n -1) -17
n2
Divisiones: + n = + 0(17 2/2)
2 2 2
173 3n 2 5n
O(n 3!3)
y sumando: - + -- - =
3 2 6
que es considerablemente inferior al coste computacional de los dos métodos anterio-
res y, para un sistema de 1Oecuaciones con 1Oincógnitas (n = 1O), el número de operacio-
nes sería de 475, que un ordenador personal tarda menos de un segundo en resolver(*) .
Podemos resumir en la siguiente tabla el cómputo del número de operaciones
al incrementar 17:
n Operaciones
3 20
10 475
50 45 ,377
100 348,250

(*) Sean/(11) y g(11) funcio nes den. Usamos la notación :

j(n) = O(g(n))

si para una constante K y n sufic ientemente grande l/(11)/g(11)I $ K. Decimos ento nces que f(n) tie ne el orden de magnitud
de g(n) co mo mucho .
Por ejemplo, si /{n) = 11 3/3 + n2rl - Sn/6 , escribimos/{n) = O (11 3/3), ya que para valores de 11 "grandes", la contribu-
ción de n 2/2 y 511/6 es despreciable frente a la de n3/3.
 E U MINACIÓN GAUSS IANA PARA SISTEMAS DE EC UACI ONES LI NEA LES • Pág . 79 a

Una variante de éste procedimiento es el llamado méto do de Gauss-Jordan ,

que consiste en eliminar no sólo X¡ de las ecuaciones E ¡ + 1, E ¡ + 2 , .. . , En, como hici-
mos en el método de Gauss, sino también de E2, ... , E ¡ _ 1, reduciendo el sistema
original a un sistema diagonal. Es decir, eliminamos en cada columna los elementos
que están situados por debajo y por encima de la diagonal principal, hallando los
multiplicadores correspondientes. Este procedimiento evita la necesidad de la sustitu-
ción hacia atrás en la eliminación gaussiana. Sin embargo, el número de operaciones
es algo mayor que en el método de Gauss y por ello poco utilizado .
En efecto, para hacer nulas las entradas por encima de la diagonal principal,
habrá que efectuar:
Productos: (n -1) · 1 + (n -2) · 2 + (n - 3) · 3 + ... + (n - (n - 1)) · (n - 1)

= ~ (n - i). i = ~ n . i - 'f ¡2 = n. ~ . n - n(n + 1)(2n + 1) =

i= I i=l i= l 2 6

Divisiones: n( n - l) = O(n 2/2), como vimos anteriormente.

2
Por tanto, el coste adicional del método de Gauss-Jordan será de O(n 3/6).
En total: O(n 3/3) + O(n3/6) = O(n 3 /2) .
Una variante de los dos métodos anteriores, Gauss y Gauss-Jordan , podría ser
el siguiente método híbrido: aplicamos la eliminación gaussiana al sistema y lo reducimos
a un sistema triangular superior. Después utilizamos la n-sima ecuación para hacer
nulos los coeficientes de xn en las primeras n - 1 ecuaciones. El coste computacional
de éste método es análogo al del método de Gauss: O(n 3!3).

Eliminación de Gauss con pivoteo parcial. ·

Si en el transcurso del anterior logaritmo resultara algún pivote nulo, no podríamos

continuar con su proceso. Para evitar esto, si es posible, utilizaremos como pivote en 1
cada columna otra entrada no nula en lugar de la entrada diagonal correspondiente.
Esto significa un intercambio de las filas de la matriz ampliada que contienen dicha
1
entrada diagonal nula y el pivote seleccionado. Este intercambio siempre será posible
si la matriz es regular. ·!
Es la práctica, hay otra razón que recomienda realizar un intercambio de
filas. Cuando tomamos como pivote un número muy pequeño, se produce inestabi-
1
lidad en la eliminación gaussiana, produciéndose un sustancial error de redondeo. Para
ilustrar esto, mostramos la resolución del siguiente sistema:
0.003000x 1 + 59.14x2 = 59.17
5.291 X¡ - 6.130 X2 = 46.78
l ¡
1

cuya solución exacta es x 1 = 10.00; x 2 = 1.000. 11

.1
Si usamos aritmética de cuatro dígitos con redondeo, tenemos:
a 11 = 0.003000; 121 = 5.291 /0.003000 = 1763.66 jl
que nos da el sistema escalonado ¡I
1
0.003000; X¡ + 59.14 X2 = 59.17
\!
-104300 X2 = -104400 11
¡l
1

¡
1
1 1
 Pág. 80 a E DUCACIÓN MATEMÁTICA a Vol. 1O No. 1 · Abril 1998 • © GEi a

que, al resolverlo por sustitución hacia atrás, se obtiene

X¡ = -10.00; Xz = 1.001

produciendo en x 1 un error del 200% .

Para evitar esto, tomaremos como pivote en cada columna la entrada de
mayor valor absoluto por debajo de la correspondiente entrada diagonal.
En definitiva, ya sea por el motivo teórico de una entrada diagonal nula o por
el práctico de una entrada casi nula, la eliminación gaussiana con intercambios consiste
en elegir como pivote, en cada momento, la entrada de cada columna que está bajo la
diagonal principal y que tiene el máximo valor absoluto; i.e., se determina p tal que:

Esta técnica se conoce como pivoteo máximo de columna o pivoteo parcial.

Escalamiento.

Otra técnica es la conocida con el nombre de pivoteo escalado de columna, que consiste
en definir para cada ecuación un factor de escala s ¡ = máx la¡¡I.A continuación, elegimos
ISJSn
el entero k que verifique:

y realizamos el intercambio de las ecuaciones E 1 y Ek.

La acción de escalar tiene por objeto asegurar que el elemento mayor de cada
ecuación tenga una magnitud relativa igual a 1 antes de empezar la comparación para
el intercambio de filas. El escalamiento se hace solamente con propósitos de compara-
ción y, por tanto, la división entre los factores de escala no produce un error de re-
dondeo en el sistema.
Veamos un ejemplo. Sea el sistema:
E 1: 30.00 x 1 + 591400 x 2 = 591700
E 2 : 5.291 x 1 - 6.130x2 = 46.78
que es el mismo de un ejemplo anterior salvo que hemos multiplicado los coeficien-
tes por 104 . Aplicando el método de Gauss obtenemos los mismos falsos resultados
anteriores x 1 = -10.00 y x 2 = 1.001, debido a la diferencia relativa entre los valo-
res de los coeficientes.
Tenemos que s 1 = 591400, s 2 = 6130. Por tanto, la11 1fs 1 = 30.00/591400 =
0.5073 x 10-4 y la 21 1fs 2 = 5.291/6.130 = 0.8631, por lo que se hace el intercam-
bio de E 1 y E 2 .
Aplicando ahora la eliminación gaussiana con pivoteo parcial, obtenemos
los resultados exactos x 1 = 10.00 y x 2 = 1.000.

Eliminación de Gauss con pivoteo total.

Otro modo de proceder es el llamado pivoteo total, que consiste en realizar los cambios
necesarios entre filas y/o columnas para situar en la posición del pivote el mayor ele-
mento en valor absoluto de la submatriz A/r) en cada paso r (ver figura):
• ELIMI NACIÓN GAUSSIANA PARA SISTEMAS DE ECUACIONES LINEALES • Pág. 81 •

e o l umna r

f- flla r
1
1

I¡

Referente al coste computacional de la eliminación de Gauss con intercam- 1

bios hemos· de decir que, dado que sólo consideramos las operaciones elementales
aritméticas (que son las que producen errores de redondeo), nos va a dar el mismo coste 1
que en el método de Gauss simple, si bien el tiempo de cálculo será algo mayor en la
práctica debido a las operaciones necesarias para el intercambio de filas o columnas.
Por tanto, el tiempo empleado en resolver el sistema será algo superior al del
método de Gauss simple y también ventajoso.

Matriz definida positiva.

n n
Una matriz cuadrada AEMnxn es definida positiva si v'Av = L ¿ v¡. ªu . vj > O,
i=l j=l

\:/vE R 11 - {O} .
Teorema.- La matriz AE Mnxn es definida positiva si y sólo si:

det=
~ll
:
...ª1.k]
: > O, k = 1, ... , n.
rªkt ... ªkk

Ver demostración en Strang (6].

Una matriz cuadrada AEMnxn nxn es estrictamente dominante diagonal-
n

mente si la¡¡¡ > ¿

j =I
laul,J * i, i = I, 2, 3, ... , n.
Teorema.- Si la matriz A es definida positiva o estrictamente dominante
diagonalmente, podemos aplicar el método de Gauss al sistema Ax = b y obtener
su solución única sin intercambio de filas o columnas y dicha solución será estable
respecto a posibles errores de redondeo.
Puede verse la demostración de Wendroff(7).

Factorización de Doolittle.
Volvamos al método de Gauss simple y consideremos de nuevo el sistema:

2x+y+z=l)
4x +y= -2
-2x + 2y + z = 7
a Pág. 82 a ED UCACIÓN MATEMÁTICA a Vol. lONo. l ·Abril 1998 a © GEi a

que podemos escribir en forma matricial Ax = b:

Aplicarnos los pasos (1 ), (2) y (3), en éste orden, de ia eliminación hacia adelante
y llegamos al sistema escalonado Ux = c siguiente:

-11 21] [x1 [- 411

[2O00-4 z -4
y e

Buscaremos ahora la relación entre Ay U, para lo cual analizaremos cada uno

de los pasos anteriores (1), (2) y (3).
La transformación elemental ( 1) (*) podernos representarla con la matriz
elemental E21 siguiente:

E21 [-; ~ g]
OO 1
Ax= b
E21AX = E2 1b
Resultando:

La transformación (2) podernos representarla con la matriz elemental E 31 :

E3¡ [1
= O o1 O
1 O 1
º] Ax
E 21 Ax
= b
= E 21 b
E 31 E 21 Ax E 31 E 21 b

[H-ífü] Hl
Asimismo, la transformación (3) podemos representarla con la matriz elemen-

Ax b
E 21 Ax E 21 b
E 31 E 21 Ax E 31 E 21 b
E32E3¡E21 Ax E32E3¡E21b
(*) Las transformaciones elementales permitidas en un sistema, que conducen a un sistema equivalente, son tres:
1) Una ecuación puede multiplicarse por una constante k, distinta de cero, y sustituir en el sistema la anterior ecuación por
la así obtenida.
2) La ec.u ación Ej puede multiplicarse por una constante k, sumarla a la ecuación E¡, y sumar la ecuación resultante en lugar
de la E¡.
3) Dos ecuaciones cualesquiera se pueden intercamb iar.
11 ELIMINACIÓN GAUSSIANA PARA SIST EMAS DE ECUACIONES LI NEALES • Pág . 83 •

Tenemos , pues :
a) E 32E 31 E 21 A == U ; E 32E 3 1E 21 b == e
b) Para sustraer de la ecuación i la ecuación j multiplicada por liJ, se constituye
la matriz elemental E¡j cuya entrada no diagonal no nula es la que ocupa la posición (i,j)
y vale -liJ.
De ésta manera, la matriz que convierte A y ben U y e, respectivamente, es:

Aún tiene mayor interés determinar la matriz E-1, que deshace los pasos de
eliminación gaussiana, i . e., regresa de U a A y de e a b. Para ello comenzamos por
invertir el paso (3), que fue el último en realizarse; el segundo paso por invertir es el
(2) y el último es el paso (1). Esto está en consonancia con la regla para calcular la
inversa del producto, es decir, las inversas están en orden opuesto a la sucesión ori-
ginal de operaciones:
t;-l == (E32E3¡E21t 1 == E21- 1E3¡-IE32-I

Evaluemos cada una de las Ei/- 1:

El inverso del paso (3) es el (4): Sumamos a la tercera ecuación la segunda
multiplicada por-3 . La matriz que representa esta operación es:

El inverso del paso (2) es el (5): Sumamos a la tercera ecuación la primera

multiplicada por -1, que viene representada por la matriz:

Finalmente, el inverso del paso (1) es el (6): Sumamos a la segunda ecuación

la primera multiplicada por 2. Siendo la matriz que representa ésta operación:

Por último, de E 32E 31 E 21 A == U, Deducimos:

A == E21-IE3¡-IE32-IU
 Pág. 84 8 EDUCACIÓN MATEMÁTICA a Vol. 1O No. 1 • Abril 1998 a © G El a

Y haciendo L == E 21-lE31- 1E 32- 1, tenemos la factorización:

A == LU
Esta matriz L es la que deshace los pasos de la eliminación gaussiana, i. e.,
regresa de U a A y de e a b.
En nuestro ejemplo, L es:

L == [~ ~
O O 1
~i [~ ~ ~] [~~ ~]
-1 O 1 O -3 1
= [ ~~
-1 -3
~]
1

Por otra parte, la matriz L tiene las siguientes características:

a) Es subtriangular.
b) Las entradas diagonales son todas iguales a 1.
c) Las entradas subdiagonales son respectivamente iguales a los multi-
plicadores de la eliminación gaussiana.
Por otra parte, la matriz U de la factorizaciónA = LU es la matriz del sistema
escalonado, al que llegaríamos después de aplicar la eliminación gaussiana, y sus en-
tradas diagonales son los pivotes.
Esta factorización se conoce con el nombre de factorización de Doolittle.

El coste computacional para resolver un sistema por el método de Doolittle

es análogo al del método de Gauss. El cálculo de la matriz U (obtenida tras las etapas
del método de Gauss) exige realizar O(n 3/3) operaciones, como vimos anteriormente.
El cálculo de L no requiere operaciones adicionales ya que se obtiene sin
más que conservar los valores de los multiplicadores de la eliminación.
Por último, la resolución de los dos sistemas triangulares requiere 2n 2
operaciones.

Notas

1) No toda matriz A admite la factorización de Doolittle.

Por ejemplo, consideremos la factorización de Doolittle para matrices 2 x 2:

La matriz:
b = o
A-[~~J-w~ a~:; ac + d == O

nos propicia un sistema incompatible, aunque la matriz A sea regular.

2) Una matriz no regular sí puede aceptar la factorización de Doolittle.
Por ejemplo:

A = [~ ~] = [~ ~][~ ~]
• E LIMINACIÓN GAUSSIANA PARA SISTEMAS DE ECUACIONES LINEALES • Pág. 85 •

De 1 y 2 se deduce que la condición de regularidad no es necesaria ni sufi- ¡¡

ciente para la fa ctorización de Doolittle .
3) Se puede demostrar que la condición suficiente para que una matriz cua-
drada A sea factorizable Doolittle es que:

ªtt ... ªtk]

det (Ak) = det : : -:t O, para k = l, .. . , n.
[
ªkl ... ªkk
4) Se puede demostrar que para toda matriz cuadrada A existe una matriz .1
de permutación P (*) tal que PA sí admite la factorización de Doolittle PA = LU
(Ver Conde y Winter (3)).

La matriz de la nota 1) [ ~ ~ Jno admite la factorización de Doolittle, sin

embargo PA = [~ ~ J · [~ ~ J = [ ~ ~ J sí la admite.
5) Si la matriz A es invertible y factorizable Doolittle, ésta factorización
es única.
'
En efecto, supongamos que existen dos factorizaciones distintas: A = ,¡
·,

L 1U 1 = L 2 U2 . Por ser A invertible, también lo son L 1, U 1, L 2 y U2 . Además,

L 2 · L 1 = U2 · u 1- 1, y como r 2- 1 · L 1 es subtriangular con diagonal unitaria y ,1
U2 · u 1- 1 es supertriangular, debe verificarse r 2- 1 · L 1 = U2 · u 1- 1 = l. Por lo
1¡11
tanto, L 1 = L2 y U 1 = U2 . ,¡
I¡L
Determinantes. I!

La factorización de Doolittle nos proporciona un método poco costoso para calcular

il
!'

eldeterminantedeunamatrizA,puestoquesiA = LU,entoncesdetA = (detL) · (det ll

U) = det U = u 11 u22 . . . unn, ya que el determinante de una matriz triangular es
:¡
h
1,
igual al producto de sus elementos diagonales. Por tanto, el determinante de A es igual ,¡
al producto de lo~ pivotes de la eliminación gauusiana .

Factorización LU.

Vamos a intentar generalizar la anterior factorización de una matriz cuadrada A en la

forma de una matriz triangular inferior L y de una triangular superior U.
Si nos imponen condiciones a las entradas diagonales de L, o bien A o bien
PA , para una cierta matriz de permutación P, se pueden expresar como el producto de
una matriz subtriangular por una matriz supertriangular. Por ejemplo, consideremos
la matriz A, cuadrada de orden n:

(* ) Una matriz de permutación es aquella que tiene un elemento cuyo valor es I en cada fila y en cada columna y todos
sus demás elementos son cero. Es decir, es aquella que se obtiene a partir de la matriz unidad intercambiando filas y co-
lumnas entre sí. Por ejemplo:

p ¡ o
= O O1
[O1 O
º]
1!
1t I!i:
1
 · 86
Pag. • EoucAcióN MATEMÁTICA • Vol. 1 O No. 1 • Abril 1998 • © GE.! •

Podemos usar los 9 elementos conocidos de A para determinar parcialmente

los 6 elementos desconocidos de L y el mismo número de los de U Si embargo, el
procedimiento no nos lleva a una solución única, se necesitan tres condiciones adi-
cionales para los elementos de L y de U Si requerimos que l 11 == 122 == 13 3 == 1, el
método se conoce con el nombre de método de Doolittle, ya estudiado . Si tomamos
el sistema del ejemplo anterior:

2x +y+ z== 11
4x + y == -2
-2x + 2y + z == 7
tenemos que:

Si identificamos coeficientes:

Ull == 2; U¡2 == 1; U13 == 1

121 U11 == 4 => 121 == 2
121 · u 12 + u 22 == 1 => u22 == -1
121 ' U13 + U23 == O => U23 == -2
131 · u 11 == -2 => 131 = -1
131 · tt12 + 132 ' U22 = 2 => -1 - l3z == 2 => 132 = -3
131 · U¡3 + 132 ' U23 + U33 = 1 => -1 + 6 + U33 = 1 => t133 = -4
Con lo que:

A == LU == [ ~ ~
-1 -3 1
~] [~ -~
O O -4
-1]
que es la factorización de Doolittle, que ya obtuvimos anteriormente.
Volviendo al principio, si requerimos que u 11 = u 22 = u 33 = 1, tenemos
la factorización conocida con el nombre de método de Crout, similar a la de Doolittle.
Procediendo de modo análogo, factorizamos por el método de Crout la anterior matriz:

Por tanto, toda matriz A, o bien PA, siendo Puna matriz de permutación, pode-
mos factorizarla con éste método de Crout.
 ,.
1

• ELIMI NACIÓN GAUSSIANA PARA SISTEMAS DE ECUACIONES LINEALES a Pág . 87 • 1,

'
Volviendo nuevamente al principio, si ahora imponemos las condiciones
111 = = u22 y 133 = u33, tenemos la llamada factorización de Cholesky.
u 11, 122
Esta factorizació n no siempre es posible. Para ello es preciso que la matriz
A sea definida positiva, para garantizar la existencia de las raíces cuadradas que apa-
recen en los cálculos.

Ejemplo. - Vamos a factorizar la matriz A = 2 S2 [4221

1 2 2

A= LU = [4
2 52221

Identificando coeficientes:
1 2 2
= [1~ o
/3
131 132 rn~ U¡2
/3
o
u,11
U23
y

A= LU [4 221
= 2 S 2 =
[ 1
2 O
2 O
O 1[2O I2 1/2
1 1
1 2 2 1/2 3/4 [Link].¡4 O O [Link].¡4
Otra consideración importante acerca de la factorizac ión de Cholesky la
constituye el hecho de que si la matriz A es simétrica, entonces la matriz U = Lt, y

factorizamos A en la Forma A = L · D . Si la matriz anterior A = [ª

ª21
11
ªª2212 J es
simétrica, entonces a 12 = a 21 y 121 = u 12, con lo que U = Lt. El coste computacional
del método de Cholesky para matrices simétricas es del orden O(n 3/6) (ver Burden [3]).
Entonces, si A es simétrica y definida positiva, Ax = b {=> (LV)x = b,
que resolveremos a través de los sistemas siguientes:
1) Resolución del sistema subtriangular Ly = b.
2) Hallado el vector y, resolución del sistema supertriangular V x = y.

Notas

1) Evidentemente, la matriz L no es igual a la matriz L del método de Doolittle.

2) La factorización LLt no es única. Combinando los signos de las raíces ob-
tenemos distintas matrices.

Ej .: Sea la matriz A = [ ~ ~ J. Podemos tomar a = f4 = -2 y obtenemos

12 1 = -1, u12 = -1 y '3 = ±2, resultando:

A = [4 2] = [-2 º] [-2 -lJ [2 º] [2 J

2 5 -1 2 O 2
=
1 2
1
O2

3) Si la matriz A no es definida positiva, las matrices L y U de la factoriza-

1
ción de Cholesky son complejas. ,¡

:1

!
 Pág. 88 • E DUCACI ÓN MATEMÁTI CA • Vol. 1 o No. l • Abril 1998 • © GEi •

4) Si la matriz es definida positiva y, por tanto, es aplicab le el método de

Ch olesky, los cálculos son estables respecto al crecimiento de los
errores de redondeo . (Wendroff [7]).

Referencias

[ l] ALLEN SMITH, W. , Añá!isis Numérico, [5] MICHAVILA, F. ET GAVETE , L., Pro-

Editorial Prentice-Hall, México, 1988. gramación y Cálculo Numérico, Edito-
(2] BURDEN, R. ET DOUGLAS , J., Aná- rial Reverté, Barcelona, 1985.
lisis Numérico, Grupo Editorial Ibero- (6] STRANG, G., Linear Algebra and its
america, México, 1985. Applications, Harcourt Brace Publishers,
[3] CONDE, C. ETWINTER, G., Métodos y San Diego, 1986.
Algoritmos Básicos del Algebra Numé- [7] WENDROFF, B., Theoretical Numeri-
rica, Editorial Reverté, Barcelona, 1990. ca!Analysis, Academic Press, NewYork,
[4] CONTE, S.D . ET DE BOOR, C. , Aná- 1966.
lisis Numérico, Editorial McGraw-
Hill, Madrid, 1974.

i
1
¡