Tema 2

Sistemas de ecuaciones lineales.

Álgebra matricial

Los sistemas de ecuaciones lineales aparecen en numerosas áreas de la ciencia. Por ejemplo, se obtienen
al plantear las ecuaciones de equilibrio de un sistema de muelles acoplados, al resolver circuitos eléctricos
mediante las reglas de Kirchhoff, al ajustar por mínimos cuadrados un conjunto de datos numéricos, al
obtener autovalores de matrices a partir de las raíces de una ecuación secular,1 al resolver numéricamente
ecuaciones diferenciales parciales por el método de diferencias finitas,2 etc.
La resolución de sistemas de ecuaciones lineales está directamente relacionada con el álgebra de matri-
ces, por lo que también discutiremos en este tema varios aspectos del álgebra matricial. Nos limitaremos
a matrices cuadradas N × N con todos sus coeficientes reales. Los problemas más habituales que suelen
tratarse en el álgebra de matrices son:

Resolución de un sistema de ecuaciones lineales con el mismo número de ecuaciones que de incógni-
tas. En este caso, la matriz de referencia está formada por los coeficientes de las incógnitas. Si se intenta
resolver un sistema de 20 ecuaciones mediante la regla de Cramer, habría que realizar unas 5 × 109
multiplicaciones, lo cual resulta prohibitivo. Por ello se han desarrollado otros procedimientos más efi-
caces para resolver este tipo de problemas: los métodos directos, que proporcionan la respuesta exacta,
excepto errores de redondeo, y los métodos iterativos, cuya solución no es exacta, sino que se obtie-
ne con una precisión preestablecida, a partir de un conjunto de valores iniciales que se va mejorando
sucesivamente.

Operaciones con matrices (suma, resta y multiplicación).

Cálculo del determinante de una matriz.

Obtención de la inversa de una matriz.

La mayoría de matrices que encontramos en cálculos científicos y técnicos pueden agruparse en dos
grandes categorías: (i) de pequeño o mediano tamaño, pero densas (es decir, con pocos elementos nulos) y
(ii) de gran tamaño, pero dispersas o ralas (la mayoría de elementos son nulos). Los métodos que se tratarán
en este tema se pueden aplicar en el primer caso (cuando N . 100 y con pocos elementos nulos). Para
las matrices grandes y dispersas no siempre son apropiados los métodos discutidos en este tema (aunque
podríamos usarlos), por lo cual conviene emplear técnicas especiales, las cuales no se tratarán en este curso.
Para resolver sistemas de ecuaciones lineales estudiaremos en primer lugar los métodos iterativos; pos-
teriormente abordaremos los métodos directos basados en el procedimiento de eliminación hacia delante
y sustitución hacia atrás. Estos últimos métodos están relacionados con el álgebra de matrices; a partir de
ellos también pueden obtenerse procedimientos para calcular el determinante y la inversa de una matriz.
1 Se discutirá en el Tema 5.
2 Se discutirá en el Tema 6.

57
58 Rafael Garcia Molina, Departamento de Física – CIOyN, Universidad de Murcia

2.1. Métodos iterativos para resolver sistemas de ecuaciones lineales

Sea el sistema de N ecuaciones lineales con N incógnitas xi (i = 1, . . . , N ):

a11 x1 + a12 x2 + ... + a1N xN = b1

a21 x1 + a22 x2 + ... + a2N xN = b2
.. .. .. .. .. (2.1)
. . . . .
aN 1 x1 + aN 2 x2 + ... + aN N xN = bN .

Este sistema de ecuaciones también puede escribirse en notación matricial como A x = b, o abreviadamente
como (A | b), donde A es la matriz formada por los coeficientes aij de las incógnitas, x es el vector columna
que contiene las incógnitas xi , y b es el vector columna formado por los términos independientes bi . Es
decir:  
a11 a12 ... ... 
x1
 
b1

 a21 a22 ... ...   x2   b2 
 
A= .. .. . ..  , x =

.

, b =  ..  .
 
(2.2)
 .
 . 


 .
.

  . 
.. .. xN bN
. . aN N
Los métodos iterativos para resolver el sistema de ecuaciones lineales Ax = b son: (i) método de Jacobi,
(ii) método de Gauss-Seidel y (iii) método de sobrerrelajación. El método de Jacobi tiene carácter didáctico
y no se usa en la práctica, pues el de Gauss-Seidel es más eficiente. El método de sobrerrelajación es una
mejora sobre el de Gauss-Seidel.

2.1.1. Método de Jacobi

Este método es análogo al método iterativo de sustitución repetida empleado en el Tema 1 para obtener
los ceros de una función, pero ahora se aplica a un sistema de ecuaciones. Para comenzar el proceso iterativo
se despeja una incógnita de cada una de las ecuaciones del sistema (2.1)

xi = f (xj6=i ) ,
(0) (0) (0)
y se proponen unas soluciones aproximadas x1 , x2 , x3 , . . . para las incógnitas, con las cuales se iniciará
el proceso iterativo.
(0)
A partir de la aproximación inicial xi (con i = 1, 2, . . . , N ), se obtiene unos valores mejorados de
cada una de las incógnitas, de acuerdo con el siguiente proceso iterativo:
 
(k+1) 1  X (k) 
xi = bi − aij xj , i = 1, 2, ..., N (2.3)
aii
j6=i

que puede reescribirse como

 
i−1 N
(k+1) 1  X (k) (k)
X
xi = bi − aij xj − aij xj  , i = 1, 2, ..., N (2.4)
aii j=1 j=i+1

El superíndice k = 0, 1, 2, ... de cada incógnita indica la iteración en que se ha obtenido.

(k+1)
Nótese que en el método de Jacobi los valores xi obtenidos en la iteración k +1 no se emplean hasta
que no ha finalizado completamente dicha iteración k+1 (y, por tanto, comienza la siguiente iteración k+2).
Por este motivo, el método de Jacobi también se denomina método de los desplazamientos simultáneos,
ya que en la iteración k + 1 se emplean simultáneamente los valores de las incógnitas obtenidos en la
iteración precedente k.
Al sustituir los valores aproximados de xi (i = 1, 2, . . . , N ) en la parte derecha de la ec. (2.4) se
obtienen nuevos valores de xi , que (presumiblemente) se aproximarán más a las soluciones del sistema
de ecuaciones. El proceso iterativo se repite partiendo en cada ocasión de los valores de xi obtenidos
F ÍSICA COMPUTACIONAL , Sistemas de ecuaciones lineales. Álgebra matricial 59

en la iteración previa. Las iteraciones finalizan cuando haya poca diferencia (con una precisión elegida
previamente) entre los valores de xi obtenidos en iteraciones sucesivas.
(k+1)
El criterio de convergencia que se elija para finalizar el proceso iterativo puede ser que | xi −
(k)
xi |< i , donde i es un conjunto de valores pequeños establecidos previamente; generalmente se toma
i =  para todas las incógnitas. Otro criterio para finalizar las iteraciones es definir un valor M (k+1) =
(k+1) (k)
máx | xi − xi | y detener el proceso iterativo si M (k+1) < . Cuando convenga, en las expresiones
(k+1) (k) (k+1) (k) (k)
anteriores puede sustituirse el error absoluto | xi − xi | por el error relativo | (xi − xi )/xi |.
Una condición suficiente para que el proceso iterativo converja es que la matriz formada por los coefi-
cientes de las incógnitas sea diagonalmente dominante, lo cual implica que:
N
X
|aii | > |aij | para todas las filas (i = 1, 2, ..., N ) . (2.5)
j=1,j6=i

Es decir, cuando el valor absoluto de cada elemento de la diagonal es mayor que la suma de los valores
absolutos del resto de elementos en la misma fila se dice que la matriz de coeficientes es diagonalmente
dominante. Esto se debe a que, según la ecuación (2.3), el error del nuevo valor de xi será la suma de los
errores del resto de incógnitas multiplicados por los cocientes de coeficientes aij /aii , y si la suma de los
coeficientes es menor que la unidad, el error disminuirá conforme avance el proceso de iteración. La condi-
ción (2.5) es una condición suficiente (pero no necesaria) para que el proceso iterativo converja, aunque en
algunas situaciones también puede converger sin que se cumpla dicha condición. Evidentemente, la veloci-
dad de convergencia también está directamente relacionada con lo “dominantes” que son los elementos de
la diagonal con respecto al resto de elementos de la misma fila.
Cuando se verifica la condición (2.5) el método de Jacobi converge independientemente de los valores
(0)
iniciales xi de las incógnitas con los cuales comience la iteración.
Los principales inconvenientes del método de Jacobi son que: (i) la convergencia es de tipo lineal (por
lo que es lento) y (ii) puede que no converja si no se verifica la condición (2.5).

2.1.2. Método de Gauss-Seidel

(k+1)
El método de Jacobi no emplea los valores xi , recién calculados en la iteración k + 1, para despejar
(k+1)
los siguientes valores de xj<i en la misma iteración. Hasta la siguiente iteración no comienzan a usarse
(k+1)
los valores obtenidos en la iteración previa. Pero en la mayoría de situaciones, los valores nuevos xi
(k)
(obtenidos en la iteración k + 1) suelen ser mejores que los valores antiguos xi (obtenidos en la iteración
previa k) y conviene usarlos tan pronto estén disponibles.
El método de Gauss-Seidel usa los valores de xi recién calculados tan pronto como se necesitan, sin
esperar a completar cada iteración. De esta forma se ahorra espacio para el almacenamiento y se emplean
mejores valores de las xi tan pronto se conocen.
Para el sistema de N ecuaciones lineales (2.1), las incógnitas xi (con i = 1, 2, . . . , N ) se calculan
mediante el procedimiento iterativo
 
i−1 N
(k+1) 1 bi −
X (k+1)
X (k)
xi = aij xj − aij xj  , i = 1, 2, . . . , N (2.6)
aii j=1 j=i+1

(k+1) (k+1)
Nótese que en el cálculo de xi ya se usan los valores precedentes xj<i obtenidos en la misma itera-
ción.
El método de Gauss-Seidel también se denomina método de desplazamientos sucesivos, pues las
incógnitas antiguas se sustituyen (o desplazan) por las más recientes tan pronto se calculan.
El procedimiento para aplicar el método de Gauss-Seidel es análogo al de Jacobi. Las iteraciones se
detienen cuando se verifique el criterio de convergencia que impongamos a las incógnitas. Si la matriz
de coeficientes es diagonalmente dominante, el método de Gauss-Seidel converge para cualquier conjunto
de valores iniciales de las incógnitas, pero lo hace más rápidamente que el método de Jacobi. Por tanto, en
general siempre es preferible usar el método de Gauss-Seidel en lugar del método de Jacobi. Cuanto mejores
60 Rafael Garcia Molina, Departamento de Física – CIOyN, Universidad de Murcia

(más próximos a la solución) sean los valores de xi propuestos para iniciar la iteración, más rápidamente
convergerá el proceso.

2.1.3. Método de sobrerrelajación

El método de sobrerrelajación3 consiste en obtener una mejora x̄i para cada incógnita mediante un
(k) (k+1)
promedio ponderado de dos iteraciones consecutivas xi y xi obtenidas por el método de Gauss-
4
Seidel:
(k+1) (k+1) (k)
x̄i = ωxi + (1 − ω)xi , i = 1, 2, ...N (2.7)

Una elección adecuada del factor de sobrerrelajación ω permite acelerar notablemente el ritmo de conver-
gencia de las iteraciones. En principio, dado el radio espectral ρ de la matriz de Jacobi,5 puede determinarse
a priori el valor teórico óptimo de ω para aplicar el método de sobrerrelajación:

2
ω= p .
1 + 1 − ρ2

