Notas Análisis Numérico Grupo 26 Unidad 3 Dr.

Guillermo Alberto Sánchez Lozano

Solución numérica de ecuaciones

lineales
Normalmente cuando resolvemos problemas en ingeniería
es normal encontrar sistemas de ecuaciones simultáneas:
𝑓 (𝑥 , 𝑥 , … , 𝑥 ) = 0
𝑓 (𝑥 , 𝑥 , … , 𝑥 ) = 0
[1]
⋮
𝑓 (𝑥 , 𝑥 , … , 𝑥 ) = 0

En ingeniería y física es común lidiar con el caso particular

en que todas las 𝑓 mencionadas, son sistemas de ecuaciones
lineales:
𝑎 𝑥 + 𝑎 𝑥 + 𝑎 𝑥 + ⋯+ 𝑎 𝑥 − 𝑏 = 0
𝑎 𝑥 + 𝑎 𝑥 + 𝑎 𝑥 + ⋯+ 𝑎 𝑥 − 𝑏 = 0
[2]
⋮
𝑎 𝑥 + 𝑎 𝑥 + 𝑎 𝑥 + ⋯+ 𝑎 𝑥 − 𝑏 = 0

Dicho sistema de ecuaciones por lo general se reescribe de

la siguiente manera:
𝑎 𝑥 + 𝑎 𝑥 + 𝑎 𝑥 + ⋯+ 𝑎 𝑥 = 𝑏
𝑎 𝑥 + 𝑎 𝑥 + 𝑎 𝑥 + ⋯+ 𝑎 𝑥 = 𝑏
⋮
𝑎 𝑥 + 𝑎 𝑥 + 𝑎 𝑥 + ⋯+ 𝑎 𝑥 = 𝑏

A lo largo de esta unidad, nos centraremos en los métodos

de resolución de este tipo de ecuaciones, además claro, de
ciertos temas de interés que atañen a las matrices, tales
como eigenvalores y eigenvectores.

Solución de sistemas de ecuaciones lineales

de orden
Los sistemas de ecuaciones lineales de orden 𝑛 ≤ 3
generalmente son fáciles de resolver mediante el esquema
clásico de Cramer, el cual consiste esencialmente en resolver
Notas Análisis Numérico Grupo 26 Unidad 3 Dr. Guillermo Alberto Sánchez Lozano

el siguiente sistema de ecuaciones mediante el uso de

determinantes.
𝑎 𝑥+𝑎 𝑦+𝑎 𝑧 =𝑏 𝑎 𝑎 𝑎 𝑥 𝑏
𝑎 𝑥+𝑎 𝑦+𝑎 𝑧 =𝑏 ⇒ 𝑎 𝑎 𝑎 𝑦 = ⎜𝑏 ⎞
⎛ ⎟ [3]
𝑎 𝑥+𝑎 𝑦+𝑎 𝑧 =𝑏 𝑎 𝑎 𝑎 𝑧 ⎝𝑏 ⎠

Cuyas soluciones
𝑏 𝑎 𝑎 𝑎 𝑏 𝑎 𝑎 𝑎 𝑏
𝑏 𝑎 𝑎 𝑎 𝑏 𝑎 𝑎 𝑎 𝑏
𝑏 𝑎 𝑎 𝑎 𝑏 𝑎 𝑎 𝑎 𝑏
𝑥= 𝑎 𝑎 𝑎 ,𝑦 = 𝑎 𝑎 𝑎 ,𝑧 = 𝑎 𝑎 𝑎 [4]
𝑎 𝑎 𝑎 𝑎 𝑎 𝑎 𝑎 𝑎 𝑎
𝑎 𝑎 𝑎 𝑎 𝑎 𝑎 𝑎 𝑎 𝑎

Nótese que para resolver un sistema de ecuaciones lineales

de 3 × 3 se requieren calcular 4 determinantes. El número
de multiplicaciones necesarias para calcular una sola
determinante de una matriz 𝐴 ∈ 𝑀 × se puede calcular con
la siguiente ecuación en diferencias:

𝑦 + = (𝑛 + 1)(𝑦 + 1) [5]

Donde 𝑦 es el número de multiplicaciones requeridas

durante el cálculo de la determinante de una matriz. Al
resolver dicha ecuación se obtiene la siguiente solución
−
𝑛! 𝑛! 𝑛! 𝑛!
𝑦(𝑛) = = 𝑛! + + + ⋯ + [6]
=
𝑘! 2! 3! (𝑛 − 1)!

Por ejemplo, para una determinante de orden 𝑛 = 4 se

requiere el siguiente número de multiplicaciones:

4! 4! 4!
𝑦(4) = = 4! + + = 40
=
𝑘! 2! 3!

Esto quiere decir, que, si deseamos resolver un sistema de

ecuaciones lineales de orden 4, se requieren calcular 5
determinantes de orden 4, lo que se traduce en 5 × 40 = 200
multiplicaciones.
Notas Análisis Numérico Grupo 26 Unidad 3 Dr. Guillermo Alberto Sánchez Lozano

Si uno desease implementar el método de Cramer para

sistemas de ecuaciones lineales grandes, llegaríamos a la
conclusión de que es un método en extremo inútil para la
resolución de sistemas algebraicos lineales. Esto se puede
apreciar en la siguiente gráfica donde se plasman el orden
del sistema de ecuaciones contra las operaciones de
multiplicación necesarias para resolverlo por Cramer:

Del previo análisis podemos apreciar que el método de

Cramer es en extremo ineficiente. Dicho esto, aunque
parezca contraintuitivo, el método de resolución más
eficiente es el denominado eliminación gaussiana, el cual es
además la antesala para el método de descomposición 𝐿𝑈 ,
el cual se emplea para factorizar una matriz.

Eliminación gaussiana
La técnica de eliminación de incógnitas se puede extender a
sistemas de cualquier orden, el cual se logra mediante un
algoritmo computacional que se encarga de realizar las
operaciones entre renglones y finalmente una etapa de
sustitución hacia atrás. El esquema más sencillo es el de
eliminación gaussiana simple. Este método contiene las
principales características de la resolución de dichos
problemas (más adelante se verá una versión mejorada del
método). Comencemos teniendo en mente el siguiente
sistema de ecuaciones:
Notas Análisis Numérico Grupo 26 Unidad 3 Dr. Guillermo Alberto Sánchez Lozano

𝑎 𝑥 + 𝑎 𝑥 + 𝑎 𝑥 + ⋯+ 𝑎 𝑥 = 𝑏
𝑎 𝑥 + 𝑎 𝑥 + 𝑎 𝑥 + ⋯+ 𝑎 𝑥 = 𝑏
𝑎 𝑥 + 𝑎 𝑥 + 𝑎 𝑥 + ⋯+ 𝑎 𝑥 = 𝑏 [7]
⋮
𝑎 𝑥 + 𝑎 𝑥 + 𝑎 𝑥 + ⋯+ 𝑎 𝑥 = 𝑏

La técnica consiste en realizar eliminación gaussiana para

colocar ceros debajo de los elementos de la diagonal
principal, para finalmente realizar la sustitución hacia atrás.

Eliminación de incógnitas hacia adelante: La primera

etapa consta en reducir el sistema de ecuaciones a uno
triangular superior. Esto se logra mediante los siguientes
pasos:

1.-Se desea eliminar la incógnita 𝑥 desde la segunda hasta

la 𝑛-ésima ecuación. Par ello, se multiplica la primera
ecuación por 𝑎 ⁄𝑎 , dando como resultado lo siguiente,
𝑎 𝑎 𝑎 𝑎
𝑎 𝑥 + 𝑎 𝑥 + 𝑎 𝑥 + ⋯+ 𝑎 𝑥 = 𝑏 (𝑖 = 2,3, … , 𝑛) [8]
𝑎 𝑎 𝑎 𝑎

La anterior ecuación se emplea de manera progresiva hasta

el último renglón, arrojando el siguiente resultado,
𝑎 𝑥 + 𝑎 𝑥 + 𝑎 𝑥 + ⋯+ 𝑎 𝑥 =𝑏
+𝑎 𝑥 + 𝑎 𝑥 + ⋯ + 𝑎 𝑥 =𝑏
+𝑎 𝑥 + 𝑎 𝑥 + ⋯ + 𝑎 𝑥 =𝑏
⋮
+𝑎 𝑥 + 𝑎 𝑥 + ⋯ + 𝑎 𝑥 =𝑏

Donde los coeficientes están dados por

𝑎
𝑎 =𝑎 − 𝑎
𝑎
𝑎 (𝑖, 𝑗 = 2,3, … 𝑛)
𝑏 =𝑏 − 𝑏
𝑎

Del procedimiento anterior, vemos que se utilizó el elemento

𝑎 para eliminar los elementos debajo de él. Este elemento
se denomina pivote, y en este caso, también su respectivo
renglón se denomina renglón pivote.
Notas Análisis Numérico Grupo 26 Unidad 3 Dr. Guillermo Alberto Sánchez Lozano

2.- El elemento pivote en esta ocasión es 𝑎 junto con su

respectivo renglón pivote, por lo que ahora se multiplicara
por 𝑎 ⁄𝑎 , es decir,
𝑎 𝑎 𝑎
𝑎 𝑥 + 𝑎 𝑥 + ⋯+ 𝑎 𝑥 = 𝑏 (𝑖 = 3,4,5, … 𝑛)
𝑎 𝑎 𝑎

La anterior ecuación se sustrae de los renglones que están

debajo de 𝑎 , lo que nos proporciona el siguiente sistema
de ecuaciones:
𝑎 𝑥 + 𝑎 𝑥 + 𝑎 𝑥 + ⋯+ 𝑎 𝑥 =𝑏
+𝑎 𝑥 + 𝑎 𝑥 + ⋯ + 𝑎 𝑥 =𝑏
+𝑎 𝑥 + ⋯ + 𝑎 𝑥 =𝑏
⋮
+𝑎 𝑥 + ⋯ + 𝑎 𝑥 =𝑏

Donde los nuevos coeficientes están dados por

𝑎
𝑎 =𝑎 − 𝑎
𝑎
(𝑖, 𝑗 = 2,3, … 𝑛)
𝑎
𝑏 =𝑏 − 𝑏
𝑎

El superíndice que aparece en las ecuaciones anteriores

indica el número de veces que el elemento (o renglón) ha
sido modificado por las operaciones de multiplicación
sustracción. Después de aplicar una metodología similar a
todo el sistema se obtiene
𝑎 𝑥 + 𝑎 𝑥 + 𝑎 𝑥 + ⋯+ 𝑎 𝑥 = 𝑏
+𝑎 𝑥 + 𝑎 𝑥 + ⋯ + 𝑎 𝑥 = 𝑏
𝑎 𝑥 + ⋯+ 𝑎 𝑥 = 𝑏 [9]
⋮
𝑎 − 𝑥 =𝑏 −

A continuación, se muestra un pseudocódigo de la

implementación de la eliminación gaussiana hacia adelante.
Notas Análisis Numérico Grupo 26 Unidad 3 Dr. Guillermo Alberto Sánchez Lozano

DOFOR 𝒌 = 𝟏, 𝒏 − 𝟏 ! Coloca ceros debajo del elemento 𝒂

DOFOR 𝒊 = 𝒌 + 𝟏, 𝒏

𝐹𝑎𝑐𝑡𝑜𝑟 = 𝑎(𝑖, 𝑘)/𝑎(𝑘, 𝑘) !Se calcula el elemento pivote

DOFOR 𝒋 = 𝒌 + 𝟏, 𝒏 ¡Se calculan los nuevos elementos de cada renglón

𝑎(𝑖, 𝑗) = 𝑎(𝑖, 𝑗) − 𝐹𝑎𝑐𝑡𝑜𝑟 ∗ 𝑎(𝑘, 𝑗)

ENDDO

𝑏(𝑖) = 𝑏 − 𝐹𝑎𝑐𝑡𝑜𝑟 ∗ 𝑏(𝑘)

ENDDO

ENDDO

Tabla 1. Pseudocódigo que realiza la eliminación hacia adelante

(basado en la figura 9.4a de Métodos Numéricos para Ingenieros 7ª
edición).

El ciclo exterior recorre el elemento pivote de la diagonal

para efectuar las operaciones hacia abajo durante la
eliminación gaussiana. El ciclo intermedio obtiene el
elemento pivote con el cual se colocarán ceros debajo de
la diagonal. El ciclo interior realiza las modificaciones de
los 𝒂(𝒊, 𝒋) y los 𝒃 sobre los demás elementos del renglón
𝒊. Lo anterior permite obtener un sistema de ecuaciones
lineales, el cual adopta la forma de un sistema triangular
superior.

Sustitución hacia atrás: De la matriz triangular superior

resultante [10] se puede inferir del último renglón que
( − )
𝑏
𝑥 = ( − )
𝑎

El resultado anterior, se sustituye en la ecuación inmediata

superior, lo cual permite obtener 𝑥 − . Este procedimiento
se repite varias veces para calcular el resto de incógnitas 𝑥
Notas Análisis Numérico Grupo 26 Unidad 3 Dr. Guillermo Alberto Sánchez Lozano

(𝑖 = 𝑛 − 1, 𝑛 − 2, … ,3,2,1). Tal proceso se sintetiza en la

siguiente fórmula:
(− ) (− )
𝑏 −∑ = +
𝑎 𝑥
𝑥 = ( − )
[10]
𝑎

La sustitución hacia atrás se puede ejecutar mediante el

siguiente pseudocódigo.

𝑥(𝑛) = 𝑏(𝑛)/𝑎(𝑛, 𝑛)
DOFOR 𝒊 = 𝒏 − 𝟏, 𝟏, −𝟏

𝑠𝑢𝑚 = 𝑏(𝑖)
DOFOR 𝒋 = 𝒊 + 𝟏, 𝒏

𝑠𝑢𝑚 = 𝑠𝑢𝑚 − 𝑎(𝑖, 𝑗) ∗ 𝑥(𝑗)

ENDDO

𝑥(𝑖) = 𝑠𝑢𝑚/𝑎(𝑖, 𝑖)
ENDDO

Tabla 2. Pseudocódigo que realiza la sustitución hacia atrás (basado en

la figura 9.4 b) de Métodos Numéricos para Ingenieros 7ª edición).
Notas Análisis Numérico Grupo 26 Unidad 3 Dr. Guillermo Alberto Sánchez Lozano

Ejercicio 1. Emplee la eliminación de Gauss simple para

resolver:
3𝑥 − 0.1𝑥 − 0.2𝑥 = 7.85
0.1𝑥 + 7𝑥 − 0.3𝑥 = −19.3
0.3𝑥 − 0.2𝑥 + 10𝑥 = 71.4

Solución: Empleando el código de la tabla 3 se obtienen

los siguientes resultados,
Notas Análisis Numérico Grupo 26 Unidad 3 Dr. Guillermo Alberto Sánchez Lozano

PROGRAM Gauss_simple
IMPLICIT NONE
INTEGER,PARAMETER:: n=3
REAL (KIND=8)::A(1:n,1:n),b(1:n),x(1:n)
INTEGER i,j
CHARACTER (LEN=6) LONGITUD_RENGLON
!El comando WRITE se puede emplear para transformar un "n" entero, en un "n" de texto
WRITE(LONGITUD_RENGLON,'(I6)') n

!Matriz de coeficientes A
A(1,1)=3.D0;A(1,2)=-0.1D0;A(1,3)=-0.2D0
A(2,1)=0.1D0;A(2,2)=7.D0;A(2,3)=-0.3D0
A(3,1)=0.3D0;A(3,2)=-0.2D0;A(3,3)=10.D0
!Vector constante
b(1)=7.85D0;b(2)=-19.3D0;b(3)=71.4D0

!IMPRESION PRELIMINAR DEL SISTEMA DE ECUACIONES

WRITE(*,*) 'Sistema de ecuaciones a resolver Ax=b'
WRITE(*,*)