Pero el esfuerzo de computación para calcular ρ suele ser prohibitivo, por lo que, en general, no es posible
conocer por adelantado el valor de ω que maximiza el ritmo de convergencia del método de sobrerrelajación.
Si es posible estimar aproximadamente ρ, se pueden obtener valores razonables para el parámetro ω.
El método de sobrerrelajación no converge si ω está fuera del intervalo (0, 2). Aunque técnicamente
debería de usarse el término subrelajación (underrelaxation) cuando 0 < ω < 1, por comodidad se emplea
el término sobrerrelajación para cualquier valor de ω ∈ (0, 2). El valor ideal de ω suele estar comprendido
entre 1 y 2, aunque conviene tantear diferentes valores de ω para ver cuál es el que funciona mejor con el
problema particular que se está resolviendo. Si ω = 1, el método de sobrerrelajación equivale al método de
Gauss-Seidel.
El método de sobrerrelajación sucesiva se usa también cuando se discretizan ecuaciones diferenciales
en derivadas parciales, tal como se discutirá en el Tema 6.

2.1.4. Estrategias para mejorar los cálculos

En ocasiones puede que no sea conveniente despejar las incógnitas siempre en el mismo orden a lo largo
del proceso iterativo. Para evitar esto se puede sortear de un modo aleatorio (mediante números elegidos
al azar, tal como se explica en el Tema 7) la ecuación por la cual se comienza la secuencia de despejar
incógnitas. Si, además, no queremos seguir siempre la misma secuencia de incógnitas a resolver (asociadas
a las correspondientes ecuaciones), podemos sortear cada vez la ecuación a resolver. Lógicamente, estos
procedimientos hacen que los cálculos sean más lentos, pero pueden servir para evitar caer en resultados que
se verían afectados por errores acumulativos o por cualquier disposición “mal condicionada” o “viciada”
del sistema de ecuaciones.

Los métodos iterativos pueden aplicarse siempre (con las precauciones debidas sobre la naturaleza do-
minante de la diagonal de la matriz de coeficientes), pero cuando la matriz de coeficientes asociada al
sistema de ecuaciones a resolver es rala son recomendables frente a los métodos directos de eliminación
que se expondrán seguidamente, ya que suelen ser más rápidos.6
3 También denominado método de sobrerrelajación sucesiva (Successive Over Relaxation, SOR).
4 En algunos textos aparece la expresión:

(k+1) 1 (k+1) λ (k)

x̄i = x + x .
1+λ i 1+λ i
que equivale a la ec. (2.7) si tenemos en cuenta que ω = 1/(1 + λ).
5 El radio espectral ρ de una matriz A(N × N ) es el máximo de sus autovalores (en valor absoluto): ρ = máx |λ |, donde λ
i i
(i = 1, ..., N ) son los autovalores de la matriz.
6 Hay métodos de eliminación que están preparados para optimizar el cálculo cuando los ceros en los elementos de la matriz

aparecen formando patrones sencillos, como es el caso de matrices tridiagonales [DeVries: p.99].
F ÍSICA COMPUTACIONAL , Sistemas de ecuaciones lineales. Álgebra matricial 61

2.2. Método de Gauss: eliminación hacia delante y sustitución hacia

atrás
Un algoritmo para resolver exactamente (es decir, sin usar métodos aproximados) un sistema de pocas
ecuaciones lineales, como el correspondiente a las ecs. (2.1), consiste en combinar adecuadamente las
ecuaciones para eliminar incógnitas hasta que se obtenga una sola ecuación con una sola incógnita, la
cual se despeja inmediatamente. Durante la manipulación de las ecuaciones, se obtiene expresiones que
permitirán calcular las incógnitas restantes a partir de las que se han hallado previamente.
A continuación se explica un procedimiento sistemático para implementar en el ordenador un algoritmo
que permita resolver un sistema de ecuaciones.
Dadas N ecuaciones, el procedimiento es el siguiente: eliminamos x1 de las ecuaciones 2 hasta N ;
procedemos de forma análoga para eliminar x2 de las ecuaciones 3 hasta N , y así sucesivamente. Es decir,
vamos eliminando xi desde la ecuación (i + 1) hasta la ecuación N , ambas incluidas. Se procede de este
modo hasta que la última ecuación contenga sólo la variable xN , y así obtenemos su valor. Este método se
denomina eliminación hacia delante.
Veamos un caso práctico con el siguiente sistema de ecuaciones:

x1 + x2 + x3 = 6
−x1 + 2x2 = 3 (2.8)
2x1 + x3 = 5 .

En la primera etapa de este proceso eliminamos x1 de la segunda ecuación (i = 2). Para ello mul-
tiplicamos la primera ecuación (que es la primera fila, i = 1) por a21 /a11 = (−1)/1 y el resultado se
lo sustraemos a la segunda ecuación (que es la segunda fila, i = 2). Seguidamente eliminamos x1 de la
tercera ecuación (que es la tercera fila, i = 3) multiplicando la primera ecuación (primera fila, i = 1) por
a31 /a11 = 2/1 y el resultado se lo restamos a la tercera ecuación:

x1 + x2 + x3 = 6
3x2 + x3 = 9 (2.9)
−2x2 − x3 = −7 .

En la siguiente etapa eliminamos x2 desde la tercera ecuación en adelante. Para ello multiplicamos la
segunda ecuación (i = 2) por a32 /a22 = (−2)/3 y el resultado se lo substraemos a la tercera ecuación
(i = 3):
x1 + x2 + x3 = 6
3x2 + x3 = 9 (2.10)
1
3 x3 = 1 .
De la última ecuación anterior obtenemos xN , que en nuestro caso es x3 .
Una vez hemos obtenido el valor de la incógnita xN tras aplicar la eliminación hacia delante, a con-
tinuación procedemos en sentido inverso: sustituimos el valor de xN en la ecuación (N − 1), de donde
obtenemos el valor de xN −1 ; luego sustituimos los valores xN y xN −1 en la ecuación (N − 2) para obtener
xN −2 ; seguimos procediendo de esta manera hasta que llegamos a la primera ecuación, de donde obtenemos
x1 . Este método se denomina sustitución hacia atrás. En el ejemplo que hemos empleado, obtendríamos
la siguiente secuencia de resultados: x3 = 3, x2 = 2 y x1 = 1.
Este procedimiento de obtener las variables por eliminación hacia delante y sustitución hacia atrás7 se
denomina también método de eliminación gaussiana o, simplemente, método de Gauss. Los pasos a se-
guir en este método se pueden programar fácilmente en un ordenador. Para N ecuaciones con N incógnitas,
el número de operaciones requerido por el método de Gauss crece como N 3 . Afortunadamente, si el sistema
es ralo (muchos de los coeficientes son nulos), puede reducirse enormemente el tiempo de cálculo.

Implementación del método de eliminación gaussiana

Para implementar el método de eliminación hacia delante y sustitución hacia atrás se comienza elimi-
nando x1 de todas las ecuaciones, excepto la primera. Consideremos la ecuación i-ésima de un sistema de
7 En inglés se denomina forward elimination - backward substitution.
62 Rafael Garcia Molina, Departamento de Física – CIOyN, Universidad de Murcia

N ecuaciones:
ai1 x1 + ai2 x2 + · · · + aiN xN = bi . (2.11)
Si a esta ecuación i-ésima le restamos la primera ecuación multiplicada por ai1 /a11 , obtenemos:
     
ai1 ai1 ai1 ai1
ai2 − a12 x2 + ai3 − a13 x3 + · · · + aiN − a1N xN = bi − b1 , (2.12)
a11 a11 a11 a11

donde i = 2, ..., N . Después de haber eliminado la incógnita x1 desde la ecuación i = 2 hasta la i = N , se

repite el procedimiento para eliminar x2 desde la ecuación i = 3 hasta la i = N . Finalmente se llega a un
conjunto de ecuaciones triangular superior que se puede resolver fácilmente por sustitución hacia atrás.8

2.2.1. Eliminación gaussiana con pivote

El método de eliminación gaussiana da problemas cuando la diagonal principal contiene un término
nulo o muy pequeño, pues hay que dividir por aii en las expresiones que se obtienen (por ejemplo, a11 en
la ec. (2.12)). Por esta razón, antes de aplicar el algoritmo de eliminación de Gauss hay que reordenar las
ecuaciones con el fin de colocar en la diagonal de la matriz de coeficientes los términos con mayor valor
absoluto posible; con esto se garantiza que el correspondiente elemento de la diagonal no sea nulo (para
evitar problemas con las divisiones) y, además, reducir los posibles errores por redondeo que, de acuerdo
con el procedimiento de la resolución, se van acumulando paso tras paso.
Se denomina pivote al término por el cual se divide al despejar la incógnita de cada ecuación, y el
método se suele denominar eliminación gaussiana con pivote. Si se procede de este modo se evita que
aparezca una solución errónea por problemas de redondeo cuando el término de la diagonal es pequeño.
Incluso cuando todos los elementos de la matriz son de magnitud comparable, el procedimiento de elimina-
ción hacia delante puede dar lugar a términos pequeños en la diagonal principal del sistema de ecuaciones
resultante. Así pues, es esencial usar el método del pivote para todas las matrices (excepto, quizá, para las
más pequeñas). El precio a pagar por la reordenación de la ecuaciones es que hay que tener más cuidado
en la programación. Incluso cuando se usa el método del pivote, no se puede garantizar que aparezcan
problemas de redondeo al manejar matrices muy grandes.

2.3. Descomposición LU de una matriz

Una matriz A de dimensión N × N se puede escribir como el producto de una matriz triangular
inferior L (del inglés lower) multiplicada por otra matriz triangular superior U (del inglés upper):

A = LU .

Para simplificar la escritura supondremos una matriz de tamaño 4 × 4:9

 
a11 a12 a13 a14
 a21 a22 a23 a24 
A=  a31 a32 a33 a34  ,


a41 a42 a43 a44

que se descompondrá en el producto de las matrices L y U dadas por:10

   
1 0 0 0 u11 u12 u13 u14
 `21 1 0 0   0 u22 u23 u24 
L=  `31 `32 1 0  ,
 U = 0
 .
0 u33 u34 
`41 `42 `43 1 0 0 0 u44
8 Puede
consultarse los pasos detallados en [DeVries: p.119, Gaussian elimination].
9 En
clase empleamos una matriz 3 × 3, por brevedad.
10 Definimos la matriz L con los elementos de su diagonal principal igual a 1, tal como se hace en Numerical recipes (cuyas

subrutinas emplearemos) y la mayoría de textos, pero [DeVries: Cap. 3, Gaussian elimination] define la matriz U con los elementos
de la diagonal principal igual a 1.
F ÍSICA COMPUTACIONAL , Sistemas de ecuaciones lineales. Álgebra matricial 63

Queremos obtener los elementos de las matrices L y U , respectivamente, en función de los elementos
de la matriz original A. Para ello hemos de resolver el sistema de ecuaciones
    
a11 a12 a13 a14 1 0 0 0 u11 u12 u13 u14
 a21 a22 a23 a24   `21 1 0 0    0 u22 u23 u24  .
 
 a31 a32 a33 a34  =  `31 `32 1 0   0
  
0 u33 u34 
a41 a42 a43 a44 `41 `42 `43 1 0 0 0 u44

Procedemos por columnas y tenemos:

Columna 1
a11 = u11 → u11 = a11
a21 = `21 u11 → `21 = a21 /u11
a31 = `31 u11 → `31 = a31 /u11
a41 = `41 u11 → `41 = a41 /u11
de donde despejamos los valores de u11 , `21 , `31 y `41 (en este orden). Cuando a11 = 0 las divisiones por
u11 ocasionan un problema, que se resuelve si se aplica el método del pivote.
Columna 2

a12 = u12 → u12 = a12

a22 = `21 u12 + u22 → u22 = a22 − `21 u12
a32 = `31 u12 + `32 u22 → `32 = (a32 − `31 u12 )/u22
a42 = `41 u12 + `42 u22 → `42 = (a42 − `41 u12 )/u22

De aquí obtenemos sucesivamente los valores de u12 , u22 , `32 , `42 , etc. (en este orden). Vuelve a surgir el
mismo problema que aparecía anteriormente si al resolver la segunda ecuación resulta que u22 = 0, porque
entonces no podremos proseguir en la forma indicada.
Columna 3

a13 = u13 → u13 = a13

a23 = `21 u13 + u23 → u23 = a23 − `21 u13
a33 = `31 u13 + `32 u23 + u33 → u33 = a33 − `31 u13 − `32 u23
a43 = `41 u13 + `42 u23 + `43 u33 → `43 = (a43 − `41 u13 − `42 u23 )/u33

Fijémonos que para obtener el conjunto de elementos ui2 y `i2 (con el índice de columna igual a 2) no
se ha empleado ningún valor ai1 , es decir, de la columna primera de la matriz A. Lo mismo sucede con los
elementos de la columna j de A, que ya no se usan tras obtener los elementos uij y `ij . Por este motivo,
podemos almacenar en la columna j de la matriz original, los elementos uij y `ij a medida que se van
obteniendo, con el consiguiente ahorro de espacio.
Al finalizar la descomposición, en lugar de los elementos aij de la matriz original A, tendremos
 
u11 u12 u13 u14
 `21 u22 u23 u24 
  .
 `31 `32 u33 u34 
`41 `42 `43 u44

Esta forma de almacenar la matriz descompuesta tiene la ventaja de ahorrar memoria, pero la desventaja es
que se pierde la matriz original. No es necesario guardar los elementos de la diagonal principal de la matriz
L porque sabemos que todos ellos son la unidad.

Algoritmo de Crout
Un procedimiento que permite construir de forma eficiente las matrices L y U es el algoritmo de Crout.
En este método se construye una nueva matriz (que se almacena justo sobre la matriz A, para ahorrar espacio
en la memoria del ordenador) cuyos elementos se obtienen por columnas como sigue.
Para cada columna j = 1, 2, 3, ..., N se ejecutan los dos siguientes pasos:
64 Rafael Garcia Molina, Departamento de Física – CIOyN, Universidad de Murcia

Primero, para i = 1, 2, ..., j (elementos en y por encima de la diagonal principal) se obtiene:

i−1
X
uij = aij − `ik ukj .
k=1

Cuando i = 1 en la ecuación anterior, el sumatorio es nulo).

Segundo, para i = j + 1, j + 2, ..., N (elementos por debajo de la diagonal principal) se obtiene:

j−1
!
1 X
`ij = aij − `ik ukj .
ujj
k=1

Cuando j = 1 en la ecuación anterior, el sumatorio es nulo.

Una vez acabados estos dos pasos, se procede a calcular los elementos de la siguiente columna. En
cada etapa del cálculo sólo hace falta emplear los elementos calculados con anterioridad; además, cada
elemento aij sólo se utiliza una vez en los cálculos. Por este motivo, los valores de `ij y uij que vayamos
calculando pueden almacenarse en el lugar de la aij empleada para obtenerlos. Los elementos diagonales
no se almacenan porque ya sabemos que valen la unidad. Así pues, el algoritmo de Crout sustituye la matriz
A por una nueva matriz que contiene los elementos de las matrices L y U almacenados como sigue:
 
u11 u12 u13 u14
 `21 u22 u23 u24 
  .
 `31 `32 u33 u34 
`41 `42 `43 u44

Aplicación del método del pivote

Tal como se dijo anteriormente, para que el método sea estable ha de realizarse la operación de pivote.11
El elemento por el que hay que dividir al calcular `ij debe elegirse como el mayor (en valor absoluto) de los
elementos ujj , `ij (i = j + 1, ..., N ) que hay en la columna j; a este elemento se le denomina pivote. Esto
significa que el pivote pasa a ocupar la posición ujj en la diagonal. El pivote parcial realiza un intercambio
entre la fila del pivote y la fila j, con lo cual se ha modificado el problema original y se ha cambiado la
paridad de la ordenación de la matriz.
Conviene guardar un conjunto de números enteros (indx(i), en la subrutina ludcmp de Numerical
recipes, que usaremos más adelante) para indicarnos cómo se han permutado las filas de la matriz a lo largo
de la descomposición.
Hay subrutinas que escogen como pivote el elemento que hubiera sido el mayor si se hubieran escalado
todas las ecuaciones originales para que tengan su mayor coeficiente normalizado a la unidad; esto se
denomina pivote implícito. Aunque es necesario un mayor esfuerzo para almacenar los factores de escala
por el cual se hubieran multiplicado las filas (tal como lo hace la subrutina ludcmp de Numerical recipes).
Otra posibilidad consiste en el pivote total, donde se elige en cada momento como pivote el de mayor
valor absoluto de todos los elementos que quedan por reordenar. El pivote parcial es más sencillo que el
pivote total, pues no es necesario guardar las permutaciones del vector solución.
La operación del pivote es necesaria por si uno de los elementos ujj fuera nulo. Esta operación también
es muy importante para minimizar los errores de redondeo, ya que si cierto ujj ha resultado ser muy peque-
ño, se debe a una cancelación producida al calcularlo y, por lo tanto, habrá perdido dígitos significativos.
Por esta razón, durante la elaboración de cada columna se elige como pivote el elemento de mayor valor
absoluto.
Debido a los inevitables errores de redondeo, el método de eliminación de Gauss sin pivote es numéri-
camente inestable, incluso cuando no se encuentre ningún pivote nulo.
11 Es recomendable leer la discusión detallada sobre el método del pivote en [Numerical recipes: Cap. 2].
F ÍSICA COMPUTACIONAL , Sistemas de ecuaciones lineales. Álgebra matricial 65

Subrutina ludcmp
En este tema empleamos las subrutinas de Numerical recipes para resolver el sistema de ecuaciones
lineales Ax = b mediante el método de la descomposición LU de la matriz de coeficientes A. En la figura
2.1 se muestra la cabecera de la subrutina ludcmp, que se usa para descomponer una matriz A en el
producto LU .

Figura 2.1: Primeras líneas de la subrutina ludcmp empleada para la descomposición LU de una
matriz [Numerical recipes: p.38]. Es importante recordar que la dimensión física np ha de ser mayor
o igual que dimensión lógica n, la cual corresponde a la dimensión N de la matriz.

Para usar bien la subrutina ludcmp hay que distinguir entre la dimensión física np de la matriz y
su dimensión lógica n. La primera es la que se declara en la sentencia de la dimensión, mientras que la
segunda es el tamaño que tiene realmente la matriz. Es útil trabajar con la matriz de dimensión física (np ≥
n) porque en ocasiones se compila el programa para una dimensión física grande y luego se puede ejecutar
el programa (ya compilado) para diferentes problemas en los que la matriz tiene una dimensión lógica
menor. La figura 2.2 ilustra la diferencia entre las dimensiones física y lógica. Por ejemplo, si se usa la
subrutina ludcmp para alojar una matriz 4 × 4 en un espacio del ordenador con dimensión física 7 × 7,
hemos de escribir
parameter (n=4,np=7)
dimension a(np,np),indx(np)
Obsérvese que tanto la matriz A como indx se declaran con la dimensión física np.
A fin de no complicarnos la vida y para el propósito de este curso es suficiente usar los mismos valores
para las dimensiones física y lógica, a los cuales asignaremos el valor de la dimensión N de la matriz.

2.3.1. Resolución de un sistema de ecuaciones lineales

El siguiente sistema de ecuaciones
    
a11 a12 a13 a14 x1 b1
 a21 a22 a23 a24 
  x2
   b2 
 = 
 a31 a32 a33 a34   x3   b3 
a41 a42 a43 a44 x4 b4
puede escribirse en notación matricial como A x = b, el cual es equivalente a LU x = b.
Las incógnitas (contenidas en el vector columna x) se obtienen al resolver la ecuación U x = y, donde
el vector columna y es la solución de L y = b. Así pues, para resolver el sistema de ecuaciones A x = b se
realizan las dos etapas siguientes:
Primero se calcula y a partir de L y = b, mediante la sustitución de arriba hacia abajo en las ecuaciones:
y1 = b1
`21 y1 + y2 = b2
`31 y1 + `32 y2 + y3 = b3
`41 y1 + `42 y2 + `43 y3 + y4 = b4
66 Rafael Garcia Molina, Departamento de Física – CIOyN, Universidad de Murcia

Figura 2.2: Una matriz de dimensión lógica m × n se almacena en una matriz de dimensión física mp
× np. Los lugares marcados con aspas (×) contienen información irrelevante que puede sobrar de
otro uso previo de la matriz. Los números dentro de círculos muestran el orden exacto de los elementos
de la matriz en la memoria del ordenador (lo que no suele ser relevante para el programador). Nótese
que los elementos de la matriz lógica no ocupan lugares consecutivos en la memoria del ordenador.
Para que una subrutina localice un elemento (i, j) correctamente, hay que darle los valores de mp y
np, no sólo i y j [Numerical recipes: Fig. 2.0.1].

A continuación se obtiene x a partir de U x = y, tras sustituir desde abajo hacia arriba:

u11 x1 + u12 x2 + u13 x3 + u14 x4 = y1

u22 x2 + u23 x3 + u24 x4 = y2
u33 x3 + u34 x4 = y3
u44 x4 = y4

Subrutina lubksb

Para resolver un sistema de ecuaciones, nosotros emplearemos la subrutina lubksb, cuya cabecera
aparece en la figura 2.3; esta subrutina ha de emplearse justo después de usar la subrutina ludcmp comen-
tada anteriormente. Los términos independientes del sistema de ecuaciones se introducen a través del vector
columna b. Es importante recordar que en el programa principal (el que llama a las subrutinas ludcmp y
lubksb se han de declarar A, indx y b con dimensión física np.
En resumen, el procedimiento para resolver el sistema de ecuaciones lineales Ax = b consiste en usar
las subrutinas ludcmp y lubksb, en el orden y forma que se indica a continuación
call ludcmp(a,n,np,indx,d) para la descomposición LU de la matriz de coeficientes A.
call lubksb(a,n,np,indx,b) para resolver el sistema de ecuaciones mediante el método de sustitu-
ción hacia atrás.
El vector de incógnitas x se devuelve a través del vector b. Tanto la matriz original A como el vector
de términos independientes b se destruyen en este proceso, por lo que conviene guardarlos previamente con
otros nombres, si tenemos interés en conservarlos.
Si se desea resolver otro conjunto de ecuaciones con la misma matriz A, pero con un vector b de términos
independientes distinto, sólo hay que volver a repetir el proceso de llamar a la subrutina lubksb
call lubksb(a,n,np,indx,b) ,
donde a e indx son los valores que ha devuelto la subrutina ludcmp. Pero hay que asegurarse de suminis-
trar a la subrutina lubksb el vector b de términos independientes correspondiente al sistema que se desea
resolver (y no el vector b que se obtiene tras la llamada anterior a la subrutina lubksb, que no es otra cosa
que la solución x del sistema de ecuaciones resuelto con anterioridad) .
F ÍSICA COMPUTACIONAL , Sistemas de ecuaciones lineales. Álgebra matricial 67

Figura 2.3: Primeras líneas de la subrutina lubksb [Numerical recipes: p.39] empleada para obtener
las incógnitas de un sistema de ecuaciones mediante el método de sustitución hacia atrás, después de
haber realizado la descomposición LU de la matriz A de coeficientes mediante la subrutina ludcmp.

2.3.2. Determinante de una matriz

Una vez realizada la descomposición LU de una matriz A

A = LU , (2.13)

es muy fácil calcular el determinante de esta matriz, pues

det A = det L · det U . (2.14)

Como el determinante de L es la unidad y el de U es el producto de sus elementos diagonales, el

determinante de la matriz LU resulta finalmente:

N
Y
det A = uii , (2.15)
i=1

donde uii son los elementos de la diagonal principal de la matriz U , que coinciden con los elementos de la
diagonal principal aii de la matriz A que devuelve la subrutina ludcmp al aplicar el método de eliminación
gaussiana.
Hay que tener en cuenta que al realizar la descomposición LU mediante la subrutina ludcmp puede
haberse producido una reordenación de filas, y cada intercambio de filas implica un cambio de signo en el
determinante. La paridad p (par o impar) del total de estos cambios de signo debe determinarse a lo largo
del proceso de descomposición, por ello el determinante de la matriz original A se obtendrá multiplicando
por (−1)p el determinante resultante de la descomposición de la matriz original:

N
Y
det A = (−1)p uii . (2.16)
i=1

La subrutina ludcmp de Numerical recipes proporciona la variable d, que vale ±1 según sea la paridad
del número de permutaciones. Así pues, si empleamos esta variable d, el determinante de A será

N
Y
det A = d uii . (2.17)
i=1
68 Rafael Garcia Molina, Departamento de Física – CIOyN, Universidad de Murcia

Determinante de una matriz tridiagonal

Si A es una matriz tridiagonal de tamaño N × N:
 
b1 c1 0 0 ...
 a2 b2 c2 0 ... 
 
A= 0 a3 b3 c3 ... (2.18)
 

 .. .. 
 . . 
bN

se puede obtener su determinante mediante la aplicación sistemática de las siguientes relaciones:12

N = 1 → D 1 = b1 ,
N = 2 → D2 = b2 D1 − a2 c1 , (2.19)
N ≥ 3 → DN = bN DN −1 − aN cN −1 DN −2 .

2.3.3. Inversión de una matriz

Para calcular la inversa A−1 de la matriz A partimos de la identidad

A · A−1 = I . (2.20)

En lo que sigue, denominaremos cij a los elementos de la matriz inversa A−1 , por lo tanto
    
a11 a12 · · · a1N c11 c12 · · · c1N 1 0 ··· 0
 a21 a22 · · · a2N   c21 c22 · · · c2N   0 1 ··· 0 
..  =  .. .. . .  . (2.21)
    
 .. .. .. ..   .. .. .. ..
 . . . .  . . . .   . . . . 
aN 1 aN 2 ··· aN N cN 1 cN 2 ··· cN N 0 0 ··· 1

A continuación descomponemos las matrices A−1 e I en columnas, de forma que resulta

       
a11 a12 · · · a1N c11 c12 c1N
 a21 a22 · · · a2N   c21   c22   c2N 
..  ·  ..  t  ..  t · · · t  .. 
       
 .. .. ..
 . . . .   .   .   . 
aN 1 aN 2 ··· aN N cN 1 cN 2 cN N
     
1 0 0
 0   1   0 
=  ..  t   t ··· t   . (2.22)
     
.. ..
 .   .   . 
0 0 1

Podemos obtener la matriz inversa mediante la resolución de N sistemas de ecuaciones lineales tales
que el sistema de ecuaciones j-ésimo se obtiene a partir de los vectores columna (cij ) y (δij ), en la forma
     
a11 a12 · · · a1N c1j δ1j
 a21 a22 · · · a2N   c2j   δ2j 
..  ·  ..  =  ..  (para cada columna j) , (2.23)
     
 .. .. ..
 . . . .   .   . 
aN 1 aN 2 ··· aN N cN j δN j

donde δij es la función delta de Kronecker:

δij = 0 (i 6= j) ,

δij = 1 (i = j) .
12 [DeVries: p. 259].
F ÍSICA COMPUTACIONAL , Sistemas de ecuaciones lineales. Álgebra matricial 69

La matriz inversa de A estará formada por los vectores columna recién calculados, dispuestos uno a
continuación de otro:      
c11 c12 c1N
 c21   c22   c2N 
A−1 =  .  t  .  t · · · t  .  . (2.24)
     
 ..   ..   .. 
cN 1 cN 2 cN N
Mediante las subrutinas de descomposición LU (ludcmp) y de sustitución hacia atrás (lubksb) es
inmediato calcular la inversa de la matriz A, mediante el procedimiento de obtener sus columnas descrito
anteriormente.
A continuación se muestra una forma para suministrarle los vectores columna de la matriz identidad I a
la subrutina lubksb. Esta forma difiere ligeramente de la propuesta en [Numerical recipes: p.40], la cual
aparece en la figura 2.4.
1 * ai(np,np) = matriz inversa; y(np,np) = matriz identidad
2 dimension ai(np,np),y(np,np) ! hay que definir la dimension de las nuevas matrices
3 [...]
4 * defino los elementos de la matriz identidad
5 do i=1,n
6 do j=1,n
7 y(i,j)=0.0 ! todos los elementos son nulos...
8 end do
9 y(i,i)=1.0 ! ...excepto los de la diagonal principal
10 end do
11
12 call ludcmp(a,n,np,indx,d) ! descomposicion LU de la matriz A original
13
14 do j=1,n ! construyo los vectores columna de la matriz identidad "y"
15 do i=1,n
16 b(i)=y(i,j) ! el vector de terminos independientes es cada columna de "y"
17 end do
18 call lubksb(a,n,np,indx,b) ! solucion del sistema para cada vector columna
19 do i=1,n
20 ai(i,j)=b(i) ! almaceno los vectores columna de la matriz inversa
21 end do
22 end do
23 [...]

Hay que declarar la matriz y con dimensión (np,np). El resultado del programa devuelve en la matriz
ai la inversa de la matriz original a (que se habrá destruido en el proceso de descomposición LU ).

Figura 2.4: Procedimiento que emplea [Numerical recipes: p.40] para calcular la matriz inversa ai
de la matriz a.

2.4. Sistemas de ecuaciones no lineales

Resolver un sistema de ecuaciones no lineales es un problema bastante más complicado que en el caso
de ecuaciones lineales. De hecho, algunos sistemas de ecuaciones no lineales carecen de soluciones reales.
70 Rafael Garcia Molina, Departamento de Física – CIOyN, Universidad de Murcia

Para resolver sistemas de ecuaciones no lineales se puede emplear un algoritmo iterativo, pero no está
garantizado que proporcione la solución. Por lo tanto, los algoritmos iterativos han de usarse con las debidas
precauciones. También se puede intentar la resolución de una versión linealizada de las ecuaciones.
Si el sistema consta de más de dos ecuaciones no lineales, es muy difícil encontrar una forma conver-
gente para poder resolver las ecuaciones.
Un criterio de convergencia (condición suficiente únicamente) es el siguiente:
El conjunto de ecuaciones

x = f (x, y, z, . . . ), y = g(x, y, z, . . . ), z = h(x, y, z, . . . ), ...

convergerá si, en un intervalo alrededor de la raíz, se cumple que

     
 ∂f   ∂f   ∂f 
 ∂x  +  ∂y  +  ∂z  + . . . < 1,
     
 ∂g   ∂g   ∂g 
 ∂x  +  ∂y  +  ∂z  + . . . < 1,
 ∂h   ∂h   ∂h 
  +   +  + ... < 1,
∂x ∂y ∂z
... ... < 1.

El tener que calcular todas las derivadas parciales y conocer dónde se halla aproximadamente cada raíz
conlleva un gran esfuerzo. Incluso, encontrar valores adecuados para comenzar con las iteraciones también
supone una gran dificultad.
2.5. Proyecto: Reglas de Kirchhoff
Un problema típico de Física en el que se resuelve un sistema de ecuaciones lineales surge cuando se
aplican las reglas de Kirchhoff para estudiar circuitos de corriente continua.
A continuación vamos a aplicar estas reglas a un caso concreto. En el circuito eléctrico representado en
la parte izquierda de la figura 2.5 se conocen los valores de las resistencias Rj (j = 1, ..., 12) y el voltaje
V aplicado. Se desea calcular las diferencias de potencial entre los extremos de cada resistencia. Como
la disposición de las resistencias no permite su agrupación en serie ni en paralelo para reemplazarlas por
su resistencia equivalente (salvo en unos pocos casos),13 aplicaremos las reglas de Kirchhoff, tal como se
explicará seguidamente.
De acuerdo con la ley de Ohm, obtendremos la diferencia de potencial Vj entre los extremos de una
resistencia Rj si conocemos la intensidad Ij que circula por ella, pues Vj = Ij Rj . Aplicaremos las reglas
de Kirchhoff para calcular las intensidades Ij (j = 1, ..., 13) que circulan por las resistencias del circuito,
de acuerdo con la nomenclatura empleada en la parte derecha de la figura 2.5.

2.5.1. Reglas de Kirchhoff

Regla de los nudos. Un nudo es un punto de ramificación de un circuito donde puede dividirse la co-
rriente. La suma de las intensidades que entran en un nudo ha de ser igual a la suma de las intensidades que
salen del mismo: X
Ij = 0 (en cada nudo).
j

Esta regla es una consecuencia de la conservación de la carga. Puesto que no hay ninguna fuente ni
sumidero de cargas en el circuito, y las cargas no se acumulan en un nudo, el número de cargas que llegan
a un nudo en un cierto intervalo de tiempo debe ser igual al número de cargas que salen de dicho nudo en
el mismo intervalo temporal.
La regla de los nudos se aplica a todos los nudos menos uno (pues la información que proporcionaría
sería redundante).14
13 Algunas resistencias (tales como las R y R o las R y R ) sí que pueden asociarse en serie, pero las mantendremos sin asociar
3 7 6 10
para calcular la diferencia de potencial entre sus extremos.
14 Si el circuito tiene n nudos y la regla se satisface para (n − 1) nudos, automáticamente se satisface para el nudo restante.
F ÍSICA COMPUTACIONAL , Sistemas de ecuaciones lineales. Álgebra matricial 71

Figura 2.5: (Izq.) Conjunto de resistencias Rj (i = 1, . . . , 12) a la que se aplica una diferencia de
potencial V entre los puntos a e i. (Der.) Propuesta de sentidos de las intensidades (líneas gruesas
en color azul) y de recorrido de las mallas (líneas delgadas en color rosa) para resolver el circuito
mediante las reglas de Kirchhoff.

Regla de las mallas. Una malla (o bucle) es una parte del circuito que puede recorrerse en una tra-
yectoria cerrada, sin pasar dos (o más) veces por el mismo lugar. La suma algebraica de las variaciones de
potencial ∆Vj a lo largo de cualquier malla del circuito debe de ser cero:
X
∆Vj = 0(en cada malla).
j

P P
Esto implica que j Ij Rj = j femj , donde Ij Rj son las caídas de potencial debido a los elementos re-
sistivos (elementos pasivos), mientras que femj corresponde a las fuentes de potencial (elementos activos).
Esta regla se basa en el principio de conservación de la energía, ya que la variación neta de energía
de una carga tras recorrer un camino cerrado ha de ser nula. Si tenemos una carga q en un punto donde el
potencial es V , la energía potencial de la carga vale qV . Cuando la carga recorre un bucle en un circuito,
pierde o gana energía al atravesar resistencias, baterías u otros elementos, pero cuando vuelve a su punto
de partida, su energía ha de ser nuevamente qV . Es decir, el cambio neto en el potencial ha de ser cero al
recorrer un bucle del circuito.
P P
Al aplicar esta regla (en la forma j Ij Rj = j femj ), hay que tener en cuenta lo siguiente. Una caída
de potencial a través de una resistencia es positiva (o negativa) según se recorra el circuito en el sentido de
la corriente (o en sentido opuesto). Cuando se pasa a través de una fem, se toma la diferencia de potencial
como positiva (o negativa) según se atraviese desde la terminal negativa a la positiva (o al revés).
La regla de las mallas se aplica a tantas mallas del circuito como sea necesario para que cada elemento
dle circuito forme parte de un tramo recorrido al menos una vez y sin que se recorra dos veces la misma
malla.

2.5.2. Aplicación de las reglas de Kirchhoff

A continuación aplicamos las reglas de Kirchhoff, de acuerdo con la nomenclatura que se muestra en la
parte derecha de la figura 2.5. Los nudos están indicados en la figura por las letras minúsculas.
Regla de los nudos. Las intensidades que circulan por cada resistencia se han representado en la figura
mediante flechas de color azul; no hemos de preocuparnos por el sentido de las intensidades propuestas,
pues el resultado que obtengamos finalmente confirmará (o cambiará) el signo de cada intensidad. Para
aplicar la regla de los nudos hemos elegido el convenio (arbitrario) de que las intensidades que entran a los
nudos son negativas y las que salen son positivas). A continuación se indican los nudos donde se aplica la
72 Rafael Garcia Molina, Departamento de Física – CIOyN, Universidad de Murcia

primera regla de Kirchhoff, así como la suma de intensidades en cada nudo:

nudo a: I1 + I2 − I13 = 0
nudo b: −I1 + I3 + I4 = 0
nudo c: −I2 − I5 + I6 = 0
nudo d: −I3 + I7 = 0
nudo e: −I4 + I5 − I8 + I9 = 0
nudo f: −I6 + I10 = 0
nudo g: −I7 + I8 + I11 = 0
nudo h: −I9 − I10 + I12 = 0
Regla de las mallas. En la figura se ha representado mediante flechas de color rosa el sentido en el
cual se recorre cada una de las mallas a las que se aplica la segunda regla de Kirchhoff. Al aplicar la regla
de las mallas, hay que tener en cuenta el sentido en que se recorren. Habría que ponerle un signo menos a
las diferencias de potencial entre los extremos de una resistencia cuando la intensidad que la recorre tenga
sentido opuesto al del recorrido de la malla. También habría que poner un signo negativo a las fuentes de
potencial que estén dispuestas para suministrar corriente en sentido contrario al propuesto para recorrer la
malla. A continuación se indican las mallas (identificadas por la secuencia de letras en que se recorren, en
sentido horario, en este caso) y la suma de variaciones de potencial en cada caso:
malla aceba: −I1 R1 + I2 R2 − I4 R4 − I5 R5 = 0
malla begdb: −I3 R3 + I4 R4 − I7 R7 − I8 R8 = 0
malla cfhec: I5 R5 + I6 R6 − I9 R9 + I10 R10 = 0
malla ehige: I8 R8 + I9 R9 − I11 R11 + I12 R12 = 0
malla iV acfhi: I2 R2 + I6 R6 + I10 R10 + I12 R12 = V
En total, tenemos un sistema de 13 ecuaciones, con 13 incógnitas (Ij , j = 1, ..., 13):
+I1 +I2 · · · · · · · · · · −I13 = 0
−I1 · +I3 +I4 · · · · · · · · · = 0
· −I2 · · −I5 +I6 · · · · · · · = 0
· · −I3 · · · +I7 · · · · · · = 0
· · · −I4 +I5 · · −I8 +I9 · · · · = 0
· · · · · −I6 · · · +I10 · · · = 0
· · · · · · −I7 +I8 · · +I11 · · = 0
· · · · · · · · −I9 −I10 · +I12 · = 0
−I1 R1 +I2 R2 · −I4 R4 −I5 R5 · · · · · · · · = 0
· · −I3 R3 +I4 R4 · · −I7 R7 −I8 R8 · · · · · = 0
· · · · +I5 R5 +I6 R6 · · −I9 R9 +I10 R10 · · · = 0
· · · · · · · +I8 R8 +I9 R9 · −I11 R11 +I12 R12 · = 0
· +I2 R2 · · · +I6 R6 · · · +I10 R10 · +I12 R12 · = V

En notación matricial, este sistema se escribe:


−1
   
1 1 0 0 0 0 0 0 0 0 0 0 
 I1 0

 −1 0 1 1 0 0 0 0 0 0 0 0 0 
 I2  
  0 


 0 −1 0 0 −1 1 0 0 0 0 0 0 0 
 I3  
  0 


 0 0 −1 0 0 0 1 0 0 0 0 0 0 
 I4  
  0 


 0 0 0 −1 1 0 0 −1 1 0 0 0 0 
 I5  
  0 


 0 0 0 0 0 −1 0 0 0 1 0 0 0 
 I6  
  0 


 0 0 0 0 0 0 −1 1 0 0 1 0 0 
 I7 =
  0 


 0 0 0 0 0 0 0 0 −1 −1 0 1 0 
 I8  
  0 


 −R1 R2 0 −R4 −R5 0 0 0 0 0 0 0 0 
 I9  
  0 


 0 0 −R3 R4 0 0 −R7 −R8 0 0 0 0 0 
 I10  
  0 


 0 0 0 0 R5 R6 0 0 −R9 R10 0 0 0 
 I11  
  0 

 0 0 0 0 0 0 0 R8 R9 0 −R11 R12 0 
 I12   0 
0 R2 0 0 0 R6 0 0 0 R10 0 R12 0  I13 V

Mediante el programa 2-09LU-rgm.f se resuelve ese sistema de ecuaciones por el método de des-
composición LU . Si empleamos los valores
R1 = R3 = R4 = R9 = R10 = R12 = 1 Ω , R2 = R5 = R6 = R7 = R8 = R11 = 2 Ω , V = 10 V ,
el resultado que se obtiene es:
F ÍSICA COMPUTACIONAL , Sistemas de ecuaciones lineales. Álgebra matricial 73

I1 (A) I2 (A) I3 (A) I4 (A) I5 (A) I6 (A) I7 (A) I8 (A) I9 (A) I10 (A) I11 (A) I12 (A) I13 (A)
3.056 1.806 1.111 1.944 −0.694 1.111 1.111 −0.694 1.944 1.111 1.806 3.056 4.861

A partir de la relación Vj = Ij Rj , obtenemos las siguientes diferencias de potencial en cada resistencia:

V1 (V) V2 (V) V3 (V) V4 (V) V5 (V) V6 (V) V7 (V) V8 (V) V9 (V) V10 (V) V11 (V) V12 (V)
3.056 3.611 1.111 1.944 −1.389 2.222 2.222 −1.389 1.944 1.111 3.611 3.056

Bibliografía
G. J. Borse, Programación en FORTRAN77, con aplicaciones de cálculo numérico en ciencias e inge-
niería (Anaya, Madrid, 1989). Cap. 12.
P. DeVries, A First Course in Computational Physics (Wiley, New York, 1994). Cap. 3 (principalmente
las secciones Tridiagonal Linear Systems y Gaussian elimination).
A. L. Garcia, Numerical methods for physics, 2nd ed. (Prentice-Hall, Upper Saddle River, NJ, 2000).
Cap. 4.
C. F. Gerald, Applied Numerical Analysis (Addison-Wesley, Reading, Mass., 1980). Cap. 2.
R. Guardiola, E. Higón y J. Ros, Mètodes numèrics per a la física (Universitat de València, Valencia,
1997). Cap. 3.

W. H. Press, S. A. Teukolsky, W. T. Vetterling y B. P. Flannery, Numerical Recipes, in Fortran. The Art

of Scientific Computing, 2nd ed. (Cambridge University Press, Cambridge, 1992). Cap. 2.
L. Vázquez, S. Jiménez, C. Aguirre y P. J. Pascual, Métodos numéricos para la Física y la Ingeniería
(McGraw-Hill, Madrid, 2009). §§6.1-6.5 y 2.8.
74 Rafael Garcia Molina, Departamento de Física – CIOyN, Universidad de Murcia

Resumen del Tema 2: Sistemas de ecuaciones lineales. Álgebra matri-

cial
Métodos iterativos
(0)
Se propone una aproximación inicial xi al conjunto de i = 1, ..., N incógnitas que aparecen en el
sistema de N ecuaciones definido por la relación matricial
     
a11 a12 · · · a1N x1 b1
 a21 a22 · · · a2N   x2   b2 
..  ·  .. =  .
     
 .. .. .. ..
 . . . .   .   . 
aN 1 aN 2 ··· aN N xN bN
(k)
Se obtiene la secuencia xi mediante un proceso iterativo, de manera que los valores obtenidos en cada
iteración k converjan hacia la solución del sistema.
A continuación aparece un resumen de los principales métodos iterativos.

Método de Jacobi
 
i−1 N
(k+1) 1  X (k)
X (k) 
xi = bi − aij xj − aij xj , i = 1, 2, ..., N
aii j=1 j=i+1

Método de Gauss-Seidel
 
i−1 N
(k+1) 1  X (k+1)
X (k) 
xi = bi − aij xj − aij xj , i = 1, 2, ..., N
aii j=1 j=i+1

Una condición suficiente para que el proceso iterativo converja es que la matriz de coeficientes sea diago-
nalmente dominante:
N
X
|aii | > |aij | para todas las filas (i = 1, 2, ..., N ) .
j=1,j6=i

Cuando se verifica esta condición, los dos métodos anteriores convergen para cualquier conjunto de valores
(0)
xi con el que se inicia la iteración. El método de Gauss-Seidel converge más rápido que el de Jacobi.

Sobrerrelajación
El método de sobrerrelajación (Successive Over Relaxation, SOR) consiste en obtener una mejora x̄i
para cada incógnita mediante un promedio ponderado de dos iteraciones consecutivas obtenidas por el
método de Gauss-Seidel:
(k+1) (k+1) (k)
x̄i = ω xi + (1 − ω) xi , k = 0, 1, 2, . . .

El factor de sobrerrelajación ω está comprendido entre 0 y 2. Si ω = 1, tenemos el método de Gauss-Seidel.

El valor óptimo de ω para aplicar el método de sobrerrelajación sucesiva es
2
ω= p ,
1 + 1 − ρ2

donde ρ es el radio espectral de la matriz de coeficientes. Pero el esfuerzo para calcular ρ suele ser prohi-
bitivo, por lo que, en general, no es posible conocer por adelantado el valor de ω que maximiza el ritmo de
F ÍSICA COMPUTACIONAL , Sistemas de ecuaciones lineales. Álgebra matricial 75

convergencia. Por ello, conviene tantear diferentes valores de ω para ver cuál es el que funciona mejor con
el problema particular que se está resolviendo.

Álgebra matricial
Descomposición LU de una matriz
Una matriz A de dimensión N × N se puede escribir como el producto de una matriz triangular inferior
L (del inglés lower) multiplicada por otra matriz triangular superior U (del inglés upper):

A = LU .

Para simplificar la escritura supondremos una matriz de tamaño 4 × 4:

 
a11 a12 a13 a14
 a21 a22 a23 a24 
A=  a31 a32 a33 a34  ,


a41 a42 a43 a44

que se descompondrá en el producto de las matrices L y U :

   
1 0 0 0 u11 u12 u13 u14
 `21 1 0 0   0 u22 u23 u24 
L=  `31 `32 1 0  ,
 U =  .
 0 0 u33 u34 
`41 `42 `43 1 0 0 0 u44

Algoritmo de Crout
Se emplea para calcular las matrices L y U (que se almacenan justo sobre la matriz A original). Los
elementos de la nueva matriz se obtienen por columnas mediante el procedimiento que se detalla a conti-
nuación.
Para cada columna j = 1, 2, 3, ..., N se ejecutan los dos pasos siguientes.
Primero, para i = 1, 2, ..., j (elementos en y por encima de la diagonal principal) se obtiene:
i−1
X
uij = aij − `ik ukj .
k=1

(Cuando i = 1 en la ecuación anterior, el sumatorio es nulo).

Segundo, para i = j + 1, j + 2, ..., N (elementos por debajo de la diagonal principal) se obtiene:
j−1
!
1 X
`ij = aij − `ik ukj .
ujj
k=1

(Cuando j = 1 en la ecuación anterior, el sumatorio es nulo).

Una vez acabados estos dos pasos, se procede a calcular los elementos de la siguiente columna. En
cada etapa del cálculo sólo hace falta emplear los elementos calculados con anterioridad; además, cada
elemento aij sólo se utiliza una vez en los cálculos. Por este motivo, los valores de `ij y uij que vayamos
calculando pueden almacenarse en el lugar de la aij empleada para obtenerlos. Los elementos diagonales
no se almacenan porque ya sabemos que valen la unidad. Así pues, el algoritmo de Crout sustituye la matriz
A por una nueva matriz que contiene los elementos de las matrices L y U almacenados como sigue:
 
u11 u12 u13 u14
 `21 u22 u23 u24 
 `31 `32 u33 u34  .
 

`41 `42 `43 u44

76 Rafael Garcia Molina, Departamento de Física – CIOyN, Universidad de Murcia

De esta forma se ahorra espacio, pero se pierde la matriz A original; si interesa conservarla, antes de emplear
el algoritmo de Crout hay que guardarla en otra matriz.
Para que el método sea estable ha de realizarse la operación de pivote: el elemento de la diagonal ujj
por el que se divide al determinar `ij debe elegirse como el mayor de los elementos `ij (i = j + 1, ..., N )
que hay en la columna j; a este elemento se le denomina pivote. Esto significa que el pivote pasa a ocupar la
posición ujj en la diagonal. Para ello habrá que realizar los intercambios de filas que sean necesarios (pivote
parcial), con lo cual se ha modificado el problema original y se ha cambiado la paridad de la ordenación de
la matriz. Conviene guardar un conjunto de números enteros que indique las permutaciones de filas que se
ha efectuado al proceder a la descomposición LU de la matriz A.

Sistema de ecuaciones lineales

   
a11 a12 a13 a14 x1 b1
 a21 a22 a23 a24   x2   b2 
  = 
 a31 a32 a33 a34   x3   b3 
a41 a42 a43 a44 x4 b4
Este sistema se escribe en notación abreviada como Ax = b, el cual equivale a LU x = b. El vector de
incógnitas x se obtiene al resolver la ecuación U x = y, donde el vector columna y es la solución de
Ly = b.
Así pues, para resolver el sistema de ecuaciones Ax = b se procede en dos etapas. Primero se calcula y
a partir de la relación Ly = b. Para ello se procede a sustituir de arriba hacia abajo en:
y1 = b1
`21 y1 + y2 = b2
`31 y1 + `32 y2 + y3 = b3
`41 y1 + `42 y2 + `43 y3 + y4 = b4
A continuación se obtiene x del sistema sistema U x = y. Ahora se procede a sustituir de abajo hacia arriba:
u11 x1 + u12 x2 + u13 x3 + u14 x4 = y1
u22 x2 + u23 x3 + u24 x4 = y2
u33 x3 + u34 x4 = y3
u44 x4 = y4

Determinante de una matriz

Una vez obtenida la descomposición LU de una matriz A, se calcula el determinante de A como sigue:
N
Y
det A = (−1)p uii ,
i=1

donde el término (−1)p tiene en cuenta la paridad (p par o impar) de los intercambios de filas que se ha
efectuado en la matriz A para descomponerla en el producto LU .
La variable d que devuelve la subrutina ludcmp de Numerical recipes vale ±1 según sea la paridad
del número de permutaciones. Así pues, si empleamos esta variable d, el determinante de A será
N
Y
det A = d uii .
i=1

Si la matriz A es tridiagonal, de tamaño N × N :

 
b1 c1 0 0 ...
 a2 b2 c2 0 ... 
 
A =  0 a3 b3 c3 ...  ,
 
 .. .. 
 . . 
bN
F ÍSICA COMPUTACIONAL , Sistemas de ecuaciones lineales. Álgebra matricial 77

se obtiene su determinante mediante la aplicación sistemática de las siguientes relaciones:

N =1 → D1 = b1 ,
N =2 → D2 = b2 D1 − a2 c1 ,
N ≥3 → DN = bN DN −1 − aN cN −1 DN −2 .

Matriz inversa
La inversa A−1 de la matriz A satisface la relación A · A−1 = I. Denominaremos cij a los elementos
de la matriz inversa.
Para calcular A−1 descomponemos las matrices A−1 e I en columnas, de forma que
             
a11 a12 · · · a1N c11 c12 c1N 1 0 0
 a21 a22 · · · a2N   c21   c22   c2N   0   1   0 
..  ·  ..  t  ..  t · · · t  ..  =  ..  t  ..  t · · · t  .. 
             
 .. .. ..
 . . . .   .   .   .   .   .   . 
aN 1 aN 2 ··· aN N cN 1 cN 2 cN N 0 0 1

Podemos obtener la matriz inversa A−1 tras resolver N sistemas de ecuaciones lineales de la forma
     
a11 a12 · · · a1N c1j δ1j
 a21 a22 · · · a2N   c2j   δ2j 
..  ·  ..  =  ..  (para cada columna j) ,
     
 .. .. ..
 . . . .   .   . 
aN 1 aN 2 ··· aN N cN j δN j

donde δij es la función delta de Kronecker:

δij = 0 (i 6= j) ,

δij = 1 (i = j) .
La matriz inversa de A se forma al disponer uno a continuación de otro los vectores columna calculados
al resolver cada uno de los sistemas de ecuaciones anteriores:
     
c11 c12 c1N
 c21   c22   c2N 
A−1 =  .  t  .  t · · · t  .  .
     
 ..   ..   .. 
cN 1 cN 2 cN N

Para resolver el sistema de ecuaciones que corresponde a cada una de las columnas de la matriz A−1 se
emplean las subrutinas de descomposición LU (ludcmp) y de sustitución hacia atrás (lubksb).
78 Rafael Garcia Molina, Departamento de Física – CIOyN, Universidad de Murcia

Ejercicios resueltos
2.1 Escribe el valor de x1 , x2 y x3 que se obtiene en cada iteración k (hasta k . 100) al resolver el
sistema de ecuaciones:
2x1 + x3 = 5 ,
−x1 + 2x2 = 3,
x1 + x2 + x3 = 6 ,
(a) por el método de Jacobi y (b) por el método de Gauss-Seidel. (c) Emplea el método de eliminación
gaussiana para resolver analíticamente el sistema de ecuaciones.
(a) Los valores de x1 , x2 y x3 que se obtienen mediante el método de Jacobi (2-01Jacobi.f) apa-
recen en la tabla adjunta, cuando se parte de los valores iniciales x1 = x2 = x3 = 1.

iteración k x1 x2 x3
1 2.00000 2.00000 4.00000
2 0.50000 2.50000 2.00000
3 1.50000 1.75000 3.00000
4 1.00000 2.25000 2.75000
5 1.12500 2.00000 2.75000
10 1.06250 2.03125 2.89063
20 1.01794 2.01013 2.96826
30 1.00527 2.00298 2.99068
40 1.00155 2.00087 2.99727
50 1.00045 2.00026 2.99920
60 1.00013 2.00008 2.99976
70 1.00004 2.00002 2.99993
80 1.00001 2.00001 2.99998
90 1.00000 2.00000 2.99999

Puede apreciarse que la solución converge lentamente a los valores x1 = 1, x2 = 2 y x3 = 3.

(b) Si se parte de los mismos valores iniciales, mediante el método de Gauss-Seidel (2-01GaussSeidel-rgm.f)
la solución converge más rápidamente, pues en la iteración k ' 45 se llega a la solución x1 = 1, x2 = 2 y
x3 = 3.
Puede comprobarse que si el sistema de ecuaciones se hubiera escrito en la forma
x1 + x2 + x3 = 6,
−x1 + 2x2 = 3,
2x1 + x3 = 5,
no se resolvería bien por el método de Jacobi ni por el de Gauss-Seidel. Este contratiempo puede resolverse
al reordenar las ecuaciones tal como aparecen en el enunciado del ejercicio:
2x1 + x3 = 5,
−x1 + 2x2 = 3,
x1 + x1 + x3 = 6.
Nótese que la nueva reordenaciónP no hace que la matriz de coeficientes sea diagonalmente dominante,
pero ahora la condición |aii | > i6=j |aij | se cumple en dos de las ecuaciones, cuando en el sistema mal
ordenado se cumplía tan sólo en una ecuación.
(c) En primer lugar, hemos de eliminar x1 de las ecuaciones segunda y tercera. Para ello procedemos
como sigue.
A continuación veremos cómo se procede para la eliminación de xi en las ecuaciones j > i.
Eliminación de x1 de la segunda ecuación. Se divide la primera ecuación por a11 = 2 y se multiplica
por a21 = −1, con lo cual se convierte en
−x1 − x3 /2 = −5/2 .
F ÍSICA COMPUTACIONAL , Sistemas de ecuaciones lineales. Álgebra matricial 79

Se resta esta ecuación a la segunda ecuación del sistema y se obtiene:

2x2 + x3 /2 = 11/2 (ésta será la nueva segunda ecuación) .

Eliminación de x1 de la tercera ecuación. Se divide la primera ecuación por a11 = 2 y se multiplica

por a31 = 1, con lo cual se convierte en

x1 + x3 /2 = 5/2 .

Tras restar esta ecuación a la tercera ecuación del sistema, se obtiene:

x2 + x3 /2 = 7/2 (ésta será la nueva tecera ecuación) .

Ahora, el nuevo sistema de ecuaciones es

2x1 + x3 = 5,
1 11
2x2 + 2 x3 = 2 ,
1 7
x2 + 2 x3 = 2.

Eliminación de x2 de la nueva tercera ecuación. Se divide la segunda ecuación por a22 = 2 y se

multiplica por a32 = 1, lo que da lugar

x2 + x3 /4 = 11/4 .

Esta ecuación se resta a la tercera ecuación que se obtuvo anteriormente, lo cual resulta en

x3 /4 = 3/4 (en esto se ha convertido la tercera ecuación) .

Así pues, tras operar sobre el sistema original de ecuaciones, éstas se han transformado en:

2x1 + x3 = 5,
2x2 + x3 /2 = 11/2 ,
x3 /4 = 3/4 .

A partir de la última ecuación se puede ir sustituyendo hacia atrás y se obtiene:

3a ec.: x3 /4 = 3/4 → x3 = 3 ,
2a ec.: 2x2 + x3 /2 = 11/2 → 2x2 + 3/2 = 11/2 → x2 = 2 ,
1a ec.: 2x1 + x3 = 5 → 2x1 + 3 = 5 → x1 = 1 .

Por lo tanto, la solución del sistema de ecuaciones es x1 = 1, x2 = 2 y x3 = 3.

2.10 Calcula el valor de x1 , x2 y x3 al resolver el sistema de ecuaciones:

2.51x1 + 1.48x2 + 4.53x3 = 0.05

1.48x1 + 9.03x2 − 1.30x3 = 1.03
2.68x1 + 0.34x2 − 1.48x3 = −0.53

(a) mediante el método de Jacobi.

(b) mediante el método de Gauss-Seidel.
(c) mediante el método de sobrerrelajación.
(d) mediante el método de descomposición LU .
Utiliza 20 iteraciones en todos los métodos iterativos.

Al aplicar los métodos iterativos al sistema de ecuaciones, los valores que se obtienen para x1 , x2 y x3
divergen.
80 Rafael Garcia Molina, Departamento de Física – CIOyN, Universidad de Murcia

Se reordenan las ecuaciones de la siguiente forma

2.68x1 + 0.34x2 − 1.48x3 = −0.53

1.48x1 + 9.03x2 − 1.30x3 = 1.03
2.51x1 + 1.48x2 + 4.53x3 = 0.05

para que la matriz de coeficientes sea diagonalmente dominante.

A continuación aparecen los resultados que se obtiene mediante cada método (con una tolerancia  =
10−4 y con un número máximo de 20 iteraciones).
(a) Método de Jacobi (2-10Jacobi-rgm.f): x1 = −0.182749, x2 = 0.152982 y x3 = 0.0622746
(tras 20 iteraciones)
(b) Método de Gauss-Seidel (2-10GaussSeidel-rgm.f): x1 = −0.182779, x2 = 0.152987 y
x3 = 0.0623296 (tras 10 iteraciones)
(c) Método de sobrerrelajación (2-10SOR-rgm.f): x1 = −0.182758, x2 = 0.152992 y x3 =
0.0623176 (tras 6 iteraciones, con ω = 0.9)
(d) Método de descomposición LU (2-10LU-rgm.f): x1 = −0.182756513, x2 = 0.152989089 y
x3 = 0.0623167716.

2.11 Emplea algún método iterativo para resolver el sistema de ecuaciones:

6x1 − 3x2 + x3 = 11
2x1 + x2 − 8x3 = −15 .
x1 − 7x2 + x3 = 10

Hay problemas de convergencia porque la matriz resultante no es dominante diagonalmente. Se reorde-

nan las ecuaciones como sigue:

6x1 − 3x2 + x3 = 11
x1 − 7x2 + x3 = 10
2x1 + x2 − 8x3 = −15

A continuación aparecen los resultados que se obtienen mediante cada método (con una precisión  =
10−4 y con un número máximo de 20 iteraciones).
(a) Método de Jacobi (2-11Jacobi-rgm.f): x1 = 1.00003, x2 = −0.999975 y x3 = 2.00003.
(b) Método de Gauss-Seidel (2-11GaussSeidel-rgm.f): x1 = 1.00001, x2 = −0.999998 y
x3 = 2.00000.
(c) Método de sobrerrelajación (2-11SOR-rgm.f): el parámetro para el cual funciona más rápido
es ω = 1, lo cual equivale al método de Gauss-Seidel y, por tanto, se obtienen los mismos resultados
x1 = 1.00001, x2 = −0.999998 y x3 = 2.00000.
(d) Método de descomposición LU (2-11LU-rgm.f): x1 = 0.99999994, x2 = −1.00000012 y
x3 = 2.00.

2.22 Un par termoeléctrico de platino-rodio tiene la siguiente dependencia entre el voltaje V y la

temperatura T :
V = x0 + x1 T + x2 T 2 .
Si disponemos de los siguientes datos:

T (o C) 630.5 960.5 1063.0

V (mV) 5.535 9.117 10.301

¿cuál es la temperatura correspondiente a una lectura de 9.000 mV?

F ÍSICA COMPUTACIONAL , Sistemas de ecuaciones lineales. Álgebra matricial 81

Al sustituir los datos en la ecuación que relaciona V y T , obtenemos un sistema de 3 ecuaciones con 3
incógnitas (x0 , x1 , x2 ):
x0 + 630.5x1 + 630.52 x2 = 5.535
x0 + 960.5x1 + 960.52 x2 = 9.117
x0 + 1063.0x1 + 1063.02 x2 = 10.301 ,
cuya solución obtenemos mediante el método de descomposición LU ; para ello empleamos el programa
2-22LU-rgm.f. El resultado que se obtiene es

x0 = −0.3333 x1 = 8.2918E − 3 x2 = 1.6108E − 6 .

Una vez conocidos los valores de x0 , x1 , x2 , el cálculo de la temperatura T correspondiente a una lectura
de V0 = 9.000 mV se reduce a resolver la ecuación de segundo grado:

x2 T 2 + x1 T + x0 − V0 = 0 ,

cuyo resultado es p
x21 − 4x2 (x0 − V0 )
−x1 ±
T = .
2x2
De aquí se obtiene que T = 950.21 o C y T = −6097.80 o C, pero sólo el primer valor de T tiene sentido
físico.

2.27 Aplica algún algoritmo iterativo para resolver el siguiente sistema de ecuaciones no lineales:
x3 + 3y 2 = 21
2 .
x + 2y + 2 = 0

p
Buscamos gráficamente la solución mediante la intersección de las curvas y = ± (21 − x3 )/3 e y =
−(2 + x2 )/2, que se obtienen de la primera
p y de la segunda ecuación, respectivamente. En la figura adjunta
se muestra que sólo la curva y = − (21 − x3 )/3 se cruza con y = −(2 + x2 )/2 en las proximidades de
los puntos (−2, −3) y (1.6, −2.3).

Mediante un procedimiento iterativo, se obtienen las siguientes soluciones: (−2.0791, −3.1613) y

(1.6431, −2.3499).
Atención, hay que despejar la x como sigue:
x=-sqrt(-2.0d0-2.0d0*y) ! para obtener la primera solución (x<0)
x=+sqrt(-2.0d0-2.0d0*y) ! para obtener la segunda solución (x>0)

2.29 Aplica algún algoritmo iterativo para resolver el siguiente sistema de ecuaciones no lineales:
x2 + x − y 2 = 1
.
y − sin(x2 ) = 0
82 Rafael Garcia Molina, Departamento de Física – CIOyN, Universidad de Murcia

La solución gráfica
√ de este sistema de ecuaciones corresponde a las intersecciones de las curvas y =
sin(x2 ) e y = ± x2 + x − 1, que están en las proximidades de (0.7, 0.5) y (−1.6, 0.4) respectivamente.
Los términos de la primera ecuación se pueden reordenar como sigue

x2 + x − (y 2 + 1) = 0 ,

de donde obtenemos que

p p
−1 ± 1 − 4[−(y 2 + 1)] −1 ± 5 + 4y 2
x= =
2 2
Tras despejar y de la segunda ecuación, llegamos a los dos siguientes sistemas de ecuaciones no lineales:
p
−1 + 5 + 4y 2
x=
2
y = sin(x2 )
y p
−1 − 5 + 4y 2
x=
2
y = sin(x2 )
Para resolver cada sistema de ecuaciones no lineales empleamos el programa 2-29-rgm.f, que aplica
el algoritmo iterativo de Gauss-Seidel. Las soluciones que obtenemos son (0.726, 0.503) y (−1.67, 0.345)
para el primer y el segundo sistema de ecuaciones, respectivamente. Pero es importante el orden en que se
resuelve las ecuaciones, así como la variable que se despeja de cada ecuación.
A continuación aparecen diferentes formas de resolver el sistema de ecuaciones, así como los resultados
que se obtiene en cada caso.

orden de las ecuaciones variable que se despeja en cada ecuación solución (x, y)
 p 
2 2
x +x−y =1 x = −1 + 5 + 4y 2 /2 (0.726, 0.503)
y − sin(x2 ) = 0 y = sin(x2 )
y − sin(x2 ) = 0 y = sin(x

2
)  (0.726, 0.503)
p
2 2
x +x−y =1 x = −1 + 5 + 4y 2 /2
 p 
x2 + x − y 2 = 1 x = −1 − 5 − 4y 2 /2 (−1.67, 0.345)
y − sin(x2 ) = 0 y = sin(x2 )
y − sin(x2 ) = 0 y = sin(x

2
)  (error, 0.345)
p
2 2
x +x−y =1 x = −1 − 5 − 4y 2 /2
F ÍSICA COMPUTACIONAL , Sistemas de ecuaciones lineales. Álgebra matricial 83

Si se despeja como sigue: √

x = √ arcsin y
y = x2 + x − 1
no se obtiene la solución, independientemente del orden en que se resuelva las ecuaciones.

2.41 (i) Resuelve los siguientes sistemas de ecuaciones lineales:

(a)
2x1 − x2 − x3 = 1.1 ,
−x1 + x2 − 3x3 = 0.05 ,
3x1 − 3x2 + 4x3 = 1.9 .

(b)
x1 − 2x2 − x4 = 3,
−2x1 + 7x2 + x3 − x5 = 2,
x2 + 2x3 + x4 + 2x5 = 1,
−x1 − x3 + 6x4 + x5 = 4,
− x2 + 2x3 + x4 + 9x5 = 0.

(c)
5x1 + 4x2 + 2x3 − 3x4 + x5 = 2,
2x1 − 6x2 + x3 + 2x4 − 3x5 = −5 ,
4x1 − x2 + 4x3 − x4 − x5 = −9 ,
x1 + x2 + x3 + x4 + x5 = −1 ,
−x1 + x2 − x3 + x4 − x5 = 7.

Calcula (ii) el determinante de la matriz de coeficientes y (iii) la correspondiente matriz inversa.

(a) El sistema de ecuaciones (A | b) que hay que resolver es:
 
2 −1 −1 1.1
 −1 1 −3 0.05  .
3 −3 4 1.9

Tal como aparecen escritas las ecuaciones no se pueden resolver mediante el método de Gauss-Seidel
(2-41aGaussSeidel-rgm.f). Pero si se reordenan las ecuaciones, permutando la 2a y la 3a ecuacio-
nes, entonces sí converge el método de Gauss-Seidel y se obtiene la siguiente solución: x1 = −0.4899166,
x2 = −1.669916 y x3 = −0.4099997, con una tolerancia de 10−4 .
Si se emplea el método de la descomposición LU (2-41aLU-rgm.f) se obtiene la solución x1 =
−0.49, x2 = −1.67 y x3 = −0.41.
Mediante el programa 2-41aLU-inversa-rgm.f se obtiene que el determinante vale det = 5 y la
matriz inversa vale  
1 −0.8 −1.4
 1 −1.4 −2.2  .
∼ 0 −0.2 −0.6
(b) El sistema de ecuaciones (A | b) que hay que resolver es:
 
1 −2 0 −1 0 3
 −2 7 1 0 −1 2 
A−1 = 
 
 0 1 2 1 2 1  .

 −1 0 −1 6 1 4 
0 −1 2 1 9 0

Mediante el método de Gauss-Seidel (2-41bGaussSeidel-rgm.f), se obtiene la solución x1 = 19.89611,

x2 = 7.137768, x3 = −6.275697, x4 = 2.620654 y x5 = 1.896501, con una tolerancia de 10−4 .
84 Rafael Garcia Molina, Departamento de Física – CIOyN, Universidad de Murcia

Si se emplea el método de la descomposición LU (2-41bLU-rgm.f) se obtiene la solución x1 =

19.89655, x2 = 7.137931, x3 = −6.275862, x4 = 2.620690 y x5 − 1.896552.
Mediante el programa 2-41bLU-inversa-rgm.f se obtiene que el determinante vale det = 116 y
la matriz inversa vale
 
4.7 1.4 −0.47 0.84 0.17
 1.6 0.66 −0.28 0.30 0.10 
A−1 = 
 
 −1.5 −0.58 0.82 −0.35 −0.21  .

 0.46 0.11 0.095 0.23 −0.034 
0.47 0.19 −0.22 0.086 0.17

(c) El sistema de ecuaciones (A | b) que hay que resolver es:

 
5 4 2 −3 1 2

 2 −6 1 2 −3 −5 


 4 −1 4 −1 −1 −9 
 .
 1 1 1 1 1 −1 
−1 1 −1 1 −1 7

Tal como está escrito el sistema, no se puede aplicar el método de Gauss-Seidel (2-41cGaussSeidel-rgm.f),
pues las iteraciones divergen. Tras probar con varias permutaciones entre las ecuaciones, tampoco consigo
resolverlo mediante el método de Gauss-Seidel.
Si se emplea el método de la descomposición LU (2-41cLU-rgm.f) se obtiene la solución x1 = 1,
x2 = 2, x3 = −3, x4 = 1 y x5 = −2.
Mediante el programa 2-41cLU-inversa-rgm.f se obtiene que el determinante vale det = 310 y
la matriz inversa vale
 
0.26 0.23 −0.23 0.14 −0.048
 0.032 −0.097 0.071 0.055 0.31 
−1
 
A =  −0.26
 −0.23 0.43 0.061 0.048  .
 −0.032 0.097 −0.071 0.45 0.19 
0.0 0.0 −0.20 0.30 −0.50

2.45 Resuelve el siguiente sistema de ecuaciones lineales mediante algún método iterativo:

8x1 + x2 − x3 = 8
2x1 + x2 + 9x3 = 12 .
x1 − 7x2 + 2x3 = −4

La solución de este sistema es x1 = x2 = x3 = 1.

Si no reordenamos las ecuaciones, tanto el método de Jacobi (2-45Jacobi-rgm.f) como el de
Gauss-Seidel (2-45GaussSeidel-rgm.f) divergen, como se aprecia en la siguiente tabla, para el caso
del método de Gauss-Seidel:

Iteración x1 x2 x3
0 0.0 0.0 0.0
1 1.000000 10.00000 32.50000
2 3.812500 −288.1250 −1012.344
3 −89.52734 9302.148 32600.28
4 2913.267 −299217.1 −1048719.
5 −93686.67 9625852. 0.3373732E8
10 0.3228039E13 −0.3316617E15 −0.1162430E16
F ÍSICA COMPUTACIONAL , Sistemas de ecuaciones lineales. Álgebra matricial 85

A continuación reordenamos las ecuaciones, intercambiando la segunda y la tercera, para que la matriz
de los coeficientes sea diagonalmente dominante:
8x1 + x2 − x3 = 8
x1 − 7x2 + 2x3 = −4
2x1 + x2 + 9x3 = 12
Ahora sí que convergen ambos métodos. En la tabla adjunta se aprecia claramente que el método de
Gauss-Seidel (segundo conjunto de resultados) converge más rápidamente que el de Jacobi (primer conjunto
de resultados), tal como era de esperar.
Iteración x1 x2 x3 x1 x2 x3
0 0.0 0.0 0.0 0.0 0.0 0.0
1 1.000000 0.5714286 1.333333 1.000000 0.7142857 1.031746
2 1.095238 1.095238 1.047619 1.039683 1.014739 0.9895440
3 0.9940476 1.027211 0.9682540 0.9968506 0.9965627 1.001082
4 0.9926304 0.9900793 0.9982993 1.000565 1.000390 0.9998311
5 1.001027 0.9984613 1.002740 0.9999301 0.9999418 1.000022
6 1.000535 1.000930 0.9999427 1.000010 1.000008 0.9999969
7 0.9998766 1.000060 0.9997779 0.9999986 0.9999989 1.000000
8 0.9999647 0.9999189 1.000021 1.000000 1.000000 0.9999999
9 1.000013 1.000001 1.000017 1.000000 1.000000 1.000000
10 1.000002 1.000007 0.9999971 1.000000 1.000000 1.000000

2.46 Resuelve el siguiente sistema de ecuaciones mediante el método de Gauss-Seidel y por el método
de descomposición LU :     
1 4 1 x1 2
 4 1 0   x2  =  1  .
0 1 4 x3 3

La matriz de coeficientes de las incógnitas ilustra un caso práctico en que el método de Gauss-Seidel
diverge para cualquier valor inicial distinto de la solución exacta, pues la matriz de coeficientes no es
diagonalmente dominante.
Partimos de los valores iniciales x1 = x2 = x3 = 0.
La tabla adjunta muestra los resultados de las iteraciones cuando la matriz de coeficientes no es diago-
nalmente dominante.
1 k x1 x2 x3
2 1 2.0000000 -7.0000000 2.5000000
3 2 27.500000 -109.00000 28.000000
4 3 410.00000 -1639.0000 410.50000
5 4 6147.5000 -24589.000 6148.0000
6 5 92210.000 -368839.00 92210.500
7 6 1383147.5 -5532589.0 1383148.0
8 7 20747210. -82988840. 20747210.
9 8 0.31120816E+09 -0.12448326E+10 0.31120816E+09
10 9 0.46681226E+10 -0.18672490E+11 0.46681226E+10
11 10 0.70021841E+11 -0.28008736E+12 0.70021841E+11
12 11 0.10503276E+13 -0.42013106E+13 0.10503276E+13
13 12 0.15754915E+14 -0.63019661E+14 0.15754915E+14
14 13 0.23632373E+15 -0.94529492E+15 0.23632373E+15
15 14 0.35448561E+16 -0.14179424E+17 0.35448561E+16
16 15 0.53172842E+17 -0.21269137E+18 0.53172842E+17
17 16 0.79759261E+18 -0.31903704E+19 0.79759261E+18
18 17 0.11963889E+20 -0.47855558E+20 0.11963889E+20
19 18 0.17945835E+21 -0.71783339E+21 0.17945835E+21
20 19 0.26918752E+22 -0.10767501E+23 0.26918752E+22
21 20 0.40378130E+23 -0.16151252E+24 0.40378130E+23

Sin embargo, intercambiando las filas primera y segunda, obtenemos una matriz que sí es diagonalmente
dominante. Por lo tanto, converge.
    
4 1 0 x1 1
 1 4 1   x2  =  2 
0 1 4 x3 3
86 Rafael Garcia Molina, Departamento de Física – CIOyN, Universidad de Murcia

En la tabla adjunta aparecen las iteraciones que se obtienen cuando la matriz de coeficientes no es
diagonalmente dominante. Puede comprobarse que el resultado converge hacia los valores 5/28, 2/7 y
19/27, que coinciden con los que proporciona el método de descomposición LU .
1 k x1 x2 x3
2 1 0.25000000 0.43750000 0.64062500
3 2 0.14062500 0.30468750 0.67382813
4 3 0.17382813 0.28808594 0.67797852
5 4 0.17797852 0.28601074 0.67849731
6 5 0.17849731 0.28575134 0.67856216
7 6 0.17856216 0.28571892 0.67857027
8 7 0.17857027 0.28571486 0.67857128
9 8 0.17857128 0.28571436 0.67857140
10 9 0.17857140 0.28571430 0.67857140
11 10 0.17857143 0.28571430 0.67857140
12 11 0.17857143 0.28571430 0.67857140
13 12 0.17857143 0.28571430 0.67857140
14 13 0.17857143 0.28571430 0.67857140
15 14 0.17857143 0.28571430 0.67857140
16 15 0.17857143 0.28571430 0.67857140
17 16 0.17857143 0.28571430 0.67857140
18 17 0.17857143 0.28571430 0.67857140
19 18 0.17857143 0.28571430 0.67857140
20 19 0.17857143 0.28571430 0.67857140
21 20 0.17857143 0.28571430 0.67857140

2.48 (a) Escribe el valor de x1 , x2 y x3 que se obtiene en cada iteración k (desde k = 1 hasta 30) al
resolver el sistema de ecuaciones:
x1 + 5x2 − x3 = −13 ,
6x1 − 2x2 + 2x3 = 8,
7x2 − x3 = 11 ,

por el método de Gauss-Seidel. Compara este resultado con el que obtienes cuanso se emplea el método de
eliminación gaussiana para resolver el sistema de ecuaciones.
La matriz de coeficientes de las incógnitas ilustra un caso práctico en que el método de Gauss-Seidel
diverge para cualquier valor inicial distinto de la solución exacta, pues la matriz de coeficientes no es
diagonalmente dominante.
Partimos de los valores iniciales x1 = x2 = x3 = 0.
La tabla adjunta muestra los resultados de las iteraciones cuando la matriz de coeficientes no es diago-
nalmente dominante.
1 k x1 x2 x3
2 1 -13.000000 -43.000000 -312.00000
3 2 -110.00000 -646.00000 -4533.0000
4 3 -1316.0000 -8485.0000 -59406.000
5 4 -16994.000 -110392.00 -772755.00
6 5 -220808.00 -1435183.0 -10046292.
7 6 -2870390.0 -18657466. -0.13060227E+09
8 7 -37314956. -0.24254714E+09 -0.16978299E+10
9 8 -0.48509414E+09 -0.31531123E+10 -0.22071785E+11
10 9 -0.63062241E+10 -0.40990458E+11 -0.28693319E+12
11 10 -0.81980899E+11 -0.53287590E+12 -0.37301314E+13
12 11 -0.10657519E+13 -0.69273869E+13 -0.48491707E+14
13 12 -0.13854772E+14 -0.90056019E+14 -0.63039215E+15
14 13 -0.18011207E+15 -0.11707284E+16 -0.81950982E+16
15 14 -0.23414562E+16 -0.15219467E+17 -0.10653627E+18
16 15 -0.30438931E+17 -0.19785306E+18 -0.13849715E+19
17 16 -0.39570619E+18 -0.25720900E+19 -0.18004630E+20
18 17 -0.51441806E+19 -0.33437173E+20 -0.23406022E+21
19 18 -0.66874355E+20 -0.43468328E+21 -0.30427831E+22
20 19 -0.86936671E+21 -0.56508832E+22 -0.39556183E+23
21 20 -0.11301768E+23 -0.73461483E+23 -0.51423037E+24

Si se reordenan las filas de la siguiente forma: 1 → 3, 2 → 1 y 3 → 2, la matriz de coeficientes sí es

diagonalmente dominante. El resultado que se obtiene aparece en la tabla adjunta.
1 k x1 x2 x3
2 1 1.3333334 1.5714285 22.190475
3 2 -5.5396824 4.7414966 31.167801
F ÍSICA COMPUTACIONAL , Sistemas de ecuaciones lineales. Álgebra matricial 87

4 3 -7.4754348 6.0239716 35.644424

5 4 -8.5401506 6.6634893 37.777298
6 5 -9.0379362 6.9681854 38.802990
7 6 -9.2782679 7.1147127 39.295296
8 7 -9.3935280 7.1850424 39.531685
9 8 -9.4488802 7.2188120 39.645180
10 9 -9.4754562 7.2350259 39.699673
11 10 -9.4882154 7.2428102 39.725838
12 11 -9.4943428 7.2465482 39.738396
13 12 -9.4972830 7.2483420 39.744427
14 13 -9.4986954 7.2492037 39.747322
15 14 -9.4993734 7.2496176 39.748714
16 15 -9.4996986 7.2498164 39.749382
17 16 -9.4998550 7.2499118 39.749702
18 17 -9.4999304 7.2499576 39.749855
19 18 -9.4999657 7.2499795 39.749931
20 19 -9.4999838 7.2499900 39.749966
21 20 -9.4999924 7.2499952 39.749985

El resultado que se obtiene mediante el método LU es: x1 = −9.50, x2 = 7.25 y x3 = 39.75.

(b) Repite el ejercicio, pero ahora considerando el siguiente sistema de ecuaciones (solo cambia el
coeficiente de x3 en la primera ecuación):

x1 + 5x2 − 8x3 = −13 ,

6x1 − 2x2 + 2x3 = 8,
7x2 − x3 = 11 ,

En este caso, el resultado es x1 = 1, x2 = 2 y x3 = 3. Pero si se emplea el método de Gauss-Seidel para

resolver el sistema de ecuaciones tal como está escrito, aparecen inmediatamente divergencias. En cambio,
si se reordenan las ecuaciones como sigue: 1 → 3, 2 → 1 y 3 → 2, entonces la matriz de coeficientes es
diagonalmente dominante y mediante el método de Gauss-Seidel se obtiene fácilmente que la solución es
x1 = 1, x2 = 2 y x3 = 3.
Al mismo resultado se llega mediante el método de descomposición LU .