DO i=1,n,1
WRITE(*,'('//LONGITUD_RENGLON//'(f17.8,x),A1,f17.8)') (a(i,j),j=1,3,1),'=',b(i) !DO IMPLICITO
END DO

WRITE(*,*)
WRITE(*,*)
!Eliminación gaussiana sobre el sistema Ax=b
CALL ELIMINACION_GAUSSIANA(A,b,n)
!Sus tución hacia atrás sobre la matriz triangular
CALL SUSTITUCION_ATRAS(A,b,x,n)
!Impresión de resultados de las x(i)
WRITE(*,*) 'Soluciones al sistema:'
WRITE(*,*)
DO i=1,n,1
WRITE(*,'(A2,I3,A2,f17.8)')'x(',i,')=', x(i)
END DO
END PROGRAM
!*****SUBRUTINA ELIMINACIÓN GAUSSIANA
SUBROUTINE ELIMINACION_GAUSSIANA(A,b,n)
IMPLICIT NONE
REAL (KIND=8)::A(1:n,1:n),b(1:n)
REAL (KIND=8) pivote
INTEGER n,k,i,j
DO k=1,n-1,1 !Mueve hacia abajo el elemento (renglón) pivote
DO i=k+1,n,1
pivote=a(i,k)/a(k,k) !se calcula el pivote
DO j=k+1,n,1
a(i,j)=a(i,j)-pivote*a(k,j)
END DO
b(i)=b(i)-pivote*b(k)
END DO
END DO
END SUBROUTINE
!************************************
!***********SUBRUTINA PARA SUSTITUCIÓN HACIA ATRÁS
SUBROUTINE SUSTITUCION_ATRAS(a,b,x,n)
IMPLICIT NONE
REAL (KIND=8)::A(1:n,1:n),b(1:n),x(1:n)
REAL (KIND=8) suma
INTEGER n,i,j

x(n)=b(n)/a(n,n)

DO i=n-1,1,-1
suma=b(i)
DO j=i+1,n,1
suma=suma-a(i,j)*x(j)
END DO
x(i)=suma/a(i,i)
END DO
END SUBROUTINE

Tabla 3. Eliminación de Gauss simple con las subrutinas de

eliminación hacia adelante y sustitución hacia atrás, basados en los
pseudocódigos de las tablas 1 y 2, para el ejercicio 1.
Notas Análisis Numérico Grupo 26 Unidad 3 Dr. Guillermo Alberto Sánchez Lozano

Técnicas para mejorar soluciones

El método de Gauss simple tiene ciertos inconvenientes. Si
un elemento pivote (o en este caso, un elemento de la
diagonal principal) 𝑎 = 0, 𝑎 ∼ 0 o 𝑎 ≪ 𝑎 se
producirán errores debido a divisiones entre cero (o
redondeos próximos a cero).

Si se llega a la conclusión de que det𝐴 ∼ 0, el sistema se

dice que está mal condicionado. Debe decirse además, que
la det𝐴 es un indicador no muy confiable para estos casos.

Nota: Una estrategia para saber si una matriz 𝐴 ∈

𝑀 × sistema está bien condicionado es calcular el
número de condición , el cual se define de la siguiente
manera:

Cond𝐴 = ‖𝐴‖ ⋅ ‖𝐴− ‖ [11]

Donde ‖𝐴‖ se define de la siguiente manera:

‖𝐴‖ = 𝑎 [12]
= =

La anterior definición se denomina Norma de Frobenius 1.

Si Cond𝐴 < 1 se dice que el sistema está bien
condicionado, por el contrario, si Cond𝐴 > 1 se dice que el
sistema está mal condicionado.

Una de las estrategias para evitar estas situaciones

consiste en cambiar el orden de los renglones evitando que
los elementos pivote sean cero o cercanos a cero, lo cual se
denomina pivoteo parcial . Para ello es conveniente

1
Existen diversas normas para las matrices adicionales a la de Frobenius, entre ellas destacan:
 La norma p.
 Magnitud máxima.
 Columna suma.
 Renglón suma.
 Norma espectral.
Notas Análisis Numérico Grupo 26 Unidad 3 Dr. Guillermo Alberto Sánchez Lozano

normalizar la matriz, es decir, dividir cada renglón (o

columna por el elemento máximo en él).

Existe además lo que se denomina pivoteo total , lo cual

consiste en cambiar el orden de los renglones y de las
columnas. No obstante, este tipo de pivoteo requiere una
cantidad ingente de operaciones lo cual puede sumar una
complejidad enorme durante lo programación.
Normalmente, el pivoteo parcial proporciona buenos
resultados y resuelve el problema del mal
condicionamiento.

Debe hacerse la nota de que realmente el escalamiento en

este caso es un indicador para determinar si el pivoteo
parcial es necesario, aunque para encontrar la solución no
es necesario escalar el sistema de ecuaciones. Es decir, el
escalamiento ayuda a determinar si es necesario el
pivoteo.

En la tabla 4 se muestra una versión mejorada del

algoritmo de la eliminación de Gauss a la previamente
mostrada en tabla 3 Las mejoras consisten en :

 Las ecuaciones no están escaladas, pero los valores

escalados de los renglones permiten determinar si es
necesario usar pivoteo parcial. Esto se logra
almacenando el elemento más grande de cada renglón
𝒔(𝒋) (𝒋 = 𝟏, 𝟐, … , 𝒏), el cual se emplea para normalizar
el 𝒋-ésimo renglón.

 Se vigila durante todo el proceso de eliminación que los

elementos pivotes (la diagonal principal) son cercanos a
cero y determinar si el sistema es singular (en el código,
si 𝑒𝑟 = −1 ∴ el sistema es singular). Para esto se prefija
una tolerancia 𝑡𝑜𝑙 con una pequeña magnitud para
detectar elementos cercanos a cero.
Notas Análisis Numérico Grupo 26 Unidad 3 Dr. Guillermo Alberto Sánchez Lozano

PROGRAM GAUSS_PIVOTEOPARCIAL_ESCALAMIENTO
IMPLICIT NONE
INTEGER, PARAMETER:: n=5 !orden del sistema de ecuaciones lineales
REAL (KIND=8)::a(1:n,1:n),b(1:n),x(1:n),s(1:n) !Asignación de las respec vas dimensiones de A=a(i,j);b(i),x(i)
INTEGER er,i,j
CHARACTER (LEN=4) CONTADOR
CHARACTER (Len=5) LONGITUD_RENGLON
!Declaración de la matriz de coeficientes A=a(i,j)
a(1,1)=6.1d0;a(1,2)=-2.4d0;a(1,3)=23.3d0;a(1,4)=-16.4d0;a(1,5)=-8.9d0
a(2,1)=-14.2d0;a(2,2)=-31.6d0;a(2,3)=-5.8d0;a(2,4)=9.6d0;a(2,5)=23.1d0
a(3,1)=10.5d0;a(3,2)=46.1d0;a(3,3)=-19.6d0;a(3,4)=-8.8d0;a(3,5)=-41.2d0
a(4,1)=37.3d0;a(4,2)=-14.2d0;a(4,3)=62.d0;a(4,4)=14.7d0;a(4,5)=-9.6d0
a(5,1)=0.8d0;a(5,2)=17.7d0;a(5,3)=-47.5d0;a(5,4)=-50.2d0;a(5,5)=29.8d0
!Declaración del vector de constantes B=b(i)
b(1)=121.7d0;b(2)=-87.7d0;b(3)=10.8d0;b(4)=61.3d0;b(5)=-27.8d0
!Se convierte el "n" númerico en texto y se almacena en Longitud_renglon para
!poder mandar a pantalla con el formato establecido

WRITE(LONGITUD_RENGLON,'(I3)') n

!Solo para visualizar los datos de entrada

WRITE(*,*)'Matriz A'
DO i=1,n,1
WRITE(*,'('//LONGITUD_RENGLON//'(f13.8,2x))') (a(i,j),j=1,n) !DO IMPLICITO
END DO
WRITE(*,*)'Vector de coeficientes B'
DO i=1,n,1
WRITE(*,'(f13.8)') b(i) !DO IMPLICITO
END DO

!Eliminación Gaussiana del sistema de ecuaciones a la forma triangular

CALL Gauss(a,b,s,n)
!Obtención del vector solución x(i) a par r de sus tución hacia atrás
CALL sus tucion_atras(a,b,x,n)
!IMPRESIÓN DE RESULTADOS
DO i=1,n
WRITE(CONTADOR,'(I4)')i
WRITE(*,'(A8,f10.5)') 'x('//CONTADOR//')=',x(i)
END DO
ENDPROGRAM

Tabla 4.1 Eliminación de Gauss con pivoteo parcial y escalamiento,

ajustado al ejemplo 5 (basado en el pseudocódigo de la figura 9.6 de
Métodos Numéricos para Ingenieros 7ª edición), para el ejercicio 2.
Notas Análisis Numérico Grupo 26 Unidad 3 Dr. Guillermo Alberto Sánchez Lozano

***************Subru na para elminación Gaussiana con pivoteo y escalamiento

!*****

SUBROUTINE Gauss(a,b,s,n)
IMPLICIT NONE
REAL (KIND=8),PARAMETER:: tol=10.D0**(-12) !Tolerancia
REAL (KIND=8)::a(1:n,1:n),b(1:n),x(1:n),s(1:n)
INTEGER i,j,n,er
!s(j) almacena los elementos los "n" elementos a(i,1), es decir, la primera columna (véase [1]),
!para posteriormente compararlos con los a(i,j) (j=2,3,4,...,n) del renglón "i"
!y determinar cual es el mayor, es decir, toma el elemento máximo de cada renglón (véase el ciclo [2])

!"er" es una variable entera que ene la función de adver r si algún elemento de la diagonal (pivote)
!es menor a una tolerancia, es decir, a(i,i)<tol, donde "tol" es un valor muy pequeño y posi vo
!que hace las veces de "cero".
er=0
DO i=1,n
s(i)=DABS(a(i,1)) ![1]
DO j=2,n ![2]
IF (DABS(a(i,j)).GT.s(i)) THEN
s(i)=DABS(a(i,j))
ENDIF
ENDDO
ENDDO
CALL Eliminacion(a,b,s,n,er,tol) !Subru na que realiza la eliminación Gaussiana
!(triangular superior)
IF (er.NE.-1) THEN !Se realiza la sus tución hacia atrás si
CALL sus tucion_atras(a,b,x,n)
ENDIF
ENDSUBROUTINE

Tabla 4.2 Eliminación de Gauss con pivoteo parcial y escalamiento,

ajustado al ejemplo 5 (basado en el pseudocódigo de la figura 9.6 de
Métodos Numéricos para Ingenieros 7ª edición).
Notas Análisis Numérico Grupo 26 Unidad 3 Dr. Guillermo Alberto Sánchez Lozano

!***********************************************************
SUBROUTINE Eliminacion(a,b,s,n,er,tol)
IMPLICIT NONE
REAL (KIND=8)::a(1:n,1:n),b(1:n),s(1:n)
REAL (KIND=8):: tol,pivote
INTEGER k,i,j,n,er
DO k=1,n-1,1 !desplaza las operaciones al k-ésimo renglón, esto se nota para i=k+1, que realmente i[2,n]
CALL PIVOTEO_PARCIAL_ESCALAMIENTO(k,a,b,s,n) !pivoteo realizado según sus elementos normalizados
IF (DABS(a(k,k)/s(k)).LT. tol) THEN !Prueba si el elemento a(k,k) NORMALIZADO por s(k) es cercano a CERO
er=-1
GOTO 100!termina el ciclo si detecta un elemento de la diagonal (pivote NORMALIZADO) cercano a CERO
END IF
DO i=k+1,n !Se realiza la eliminación Gaussiana volviendo la matriz de coeficientes en
!una triangular superior
pivote=a(i,k)/a(k,k)
DO j=k+1,n,1
a(i,j)=a(i,j)-pivote*a(k,j)
ENDDO
b(i)=b(i)-pivote*b(k)
ENDDO
ENDDO
100 IF (DABS(a(n,n)/s(n)).LT.tol) THEN !
er=-1
WRITE(*,*) 'Matriz singular'
END IF
ENDSUBROUTINE

Tabla 4.3 Eliminación de Gauss con pivoteo parcial y escalamiento,

ajustado al ejemplo 5 (basado en el pseudocódigo de la figura 9.6 de
Métodos Numéricos para Ingenieros 7ª edición).
Notas Análisis Numérico Grupo 26 Unidad 3 Dr. Guillermo Alberto Sánchez Lozano

!*****************************Subru na para pivoteo parcial con Escalamiento

SUBROUTINE PIVOTEO_PARCIAL_ESCALAMIENTO(k,a,b,s,n) !pivoteo realizado según sus elementos normalizados
IMPLICIT NONE
REAL (KIND=8):: a(1:n,1:n),b(1:n),s(1:n)
REAL (KIND=8) a_mayor, tere
INTEGER ii,jj,k,n,p
p=k !Posición del pivote en la columna k
a_mayor=ABS(a(k,k)/s(k)) !Elemento pivote NORMALIZADO en la columna k

!Este ciclo compara el pivote NORMALIZADO a(k,k)/s(k) con los elementos también NORMALIZADOS debajo de él
!y determina cual de todos es el mayor
DO ii=k+1,n
tere=ABS(a(ii,k)/s(ii))
IF ( tere .GT. a_mayor) THEN
a_mayor= tere
p=ii
ENDIF
ENDDO
!La siguiente condición intercambia los renglones respecto al elemento más grande
!almacenado al final del ciclo anterior (junto con su posición), si es que el pivote original
!no es el de mayor magnitud
IF (p .NE. k) THEN !intercambio del renglón "k" con el renglón "j"

DO jj=k,n !justo en este ciclo se intercambian cada uno de los elementos

!del renglón a(k,jj) con el mayor a(p,jj)
tere=a(p,jj)
a(p,jj)=a(k,jj)
a(k,jj)= tere
END DO
tere=b(p) !Se intercambian b(p) y b(k) de posición
b(p)=b(k)
b(k)= tere
!En este momento se intercambian los elementos s(k) por s(p)
!(recuerde que s(j) (j=1,2,...,n) es la magnitud
!el elemento más grande de cada renglón
tere=s(p)
s(p)=s(k)
s(k)= tere
END IF
ENDSUBROUTINE
!********************Subru na para sus tución hacia atrás
SUBROUTINE sus tucion_atras(a,b,x,n)
IMPLICIT NONE
REAL (KIND=8)::a(1:n,1:n),b(1:n),x(1:n)
REAL (KIND=8) sum
INTEGER i,j,n
x(n)=b(n)/a(n,n)
DO i=n-1,1,-1
sum=0.d0
DO j=i+1,n,1
sum=sum+a(i,j)*x(j)
ENDDO
x(i)=(b(i)-sum)/a(i,i)
ENDDO
ENDSUBROUTINE

Tabla 4.4 Eliminación de Gauss con pivoteo parcial y escalamiento,

ajustado al ejemplo 5 (basado en el pseudocódigo de la figura 9.6 de
Métodos Numéricos para Ingenieros 7ª edición).
Notas Análisis Numérico Grupo 26 Unidad 3 Dr. Guillermo Alberto Sánchez Lozano

Ejercicio 2. Emplee la eliminación de Gauss con pivoteo

parcial para resolver el siguiente sistema
6.1 −2.4 23.3 −16.4 −8.9 𝑥 121.7
⎛−14.2 −31.6 −5.8 9.6 23.1 ⎞ ⎛𝑥 ⎞ ⎛−87.7 ⎞
⎜
⎜ ⎟
⎟ ⎜
⎜ ⎟
⎟ ⎜
⎜ ⎟
⎟
⎜ 46.1 −19.6 −8.8 −41.2⎟ 𝑥 =
⎜ 10.5
⎜ 37.3 −14.2 ⎟⎜⎜ ⎟
⎟ ⎜ 10.8
⎟ ⎜ 61.3 ⎟
⎜ ⎟
62 14.7 −9.6 ⎟ ⎜𝑥 ⎟
⎝ 0.8 17.7 −47.5 −50.2 29.8 ⎠ ⎝𝑥 ⎠ ⎝−27.8⎠

Solución

El sistema de ecuaciones fue resuelto mediante el código de la

tabla 4, cuyos resultados son:
Notas Análisis Numérico Grupo 26 Unidad 3 Dr. Guillermo Alberto Sánchez Lozano

Descomposición
La eliminación gaussiana tradicional requiere una cantidad
considerable de operaciones entre renglones. A continuación,
se estudiará una técnica denominada descomposición LU la
cual efectúa principalmente operaciones sobre la matriz de
coeficientes permitiendo obtener una factorización de esta.
De manera tradicional, el método de eliminación gaussiana
se encarga de resolver sistemas lineales de ecuaciones

𝐴𝑋 = 𝐵 [13]

Donde 𝐴 ∈ 𝑀 × y 𝑋, 𝐵 ∈ 𝑀 × . Al final de la aplicación

gaussiana sobre [9] se obtiene un sistema triangular superior
𝑢 𝑢 𝑢 … 𝑢 𝑥 𝑑
⎛ 0 𝑢 𝑢 … 𝑢 ⎞ ⎛ 𝑥 ⎞ ⎛
⎜ 𝑑 ⎞
⎟
⎜
⎜ ⎟
⎟ ⎜
⎜ ⎟
⎟ ⎜ ⎟
⎜ 0 0 𝑢 … 𝑢 ⎟ ⎜ 𝑥 ⎟ ⎜
= ⎜𝑑 ⎟
⎟ [14]
⎜
⎜ ⋮ ⎟
⎟ ⎜
⎜ ⎟
⎟ ⎜ ⎟
⋮ ⋮ ⋱ ⋮ ⋮ ⎜ ⋮ ⎟
⎝ 0 0 0 … 𝑢 ⎠ ⎝𝑥 ⎠ ⎝𝑑 ⎠

Lo cual de manera sintetizada se puede expresar como

𝐴𝑋 = 𝐵 ⇒ 𝑈 𝑋 = 𝐷

Suponga además que existe una matriz triangular inferior

con la siguiente estructura
1 0 0 0 0
⎛
⎜ 𝑙 1 0 0 0⎞
⎟
⎜
⎜ ⎟
𝐿 = ⎜𝑙 𝑙 1 0 0⎟
⎟ [15]
⎜
⎜ ⋮ ⎟
⋮ ⋮ ⋱ ⋮⎟
⎝𝑙 𝑙 𝑙 … 1⎠

Se puede demostrar que toda matriz no singular se puede

descomponer como el siguiente producto matricial

𝐴 = 𝐿𝑈 [16]

La matriz 𝑈 se obtiene de la siguiente manera

𝑈 =𝐸 𝐸 − …𝐸 𝐸 𝐴 [17]
Notas Análisis Numérico Grupo 26 Unidad 3 Dr. Guillermo Alberto Sánchez Lozano

Donde 𝐸 (𝑖 = 1,2, … , 𝑛) son matrices elementales 2. De la

anterior ecuación

𝐸 − 𝐸 − … 𝐸 −− 𝐸 − 𝑈 = 𝐸 − 𝐸 − … 𝐸 −− 𝐸 − 𝐸 𝐸 − …𝐸 𝐸 𝐴

Simplificando

𝐸 − 𝐸 − … 𝐸 −− 𝐸 − 𝑈 = 𝐴

Es decir,

𝐿 = 𝐸 − 𝐸 − … 𝐸 −− 𝐸 − [18]

Lo anterior significa que la matriz 𝐿 se obtiene almacenando

las operaciones inversas que se utilizaron para la obtención
de la matriz 𝑈.

Ejercicio 3. Obtenga la factorización 𝐿𝑈 de la siguiente

matriz
1 2 3
𝐴= 4 5 6 [𝐸1]
7 8 4
Solución

Debemos obtener la factorización 𝐿𝑈 de la matriz dada en

[E1], para ello debemos convertir esta en una triangular
superior mediante eliminación gaussiana:
→ −
1 2 3 → − 1 2 3 → − 1 2 3
4 5 6 0 −3 −6 0 −3 −6
7 8 4 0 −6 −17 0 0 −5
De lo anterior, deducimos que

2
Una matriz elemental es aquella que codifica en sus elementos operaciones en matrices tales como:
permutación, sustracción (o adición) o escalamiento de renglones. Una bibliografía recomendada es
Álgebra Lineal de Stanley Grossman.
Notas Análisis Numérico Grupo 26 Unidad 3 Dr. Guillermo Alberto Sánchez Lozano

1 0 0 1 0 0 1 0 0 1 2 3 1 2 3
0 1 0 0 1 0 −4 1 0 4 5 6 = 0 −3 −6
0 −2 1 −7 0 1 0 0 1 7 8 4 0 0 −5

Por lo tanto, la matriz 𝐿 se obtiene a partir del producto

de las inversas de las matrices elementales.

𝐿 = 𝐸− 𝐸− 𝐸−

De este modo:

1 0 0 1 0 0 1 0 0 1 0 0
𝐿= 4 1 0 0 1 0 0 1 0 = 4 1 0
0 0 1 7 0 1 0 2 1 7 2 1

Ejercicio 3.1 Obtenga la descomposición 𝐿𝑈 de la siguiente

matriz:
1 2 3 5
⎛−3 3 2 −3⎞
𝐴=⎜
⎜4 ⎟
6 4 8⎟
⎝−5 −10 6 11 ⎠

Solución

Como recomendación coloque la matriz identidad a la

izquierda de 𝐴:
1 0 0 0 1 2 3 5
⎛
⎜0 1 0 ⎞ ⎛
0⎟ ⎜−3 3 2 −3⎞⎟
⎜0 0 1 0⎟ ⎜ 4 6 4 8⎟
⎝0 0 0 1⎠ ⎝−5 −10 6 11 ⎠
Los elementos 𝑎 , 𝑎 y𝑎 se utilizan para calcular los
ceros debajo de estos. Para ello simplemente se divide
𝑎
𝑓 = 𝑖>𝑗
𝑎

Posteriormente, colocamos ceros en 𝐴 de la siguiente forma

Notas Análisis Numérico Grupo 26 Unidad 3 Dr. Guillermo Alberto Sánchez Lozano

1 0 0 0 1 2 3 5
⎛
⎜ ⎞ ⎛
−3 1 0 0 ⎜0 9 11 12 ⎞
⎜ 4 0 1 0⎟
⎟ ⎜0 −2 −8
⎟
−12⎟
⎝−5 0 0 1⎠ ⎝0 0 21 36 ⎠
Como ahora 𝑎 = 9 es el pivote, usamos este para colocar
ceros
1 0 0 0 1 2 3 5
⎛
⎜−3 1 0 ⎞ ⎛
0⎟ ⎜0 9 11 12 ⎞⎟
⎜
⎜4 ⎟⎜ ⎟
− 2⁄9 1 0⎟ ⎜0 0 − 50⁄9 − 28⁄3⎟
⎝−5 0 0 1⎠ ⎝0 0 21 36 ⎠
Ahora se emplea el elemento pivote 𝑎 = − 50⁄9, de este
modo:
1 0 0 0 1 2 3 5
⎛
⎜−3 1 0 0⎞⎟⎛
⎜0 9 11 12 ⎞⎟
⎜
⎜4 − 2⁄9 1 0⎟⎟⎜
⎜0 0 − 50⁄9 − 28⁄3⎟
⎟
⎝−5 0 − 189⁄50 1⎠ ⎝ 0 0 0 18⁄25 ⎠

Eliminación de Gauss para calcular la

descomposición 𝑨 = 𝑳𝑼
La eliminación de Gauss puede usarse directamente para
descomponer 𝐴 en 𝐿𝑈 . Recuerde qué 𝐴 puede
transformarse mediante operaciones entre renglones como
𝑎 𝑎 𝑎 … 𝑎 ( − ) 𝑎
⎛
⎜ 0 𝑎 𝑎 … 𝑎 𝑎 ⎞
⎟
( − )
⎜
⎜ ⎟
⎟
⎜ 0 0 𝑎 … 𝑎( −) 𝑎 ⎟
𝑈=⎜
⎜
⎜ ⋮ ⋮ ⋮ ⋱ ⋮ ⋮ ⎟
⎟
⎟ [19]
⎜
⎜ ( − ) ( − ) ⎟
⎟
⎜ 0
⎜ 0 0 … 𝑎( − )( − )
𝑎( − ) ⎟⎟
( − )
⎝ 0 0 0 … 0 𝑎 ⎠

Donde el superíndice que aparece indica las veces que

fue modificado el renglón debido a las operaciones de
eliminación de Gauss . La matriz 𝐿 se calcula de la
siguiente forma
Notas Análisis Numérico Grupo 26 Unidad 3 Dr. Guillermo Alberto Sánchez Lozano

1 0 0 … 0 0
⎛
⎜ 𝑓 1 0 … 0 0⎞⎟
⎜
⎜ ⎟
𝑓 𝑓 1 … 0 0⎟
𝐿=⎜
⎜
⎜ ⋮ ⋮ ⋮ ⋱ ⋮ ⋮⎟
⎟
⎟ [20]
⎜
⎜ ⎟
𝑓
⎜ ( − ) 𝑓( − ) 𝑓( − ) … 1 0⎟⎟
⎝ 𝑓 𝑓 𝑓 … 𝑓 ( − ) 1⎠

Donde 𝑓 se obtiene del elemento diagonal de la matriz

original (elemento pivote) que fue usado para colocar un
cero en la posición 𝑖𝑗.
𝑎
𝑓 =
𝑎 𝑎
𝑎 𝑓 =
𝑓 = 𝑎 ⋮
𝑎 ⋮
⋮ , ,𝑓 𝑎( − ) ,…,
𝑎( − 𝑎( − ) ( − ) =
𝑓( =
) 𝑓( − ) = 𝑎 𝑎
( − )
− )
𝑎 𝑎 ( − )
𝑎 𝑓 ( − ) = ( − )
𝑎 𝑎 𝑓 = 𝑎( − )( − )
𝑓 = 𝑓 = 𝑎
𝑎 𝑎

Simplificando, los factores anteriores pueden obtenerse

mediante la siguiente expresión:
( − )
𝑎 𝑖 = 2,3, … , 𝑛
𝑓 = con y 𝑖>𝑗 [21]
𝑎
( − ) 𝑗 = 1,2, … , 𝑛 − 1

Cabe aclarar que cada operación del factor 𝑓 para colocar

ceros en 𝐴 es en síntesis equivalente a multiplicar por
una matriz elemental que sustrae el renglón 𝒊 a los
consecuentes.

El almacenamiento eficiente (computacionalmente) de la

descomposición 𝐿𝑈 se sintetiza en la siguiente matriz:
Notas Análisis Numérico Grupo 26 Unidad 3 Dr. Guillermo Alberto Sánchez Lozano

𝑎 𝑎 𝑎 … 𝑎 ( − ) 𝑎
⎛
⎜ 𝑓 𝑎 𝑎 … 𝑎 𝑎 ⎞
⎟
( − )
⎜
⎜ ⎟
⎟
⎜ 𝑓 𝑓 𝑎 … 𝑎 ( − ) 𝑎 ⎟
⎜
⎜ ⎟ [22]
⎜ ⋮ ⋮ ⋮ ⋱ ⋮ ⋮ ⎟ ⎟
⎜
⎜ ( − ) ⎟
𝑎( − ) ⎟
( − )
⎜𝑓 𝑓( 𝑓( … 𝑎( ⎟
⎜ ( − ) − ) − ) − )( − ) ⎟
( − )
⎝ 𝑓 𝑓 𝑓 … 𝑓 ( − ) 𝑎 ⎠

La matriz anterior puede conseguirse mediante el siguiente

pseudocódigo para la fase de descomposición (análoga a la
tabla 1), donde se almacenan al mismo tiempo los
elementos 𝑓 de la matriz 𝐿.

Tipos de descomposición LU
La descomposición 𝐿𝑈 con la eliminación de Gauss requiere
que la matriz triangular inferior [𝐿] tenga en su diagonal
únicamente 1s. Formalmente, esta descomposición se conoce
como descomposición de Doolitle. El método alternativo
consiste en colocar 1s en la diagonal principal de la matriz
triangular superior [𝑈 ], a esta descomposición se le conoce
como descomposición de Crout . Existen algunas diferencias
entre ambos métodos (algoritmo) pero su funcionamiento es
muy similar. Lo anterior se resume de la siguiente manera:

Descomposición Doolitle
1 0 0 … 0 𝑢 𝑢 𝑢 … 𝑢
⎡𝑙 1 0 … 0⎤ ⎡ 0 𝑢 𝑢 … 𝑢 ⎤
⎢ … 0⎥ ⎢ ⎥
[𝐿][𝑈 ] = ⎢ 𝑙 𝑙 1 ⎥⎢ 0 0 𝑢 … 𝑢 ⎥
⎢ ⋮ ⋮ ⋮ ⋱ 0⎥ ⎢ ⋮ ⋮ ⋮ ⋱ 𝑢 − ⎥
⎣𝑙 𝑙 𝑙 … 1⎦ ⎣ 0 0 0 … 𝑢 ⎦

Descomposición Crout
𝑙 0 0 … 0 1 𝑢 𝑢 … 𝑢
⎡𝑙 𝑙 0 … 0 ⎤ ⎡0 1 𝑢 … 𝑢 ⎤
⎢ ⎥⎢ ⎥
[𝐿][𝑈 ] = ⎢ 𝑙 𝑙 𝑙 … 0 ⎥ ⎢0 0 1 … 𝑢 ⎥
⎢ ⋮ ⋮ ⋮ ⋱ 0 ⎥⎢⋮ ⋮ ⋮ ⋱ 𝑢 − ⎥
⎣𝑙 𝑙 𝑙 … 𝑙 ⎦ ⎣0 0 0 … 1 ⎦
Notas Análisis Numérico Grupo 26 Unidad 3 Dr. Guillermo Alberto Sánchez Lozano

Ejercicio 4. Obtenga la descomposición 𝐿𝑈 de la matriz

2 −7 1
8 −6 9 [𝐸1]
−1 3 5
Solución

Realizamos operaciones entre renglones para llevar a una

matriz triangular superior
→ −
2 −7 1 2 −7 1
2 −7 1 → + → +
⎛
⎛0 22 5 ⎞ 0 22 5 ⎞
8 −6 9 ⎜
⎜ ⎟ 𝑈 =⎜ ⎟
1 11⎟ ⎜ 247⎟
−1 3 5 0 − 0 0
⎝ 2 2⎠ ⎝ 44 ⎠

La matriz 𝐿 codifica las operaciones entre renglones, pero

con los signos contrarios, es decir es la siguiente

1 0 0
⎛ 4 1 0⎞
𝐿=⎜
⎜ 1 ⎟
⎟
1
− − 1
⎝ 2 44 ⎠

Ejercicio 5 Obtenga la descomposición 𝐿𝑈 de Doolitle

para la siguiente matriz,
3 −0.1 −0.2
𝐴= 0.1 7 −0.3
0.3 −0.2 10
En la tabla 5 se presenta un código para realizar la
descomposición de Doolitle de la matriz 𝐴, cuya
descomposición de muestra a continuación:
Notas Análisis Numérico Grupo 26 Unidad 3 Dr. Guillermo Alberto Sánchez Lozano
PROGRAM DescomposicionLU
IMPLICIT NONE
INTEGER, PARAMETER::n=3 !Orden de la matriz
REAL (KIND=8)::a(1:n,1:n)
INTEGER i,j
CHARACTER (LEN=5):: TAMANO

WRITE(TAMANO,'(I5)') n !convierte el "n" numérico en "n" texto

!Declaración de la matriz
a(1,1)=3.d0;a(1,2)=-0.1d0;a(1,3)=-0.2d0
a(2,1)=0.1d0;a(2,2)=7.d0;a(2,3)=-0.3d0
a(3,1)=0.3d0;a(3,2)=-0.2d0;a(3,3)=10.d0
!Impresión de la matriz original
WRITE(*,*) ' Matriz [A] original'

DO i=1,n,1
WRITE(*,'('//TAMANO//'(F12.7,2x))') (a(i,j),j=1,n,1)
END DO

!Se ob ene la descomposicion A=LU

CALL DESCOMPOSICION(a,n)

!Impresion de A=LU
WRITE(*,*)
WRITE(*,*) 'Matriz [A] con [L] y [U] codificados'

DO i=1,n,1
WRITE(*,'('//TAMANO//'(F12.7,2x))') (a(i,j),j=1,n,1)
END DO

END PROGRAM
!********************
SUBROUTINE DESCOMPOSICION(a,n)
IMPLICIT NONE
REAL (KIND=8)::a(1:n,1:n),factor
INTEGER k,i,j,n

DO k=1,n-1,1 !Mueve la columna "k" y pone ceros debajo de a(k,k)

DO i=k+1,n,1 !ob ene los elementos f(i,k)
factor=a(i,k)/a(k,k) !f(i,k)
a(i,k)=factor !se almacena f(i,k) en la posición original de a(i,k)
DO j=k+1,n,1 !Realiza las modificaciones sobre el renglón i
a(i,j)=a(i,j)-factor*a(k,j)
END DO
END DO
END DO
END SUBROUTINE

Tabla 5. Programa en Fortran90 para la descomposición LU.

Implementado en el ejercicio 5.
Notas Análisis Numérico Grupo 26 Unidad 3 Dr. Guillermo Alberto Sánchez Lozano

Método para resolver ecuaciones lineales

mediante el método de descomposición 𝑳𝑼
La descomposición 𝐿𝑈 es ad hoc a situaciones donde se
necesita resolver reiteradamente 𝐴𝑋 = 𝐵 pero con distintos
𝐵. Como ya se mencionó antes, la ecuación 𝐴𝑋 = 𝐵
(expandida en [7]) se puede llevar mediante eliminación
gaussiana a una forma triangular superior, la cual es

𝑈𝑋 = 𝐷 [23]

Como se detalló en la ecuación [16] es posible factorizar la

matriz 𝐴 en el producto de una matriz triangular superior
e inferior, es decir, 𝐴 = 𝐿𝑈 (mostradas en las ecuaciones
[14] y [15]). Si multiplicamos [19] por la matriz [𝐿] en ambos
lados, obtenemos:

𝐿𝑈 𝑋 = 𝐿𝐷 [24]

El trabajo anterior parece inútil ya aparentemente se regresa

al sistema original, sin embargo, las ecuaciones [19] y [20]
nos permite obtener dos sistemas de ecuaciones lineales:

𝐿𝐷 = 𝐵 [25]
𝑈𝑋 = 𝐷
Ambas ecuaciones matriciales nos permitirán obtener un
esquema para resolver el sistema 𝐴𝑋 = 𝐵 lineales tal como
se describe en el siguiente algoritmo:
Notas Análisis Numérico Grupo 26 Unidad 3 Dr. Guillermo Alberto Sánchez Lozano

Algoritmo de resolución del sistema 𝑨𝑿 = 𝑩 mediante

factorización 𝑳𝑼 .

1. -Descomposición 𝑳𝑼: Del sistema de ecuaciones 𝐴𝑋 =

𝐵 se selecciona a la matriz 𝐴 y se descompone en sus
matrices 𝐿 y 𝑈 .

2. -Sustitución: Las matrices 𝐿 y 𝑈 se emplean para

calcular 𝑋 dado un cierto 𝐵. Este proceso a su vez se
divide en dos etapas:

2.1-Se resuelve la primera ecuación de [25], es decir,

𝐿𝐷 = 𝐵, cuya solución determina 𝐷 mediante sustitución
hacia adelante .

2.2-Posteriormente se sustituye 𝐷 en la segunda

ecuación de [25], es decir, 𝑈𝑋 = 𝐷 cuya solución
determina nuestra anhelada solución 𝑋 mediante
sustitución hacia atrás.

La gran valía del método de descomposición 𝐴 = 𝐿𝑈

consiste en que dicha factorización se realiza solo una vez,
y permite resolver el mismo sistema de ecuaciones
fácilmente para distintas 𝐵.
Notas Análisis Numérico Grupo 26 Unidad 3 Dr. Guillermo Alberto Sánchez Lozano

SUB DESCOMPOSICION(a,n)

DOFOR 𝑘 = 1, 𝑛 − 1 ¡recorre todos los elementos de la diagonal para calcular

¡los factores 𝑓

DOFOR 𝑖 = 𝑘 + 1, 𝑛 ¡Se ob enen todos los 𝑓 debajo del elemento 𝑎(𝑘, 𝑘)

𝑓𝑎𝑐𝑡𝑜𝑟 = 𝑎(𝑖, 𝑘)/𝑎(𝑘, 𝑘)

𝑎(𝑖, 𝑘) = 𝑓𝑎𝑐𝑡𝑜𝑟 ¡Factor 𝑓 de la matriz [𝐿] almacenado en matriz
( )
DOFOR 𝑗 = 𝑘 + 1, 𝑛 ¡recorre y ob ene todos los elementos 𝑎

¡del renglón 𝑖

𝑎(𝑖, 𝑗) = 𝑎(𝑖, 𝑗) − 𝑓𝑎𝑐𝑡𝑜𝑟 ∗ 𝑎(𝑘, 𝑗)

ENDDO

ENDDO

ENDDO

END DESCOMPOSICION

Tabla 5. Código para obtener la descomposición 𝐿𝑈 de un sistema de

ecuaciones. No considera pivoteo (ni escalamiento). Este código está
implementando en Fortran en la tabla 6.
Notas Análisis Numérico Grupo 26 Unidad 3 Dr. Guillermo Alberto Sánchez Lozano

Algoritmo numérico para resolución de un

sistema de ecuaciones mediante la
descomposición
En la tabla 6 se presenta un algoritmo desarrollado en
Fortran90 en el cual se implementa la descomposición 𝐿𝑈
a través de la eliminación de Gauss (es recomendable
comparar con el código de la tabla 4). Los puntos
relevantes del código son los siguientes:

 Los factores necesarios para la obtención de 𝑳 se

almacenan en la parte inferior de la matriz modificada .
Esto puede llevarse a cabo ya que durante la eliminación
de Gauss se colocan ceros, siendo estos espacios
aprovechables para almacenar las 𝑓 de 𝐿. Dicho
almacenamiento permite ahorrar espacio en memoria .
 El algoritmo usa un vector de orden (ordenamiento) 𝑶.
Esto acelera notablemente el algoritmo ya que evita
pivotear el vector 𝑶 (y no el renglón completo).
 Las ecuaciones no están escaladas, pero se usan los valores
escalados de los elementos para determinar si es necesario
el pivoteo.
 El término de la diagonal se verifica durante la fase de
pivoteo para detectar ocurrencias cercanas a cero, esto
para advertir de sistemas singulares . Si se devuelve un
valor de 𝑒𝑟 = −1, entonces se detectó una matriz singular
y se termina el cálculo. Lo anterior se controla mediante
un parámetro 𝑡𝑜𝑙 con un valor pequeño, para detectar
valores cercanos a cero.

El algoritmo descrito mencionado se implementa en el

código de la tabla 6.
Notas Análisis Numérico Grupo 26 Unidad 3 Dr. Guillermo Alberto Sánchez Lozano

PROGRAM ResolucionLU !Resolución mediante descomposición LU con pivoteo y escalamiento

IMPLICIT NONE
INTEGER,PARAMETER:: n=3 !orden del sistema de ecuaciones
INTEGER:: o(1:n),i,j,er
CHARACTER (LEN=3):: TAMANO,INDICE
REAL(KIND=8)::a(1:n,1:n),b(1:n),x(1:n)
!Declaración de la matriz [A] a la cual se obtendrá su descomposición [L][U]
a(1,1)=3.d0;a(1,2)=-0.1d0;a(1,3)=-0.2d0
a(2,1)=0.1d0;a(2,2)=7.d0;a(2,3)=-0.3d0
a(3,1)=0.3d0;a(3,2)=-0.2d0;a(3,3)=10.d0
!Declaración del vector {B}
b(1)=7.85d0;b(2)=-19.3d0;b(3)=71.4d0
!impresión de la matriz original
WRITE(*,*) ' Matriz [A] '
WRITE(TAMANO,'(I3)')n
DO i=1,n
WRITE(*,'('//TAMANO//'(f12.7,2x))') (a(i,j),j=1,n) !se usa DO IMPLICITO para escribir la matriz
ENDDO
!Obtención de la matriz que con ene cifradas [L] y [U]
CALL DESCOMPOSICION_LU(a,n,o,er)
!Obtención mediante sus tución adelante y hacia atrás del vector solución {X}
IF (er.NE.-1)THEN
CALL SUSTITUCION(a,o,n,b,x)
ENDIF
!Impresión de resultados para {X}
WRITE(*,*)
WRITE(*,*) ' Solución {X}'
DO i=1,n
WRITE(INDICE,'(I3)')i
WRITE(*,*) 'x('//INDICE//')=',x(i)
ENDDO
ENDPROGRAM
Tabla 6.1 Resolución de sistema por medio de descomposición LU con
pivoteo parcial en Fortran90. Programa principal donde se declara 𝐴 y
𝐵. Todo el código está basado en el pseudocódigo de la figura 10.2 de
Métodos Numéricos para Ingenieros 7ª edición, ejercicio 6.
Notas Análisis Numérico Grupo 26 Unidad 3 Dr. Guillermo Alberto Sánchez Lozano

!*****************************************************************************
!**************Subru na para descomposicion LU, en la que se llama a*********
!************************** subru na PIVOTEO*********************************
SUBROUTINE DESCOMPOSICION_LU(a,n,o,er)
IMPLICIT NONE
REAL (KIND=8),PARAMETER::tol=10.d0**(-8)
REAL(KIND=8)::a(1:n,1:n),s(1:n),factor
INTEGER:: o(1:n),i,j,k,n,er

DO i=1,n
o(i)=i !este vector lleva la cuenta el número de pivoteos.
s(i)=DABS(a(i,i))
DO j=2,n !Se determina cual es el elemento máximo de cada renglón en [A], es decir,s(i)
IF (DABS(a(i,j)).GT.s(i)) THEN
s(i)=DABS(a(i,j))
ENDIF
ENDDO
ENDDO

DO k=1,n-1,1
CALL PIVOTEO(a,o,s,n,k) !Se reordena la matriz [A] usando pivoteo parcial
IF (DABS(a(o(k),k)/s(o(k))).LT.tol) THEN !Determina si la matriz [A] es singular
er=-1
WRITE(*,*) a(o(k),k)/s(o(k))
WRITE(*,*) 'Matriz singular'
GOTO 100 !ESCAPA DEL CICLO Y CALCULOS SI [A] ES SINGULAR
ENDIF

DO i=k+1,n
factor=a(o(i),k)/a(o(k),k)
a(o(i),k)=factor
DO j=k+1,n
a(o(i),j)=a(o(i),j)-factor*a(o(k),j)
ENDDO
ENDDO
ENDDO
100 IF (DABS(a(o(n),n)/s(o(n))).LT.tol) THEN !Determina si la matriz [A] es singular
er=-1
WRITE(*,*) a(o(n),k)/s(o(n))
ENDIF
ENDSUBROUTINE DESCOMPOSICION_LU

Tabla 6.2 Subrutina DESCOMPOSICION_LU donde se almacena de

manera codificada 𝐿 y 𝑈 .
Notas Análisis Numérico Grupo 26 Unidad 3 Dr. Guillermo Alberto Sánchez Lozano
!*****************************************************************
!***********Subru na para realizar pivoteo parcial***************
SUBROUTINE PIVOTEO(a,o,s,n,k)
IMPLICIT NONE
REAL(KIND=8)::a(1:n,1:n),s(1:n),factor
REAL(KIND=8) a_mayor, tere
INTEGER:: o(1:n),i,j,k,n,p
p=k !Posición del pivote en la columna "k"
a_mayor=DABS(a(o(k),k)/s(o(k))) !elemento pivote normalizado

DO i=k+1,n
tere=DABS(a(o(i),k))/s(o(i))
IF ( tere.GT.a_mayor) THEN
a_mayor= tere
p=i
ENDIF
ENDDO
!En el siguiente proceso solo se intercambian las posiciones de o(p) y o(k), eso evita tener
!que intercambiar todos los elementos a(p,j) por a(k,j) (j=1,2,...,n)
tere=o(p)
o(p)=o(k)
o(k)= tere
ENDSUBROUTINE PIVOTEO
!*****************************************************************************
!Subru na que realiza sus tución hacia adelante y hacia atrás para encontrar
!**********resolver [L]{D}={B}, y posteriormente resolver [U]{X}={D}*********
SUBROUTINE SUSTITUCION(a,o,n,b,x)
IMPLICIT NONE
REAL(KIND=8)::a(1:n,1:n),b(1:n),x(1:n),sum
INTEGER:: o(1:n),i,j,k,n
!SUSTITUCIÓN HACIA ADELANTE
DO i=2,n
sum=b(o(i))
DO j=1,(i-1),1
sum=sum-a(o(i),j)*b(o(j))
ENDDO
b(o(i))=sum
ENDDO
!SUSTITUCIÓN HACIA ATRÁS
x(n)=b(o(n))/a(o(n),n)
DO i=n-1,1,-1
sum=0.d0
DO j=i+1,n
sum=sum+a(o(i),j)*x(j)
ENDDO
x(i)=(b(o(i))-sum)/a(o(i),i)
ENDDO
ENDSUBROUTINE SUSTITUCION

Tabla 6.3 Subrutinas PIVOTEO y SUSTITUCION. La primera realiza el pivoteo

parcial ahorrando espacio en memoria a través del vector 𝑜(𝑗) y qué evita pivotear
todos los elementos del 𝑎(𝑗, 𝑘) del correspondiente renglón. El segundo, realiza la
sustitución hacia adelante y después la sustitución hacia atrás para finalmente
obtener 𝑋 .
Notas Análisis Numérico Grupo 26 Unidad 3 Dr. Guillermo Alberto Sánchez Lozano

Ejercicio 6. Resuelva el siguiente sistema de ecuaciones,

3𝑥 − 0.1𝑥 − 0.2𝑥 = 7.85
0.1𝑥 + 7𝑥 − 0.3𝑥 = −19.3
0.3𝑥 − 0.2𝑥 + 10𝑥 = 71.4

Solución: Se aplica el código de la tabla 7 para obtener los

siguientes resultados,
Notas Análisis Numérico Grupo 26 Unidad 3 Dr. Guillermo Alberto Sánchez Lozano

Método de Jacobi
Hasta ahora, los métodos que hemos analizado se basan
únicamente en eliminación gaussiana, los cuales son
principalmente:

 Eliminación gaussiana con pivoteo parcial.

 Resolución mediante factorización de la matriz 𝐿𝑈 .

Los métodos anteriores consisten en realizar operaciones

entre renglones simplificando las ecuaciones en estos.

El método de Jacobi es el primer método iterativo para

encontrar la solución de sistemas de ecuaciones lineales.
Suponga que tiene el sistema de ecuaciones siguiente:
𝑎 𝑎 … 𝑎 𝑥 𝑏
⎛ 𝑎 𝑎 … 𝑎 ⎞ ⎛𝑥 ⎞ ⎛
⎜𝑏 ⎞
⎟
⎜
⎜ ⋮ ⎟
⎟⎜⎜ ⋮ ⎟
⎟=⎜ ⎟ [26]
⋮ ⋱ ⋮ ⎜⋮ ⎟
⎝𝑎 𝑎 … 𝑎 ⎠ ⎝𝑥 ⎠ ⎝𝑏 ⎠

Si los elementos de la diagonal principal son mayores que

la suma de los elementos de su respectivo renglón (en valor
absoluto), entonces se puede despejar una incógnita de
cada renglón. Tal condición se resume a la siguiente
ecuación

|𝑎 | > 𝑎 [27]
=
≠

Si un sistema de ecuaciones satisface la condición anterior,

se dice que la matriz es diagonalmente dominante .
Posteriormente a comprobar que se satisface la ecuación
[27], podemos obtener el siguiente sistema de ecuaciones
recursivas:

⎛ ⎞
𝑥 =⎜
⎜𝑏 −
⎜ 𝑎 𝑥 ⎟
⎟ 𝑎
⎟ [28]
=
⎝ ≠ ⎠
Notas Análisis Numérico Grupo 26 Unidad 3 Dr. Guillermo Alberto Sánchez Lozano

Las cuales se resuelven de manera iterativa y cuya

condición de convergencia se satisface cuando el error
relativo
𝑥 −𝑥 −
𝜀 = × 100% < 𝜀 [29]
𝑥

Donde 𝜀 es el error relativo asociado a la variable 𝑥 .

Ejercicio 7. Resuelva el siguiente sistema de ecuaciones

mediante el método de Jacobi:
5 −1 3 𝑥 8
1 −4 −1 𝑥 = −9 [𝐸1]
3 5 11 𝑥 13
Solución

En primer lugar, verificamos que la matriz es

diagonalmente dominante:
|𝑎 | = |5| > |−1| + |3| = 4
|𝑎 | = |−4| > |1| + |−1| = 2
|𝑎 | = |11| > |3| + |5| = 8

Como la matriz si satisface la condición [27], podemos

despejar de cada renglón su respectiva variable, lo cual nos
permite obtener las siguientes ecuaciones de recurrencia:
8 + 𝑥 − 3𝑥
𝑥 =
5
−9 − 𝑥 + 𝑥
𝑥 =
(−4)
13 − 3𝑥 − 5𝑥
𝑥 =
11
El ciclo de iteraciones se muestra a continua
Notas Análisis Numérico Grupo 26 Unidad 3 Dr. Guillermo Alberto Sánchez Lozano

𝑖 𝑥 , 𝑥 , 𝑥 , 𝑥 , 𝑥 , 𝑥 , 𝜀 , 𝜀 , 𝜀 ,
0 -10 12 30 -14 -7.75 -1.54545455 28.5714286 254.83871 2041.17647
1 -14 -7.75 -1.54545455 0.97727273 -0.86363636 8.52272727 1532.55814 797.368421 118.133333
2 0.97727273 -0.86363636 8.52272727 -3.68636364 0.36363636 1.30785124 126.510481 337.5 551.658768
3 -3.68636364 0.36363636 1.30785124 0.88801653 1.00144628 2.02190083 515.123313 63.6888797 35.3157572
4 0.88801653 1.00144628 2.02190083 0.58714876 1.96652893 0.484429 51.2421705 49.075436 317.378155
5 0.58714876 1.96652893 0.484429 1.70264838 2.27567994 0.12780992 65.51556 13.5849954 279.023014
6 1.70264838 2.27567994 0.12780992 1.97845004 2.64370962 -0.31694044 13.940289 13.9209569 140.326162
7 1.97845004 2.64370962 -0.31694044 2.31890619 2.82384762 -0.55944529 14.6817561 6.37916868 43.3473752
8 2.31890619 2.82384762 -0.55944529 2.5004367 2.96958787 -0.73417788 7.25995224 4.90776015 23.799762
9 2.5004367 2.96958787 -0.73417788 2.6344243 3.05865364 -0.84993177 5.08602971 2.91192743 13.6191979
10 2.6344243 3.05865364 -0.84993177 2.72168979 3.12108902 -0.92695828 3.20629811 2.00043551 8.30959901
11 2.72168979 3.12108902 -0.92695828 2.78039277 3.16216202 -0.97913768 2.11131985 1.29888984 5.32911712
12 2.78039277 3.16216202 -0.97913768 2.81991501 3.18988261 -1.01381713 1.40153999 0.86901613 3.42068109
13 2.81991501 3.18988261 -1.01381713 2.8462668 3.20843303 -1.03719619 0.92583694 0.578177 2.25406362
14 2.8462668 3.20843303 -1.03719619 2.86400432 3.22086575 -1.05281505 0.61932594 0.38600531 1.48353327
15 2.86400432 3.22086575 -1.05281505 2.87586218 3.22920484 -1.06330379 0.41232364 0.25823991 0.98642919
16 2.87586218 3.22920484 -1.06330379 2.88382324 3.23479149 -1.07032825 0.27605931 0.17270509 0.65629023
17 2.88382324 3.23479149 -1.07032825 2.88915525 3.23853787 -1.07503884 0.18455242 0.11568123 0.43817813
18 2.88915525 3.23853787 -1.07503884 2.89273088 3.24104852 -1.07819592 0.12360732 0.07746406 0.29281169
19 2.89273088 3.24104852 -1.07819592 2.89512726 3.2427317 -1.08031229 0.08277286 0.05190617 0.19590396
20 2.89512726 3.2427317 -1.08031229 2.89673372 3.24385989 -1.08173093 0.05545764 0.0347792 0.13114526
21 2.89673372 3.24385989 -1.08173093 2.89781054 3.24461616 -1.08268187 0.03715982 0.02330861 0.08783176
22 2.89781054 3.24461616 -1.08268187 2.89853236 3.2451231 -1.08331931 0.02490288 0.01562159 0.05884135
23 2.89853236 3.2451231 -1.08331931 2.89901621 3.24546292 -1.0837466 0.01669021 0.01047045 0.03942681
Notas Análisis Numérico Grupo 26 Unidad 3 Dr. Guillermo Alberto Sánchez Lozano

Métodos para calcular eigenvalores

eigenvectores
En ingeniería y en física la determinación de eigenvalores
juega un papel importante durante la solución de
problemas, entre los que destacan:

 Resolución de sistemas de ecuaciones diferenciales.

 Máximos y mínimos de sistemas de ecuaciones.
 Teoría de estabilidad.

El problema original consiste en encontrar las constantes 𝜆

que satisfagan

𝐴𝑣⃗ = 𝜆𝑣⃗ (𝑀 × ) [30]

Donde 𝐴 ∈ 𝑀 × , 𝑣⃗ ∈ 𝑀 × y 𝜆 es una constante. Una

forma de calcular las constantes que satisfacen [30] es
mediante la siguiente ecuación,

det(𝐴 − 𝜆𝐼) = 0 [31]

La cual al momento de expandirse proporciona el siguiente

polinomio
−
𝜆 +𝑏 − 𝜆 + ⋯+ 𝑏 𝜆 + 𝑏 = 0 [32]

Una vez calculados los eigenvalores se deben obtener los

eigenvectores asociados. Dicha tarea es bastante tortuosa
para sistemas grandes, no obstante, existen metodologías
para calcular el polinomio característico dado en [32] y en
ciertos casos hasta el mismísimo eigenvector asociado. Tal
metodología está compuesta por una serie de algoritmos
conocidos como Métodos de Krylov .
Notas Análisis Numérico Grupo 26 Unidad 3 Dr. Guillermo Alberto Sánchez Lozano

Método de las potencias

El método de las potencias es un método iterativo para
calcular el máximo eigenvalor de una matriz (y también
permite encontrar a la par su eigenvalor asociado)3. Para
ello se emplea la siguiente fórmula de recurrencia:

𝐴𝑣⃗( ) = 𝜆 + 𝑣⃗( + )
[33]

Es decir, el producto de 𝐴𝑣⃗( )

produce un vector, el cual se
normaliza mediante la factorización del elemento con la
magnitud más grande (el cual se desiganará con 𝜆 + y el
cual será nuestra aproximación al eigenvalor) y el vector
restante será nuestro eigenvector asociado 𝑣⃗( + )
.

El vector inicial 𝑣⃗( )

su elemento máximo debe ser a lo
más igual a 1.

Ejercicio 8. Obtenga el mayor eigenvalor de la siguiente

matriz junto con su eigenvector asociado:
−4 14 0
𝐴= −5 13 0 [𝐸1]
−1 0 2
Solución

A continuación, se comenzará el ciclo de iteraciones con el

vector inicial 𝑣⃗( ) = (0,1,1)

3
El método de las potencias asume que todos los eigenvalores son reales y existe un eigenvalor
dominante de tal manera que estos satisfacen la siguiente cadena de desigualdades:
|𝜆 | > |𝜆 | > |𝜆 | > ⋯ > |𝜆 |
Esto asume también que los eigenvectores forman una base independiente de vectores.
Notas Análisis Numérico Grupo 26 Unidad 3 Dr. Guillermo Alberto Sánchez Lozano

−4 14 0 0 14 1
𝐴𝑣⃗ =( )
−5 13 0 =1 = 14 ⎜13⁄14⎞
⁠ ⁠
13 ⎛ ⎟
−1 0 2 1 2 1
⎝ ⁄ ⎠ 7
9 1
−4 14 0 1 ⎛ 99 ⎞ ⎛ 11 ⎞
⎜
⎜ ⎟
⎟ ⎜
⎜ ⎟
⎛
𝐴𝑣⃗ = −5 13 0 ⎜13⁄14⎟ = ⎜
( ) ⎞
⎜ ⁠ 14 ⎟ ⁠⎟ = 9 ⎜
⎜ ⁠ 14 ⎟ ⁠⎟
⎟
−1 0 2 ⎝ 1⁄7 ⎠ ⎜ 5 ⎟ ⎜ 5 ⎟
− −
⎝ 7⎠ ⎝ 63⎠
1 7 1
−4 14 0 ⎜ ⎛ 11 ⎞ ⎛
⎜ 73 ⎞
⎟ ⎛ 73 ⎞
⎜ ⎟
⎟ ⎜ ⎟ ⎜
⎜ ⎟
( )
𝐴𝑣⃗ = −5 13 0 ⎜ ⎜ ⁠ 14 ⎟ ⁠⎟ = ⎜
⎜ ⁠ 14 ⎟ ⁠⎟ = 7 ⎜ ⎜ ⁠ 98 ⎟ ⁠⎟
⎟
−1 0 2 ⎜ 5 ⎟ ⎜ 73 ⎟ ⎜ 73 ⎟
− − −
⎝ 63⎠ ⎝ 63⎠ ⎝ 441⎠
45
1 ⎛ 1
⎛ 73 ⎞ ⎜ 7 ⎞ ⎟ ⎛ 51 ⎞
−4 14 0 ⎜ ⎟ ⎜
⎜ ⎟
⎟ 45 ⎜ ⎟
⎜ ⎟ 459 ⎜
( )
𝐴𝑣⃗ = −5 13 0 ⎜ ⎜ ⁠ 98 ⎟ ⁠⎟ = ⎜⎜ ⁠ ⁠⎟
⎟ = ⎜ ⁠ 70 ⎟ ⁠ ⎟
⎜ 73 ⎟ ⎜
⎜ 98 ⎟
⎟ 7 ⎜ 587 ⎟
⎜ ⎟
−1 0 2 ⎜ 587⎟
− −
⎝ 441⎠ − ⎝ 2835⎠
⎝ 441⎠
31
1 ⎛ ⎞ 1
51 ⎜ 5 ⎟ ⎛ 313 ⎞
−4 14 0 ⎛ ⎞ ⎜
⎜ ⎟ ⎜ 313 ⎟ ⎟ 31 ⎜ ⎟
𝐴𝑣⃗( ) = −5 13 0 ⎜ ⎜
⎜ ⁠ 70 ⎟
⁠ ⎟
⎟ = ⎜
⎜ ⁠ ⎟
⁠ ⎟ = ⎜
⎜
⎜ ⁠ 434 ⎟
⁠⎟
⎟
⎜ 587 ⎟ ⎜
⎜ 70 ⎟
⎟ 5 ⎜ 4009 ⎟
−1 0 2 ⎜ 4009⎟
− −
⎝ 2835⎠ − ⎝ 17577⎠
⎝ 2835⎠

De los resultados anteriores, determinamos en la quinta

iteración que

31 313 4009
𝜆≈ = 6.2 con 𝑣⃗ = 1, ,−
5 434 17577

El método de las potencias se puede emplear para calcular

a su vez el menor valor característico junto con su
eigenvector asociado mediante la siguiente ecuación de
recurrencia:
1
𝐴− 𝑣⃗( ) = 𝑣⃗( + )
[34]
𝜆 +

Donde 𝐴− designa la matriz inversa. El producto 𝐴− 𝑣⃗( )

produce un vector, el cual se normaliza mediante la

factorización del elemento más grande (en cuanto
Notas Análisis Numérico Grupo 26 Unidad 3 Dr. Guillermo Alberto Sánchez Lozano

magnitud, el cual se designa con 1/𝜆 + , y nos aproxima al

recíproco del menor eigenvalor), y el vector restante es
nuestro eigenvector asociado 𝑣⃗( + )
.

Ejercicio 9. Obtenga el menor eigenvalor de la matriz

siguiente
−4 14 0
𝐴= −5 13 0 [𝐸1]
−1 0 2
Solución

Para encontrar el menor eigenvalor asociado debemos

calcular primeramente la matriz inversa:
13 7
⎛ − 0⎞
⎜ 18 9 ⎟
⎜
⎜ 5 2 ⎟
⎟
𝐴− =⎜
⎜ − 0⎟
⎟
⎜
⎜ 18 9 ⎟
⎟
⎜13 7 1⎟
−
⎝36 18 2⎠
Es a partir de esta matriz que se obtiene el siguiente ciclo
de iteraciones con 𝑣⃗( ) = (0,1,1) :
13 7 7
⎛ − 0⎞ ⎛ − ⎞ 1
⎜ 18 9 ⎟ ⎜ 9⎟ ⎛ 2 ⎞
⎜ ⎟ 0 ⎜
⎜ 5 2 ⎟ ⎜ 2⎟⎟ 7⎜⎜ ⎟
⎟
𝐴− 𝑣⃗( ) = ⎜
⎜ − 0⎟
⎟ 1 = ⎜− ⎜ ⁠ ⎟ ⎟
⁠ =− ⎜ ⎜ ⁠ 7 ⁠⎟
⎟
⎜
⎜ 18 9 ⎟
⎟ ⎜
⎜ 9 ⎟
⎟ 9 ⎜ 1 ⎟
⎜13 1 ⎜ 1 ⎟
7 1⎟ ⎝ 7⎠
−
−
⎝36 18 2⎠ ⎝ 9 ⎠
13 7 1
⎛ − 0⎞ 1 ⎛ 1
⎜ 18 9 ⎟ ⎛ 2 ⎞ ⎜ 2⎞⎟ ⎛ 3⎞
⎜
⎜ 5 2 ⎟
⎟ ⎜ ⎟ ⎜
⎜ 3⎟⎟ 1 ⎜ ⎟
𝐴− 𝑣⃗( ) = ⎜ − ⎟ ⎜
0⎟ ⎜ ⁠ 7 ⎟ ⎟ ⎜ ⎟
⁠⎟ = ⎜⁠ ⁠⎟ = ⎜ ⎜ ⁠7⎟ ⁠ ⎟
⎜ ⎜ 2 ⎜5⎟
⎜
⎜
⎜ 18 9 ⎟
⎟ ⎜ 1⎟ ⎜ ⎜ 14⎟ ⎟ ⎟
⎜13 7 1⎟ ⎝− ⎠ ⎜ 5 ⎟ ⎝14⎠
− 7
⎝36 18 2⎠ ⎝28⎠
13 7 7
⎛ − 0⎞ 1 ⎛ 1
⎜ 18 9 ⎟ ⎛ 3 ⎞ ⎜ 18 ⎞ ⎟ ⎛ 23 ⎞
⎜
⎜ 5 2 ⎟
⎟ ⎜ ⎟ ⎜
⎜ 23 ⎟ ⎟ 7 ⎜ ⎟
𝐴− 𝑣⃗( ) = ⎜ − ⎟ ⎜
0⎟ ⎜ ⁠7⎟⎟ ⎜
⁠⎟ = ⎜⁠ ⁠⎟ = ⎜⎟ ⎜ ⁠ ⎟
49 ⁠ ⎟
⎜
⎜ ⎟ ⎜ 18 ⎜47⎟
⎜
⎜ 18 9 ⎟ ⎜5⎟ ⎜ ⎜126⎟ ⎟ ⎟
⎜13 7 1⎟ ⎝ ⎠ ⎜ 47 ⎟
− 14 ⎝49⎠
⎝36 18 2⎠ ⎝126⎠

Para la iteración 21
Notas Análisis Numérico Grupo 26 Unidad 3 Dr. Guillermo Alberto Sánchez Lozano

1468057
⎛39211644449⎞
39211644449 ⎜
⎜ ⎟
𝐴− 𝑥⃗( )
= ⎜
⎜ 7340027 ⎟⎟
⎟
78408608841 ⎜ ⎟
39211644449
⎝ 1 ⎠

De la ecuación [34] deducimos que

78408608841
𝜆= ≈ 1.99963
39211644449
Y el eigenvector asociado es
1468057
⎛
⎜39211644449⎞⎟ 0.00037438
⎜
⎜ 7340027 ⎟ ⎟ ≈ 0.00018719
⎜
⎜ ⎟
⎟
39211644449 1
⎝ 1 ⎠

Polinomio característico
Toda matriz 𝐴 ∈ 𝑀 × que tiene la ecuación característica
−
𝜆 +𝑏 − 𝜆 + ⋯+ 𝑏 𝜆 + 𝑏 = 0 [35]

Satisface su propia ecuación matricial asociada4:

𝐴 +𝑏 − 𝐴 −
+ ⋯ + 𝑏 𝐴 + 𝑏 𝐼 = 0⃗ [36]

La ecuación anterior puede reorganizarse de la siguiente

manera al ser multiplicada por un vector arbitrario 𝑦⃗ ∈
𝑀 × :

𝐴 𝑦⃗ + 𝑏 − 𝐴 −
𝑦⃗ + ⋯ + 𝑏 𝐴𝑦⃗ + 𝑏 𝑦⃗ = 0⃗

Que al reescribirse proporciona el siguiente sistema de

ecuaciones
−
𝑏 − 𝐴 𝑦⃗ + ⋯ + 𝑏 𝐴𝑦⃗ + 𝑏 𝑦⃗ = −𝐴 𝑦⃗ [37]

4
El hecho de que las ecuaciones [35] y [36] sean iguales, es debido a el Teorema de Cayley-Hamilton.
Notas Análisis Numérico Grupo 26 Unidad 3 Dr. Guillermo Alberto Sánchez Lozano

Lo anterior nos permite obtener un sistema de ecuaciones

lineales para los coeficientes 𝑏 (𝑖 = 0,1,2, … , 𝑛 − 1). De
esta forma, se puede determinar el polinomio característico
[35], el cual nos permitirá calcular los eigenvalores.

Ejercicio 10. Obtenga el polinomio característico de la

siguiente matriz
1 3 7 9 −1
⎛
⎜−3 2 2 7 8 ⎞⎟
𝐴=⎜
⎜
⎜ 4 5 1 1 −9⎟⎟
⎟ [𝐸1]
⎜1 −3 −4 3 0 ⎟
⎝8 6 3 2 2 ⎠
Solución

Debemos determinar los coeficientes del polinomio

característico de la matriz dada en [E1], es decir,

𝜆 +𝑏 𝜆 +𝑏 𝜆 +𝑏 𝜆 +𝑏 𝜆+𝑏 =0 [𝐸2]

Por el teorema de Cayley-Hamilton, sabemos que se

satisface la ecuación matricial asociada a [E2]

𝐴 + 𝑏 𝐴 + 𝑏 𝐴 + 𝑏 𝐴 + 𝑏 𝐴 + 𝑏 𝐼 = 0⃗ (𝑀 × ) [𝐸3]

Si multiplicamos la anterior ecuación por un vector

arbitrario 𝑦⃗ ∈ 𝑀 × , la ecuación anterior se puede
reorganizar así

𝐴 𝑦⃗ + 𝑏 𝐴 𝑦⃗ + 𝑏 𝐴 𝑦⃗ + 𝑏 𝐴 𝑦⃗ + 𝑏 𝐴𝑦⃗ + 𝑏 𝑦⃗ = 0⃗ (𝑀 × )

Al desarrollar productos
𝑏 + 𝑏 + 21𝑏 − 362𝑏 + 1667𝑏 − 19457 0
⎛
⎜ −3𝑏 + 70𝑏 + 60𝑏 + 5454𝑏 − 17147 ⎞⎟ ⎛0⎞
⎜ ⎟ ⎜
⎜ ⎟
⎟
⎜
⎜⁠ 4𝑏 − 78𝑏 + 173𝑏 − 4353𝑏 + 37810 ⎟
⁠ =⎜
⎟ 0 ⎟
⎜ ⎟ ⎜ ⎟
⎜ 𝑏 − 3𝑏 + 114𝑏 − 892𝑏 + 41 ⎟ ⎜0⎟
⎝ 8𝑏 + 20𝑏 + 388𝑏 − 1013𝑏 + 29191 ⎠ ⎝0⎠

Lo anterior puede reescribirse como el siguiente sistema

matricial:
Notas Análisis Numérico Grupo 26 Unidad 3 Dr. Guillermo Alberto Sánchez Lozano

1 1 21 −362 1667 𝑏 19457

⎛0 −3 70 60 ⎞ ⎛
5454 ⎟ ⎜ 𝑏 ⎞
⎟ ⎛ 17147 ⎞
⎜
⎜ ⎟ ⎜ ⎟ ⎜
⎜ ⎟
⎜
⎜0 4 −78 173 −4353 ⎟
⎟
⎜
⎜ ⎟
𝑏 ⎟=⎜
⎜−37810⎟⎟
⎟
⎜0 1 ⎟ ⎜ ⎟ ⎜
−3 114 −892 ⎜ 𝑏 ⎟ −41 ⎟
⎝0 8 20 388 −1013 ⎠ ⎝𝑏 ⎠ ⎝−29191⎠

Resolviendo este sistema se obtienen los siguientes

coeficientes para el polinomio característico:

𝜆 − 9𝜆 + 5𝜆 + 656𝜆 − 6671𝜆 + 29165 = 0