CAPITULO 3 – SISTEMAS DE ECUACIONES

3.1 SISTEMAS DE ECUACIONES LINEALES

Un gran número de problemas prácticos de ingeniería se reduce al problema de

resolver un sistema de ecuaciones lineales. Por ejemplo, puede citarse la solución
de sistemas de ecuaciones no lineales, la aproximación polinomial, la solución de
ecuaciones diferenciales parciales, entre otros.
Un sistema de m ecuaciones lineales con n incógnitas tiene la forma general
a11 p1  a12 p 2  ...  a1n p n  q1
a 21 p1  a 22 p 2  ...  a 2 n p n  q 2 3.1
... ... ... .
a m1 p1  a m 2 p 2  ...  a mn p n  q m
Con la notación matricial se puede escribir la ecuación anterior como
 a11 a12 ... a1n   p1   q1 
a     
 21 a 22 ... a 2 n  x  p 2  =  q 2 
 ... ... ... ...   ...   ... 
     
a m1 a m 2 ... a mn   pn  q m 
Y concretamente como A p = q., donde A es la matriz coeficiente del sistema, p
el vector incógnita y q el vector de términos independientes.
Dados A y q, se entiende por resolver el sistema (Ec. 3.1) encontrar los vectores
p que lo satisfagan. A continuación estudiaremos las técnicas que permiten
encontrar p mediante los métodos directos y los métodos iterativos.

3.2 MÉTODOS DIRECTOS DE SOLUCIÓN

Son aplicables a sistemas de ecuaciones lineales de tamaño pequeño o

mediano. La particularidad es que obtienen las soluciones exactas mediante
operaciones algébricas y trabajando con todo el sistema a la vez. Los sistemas
grandes presentan la inconveniencia de requerir mucha memoria del computador, la
cual a veces puede ser insuficiente. Han sido planteados diversos métodos de
solución y la literatura sobre solución de sistemas lineales está enriquecida con
interesantes aportes.
A continuación se describirán algunos de los métodos mas aplicables a la
simulación matemática de reservorios con una relativamente pequeña cantidad de
bloques.
El prototipo de todos estos métodos se conoce como la eliminación de Gauss
y se presenta a continuación.

3.2.1 ELIMINACIÓN DE GAUSS

Considérese un sistema general de tres ecuaciones lineales con tres

incógnitas
a11p1 + a12p2 + a13p3 = q1
a21p1 + a22p2 + a23p3 = q2 3.2
a31p1 + a32p2 + a33p3 = q3
Que es representado de la siguiente manera:

Ing. Hermas Herrera Callejas Página : 1 de 33

 Programación Aplicada Capítulo 3–Sistemas de Ecuaciones

 a11  p1 
a12 a13 
 q1 
a a 22   a 23 
 
 21  p 2  = q 2 
a31a32  p3  a33 
 q3 
A p = q
Básicamente este método tiene el objetivo de convertir la matriz A en una
matriz triangular, cuyos elementos inferiores son ceros; para ello hay que multiplicar
convenientemente una fila por un elemento y sumarla a otra.
Como primer paso, se remplaza la segunda ecuación con lo que resulte de
sumarle la primera ecuación multiplicada por (-a21/a11), con lo que se obtiene un 0 en
la posición de a21. Similarmente se sustituye la tercera ecuación con el resultado de
sumarle la primera ecuación multiplicada por (-a31/a11).
Esto da lugar al nuevo sistema:
a11 a12 a13   p1   q1 
 0 a a a /a a 23  a13 a 21 / a11   p 2  = q 2  q1 a 21 / a11 
   
 22 12 21 11

 0 a32  a12 a31 / a11 a33  a13 a31 / a11   p3   q 3  q1 a 31 / a11 

Simbólicamente los nuevos valores pueden representarse más simplificadamente así:
a11a12 a13   p1   q1 
0 '
a 22 '  
a 23 p   '
   2  = q 2  3.3
 0a32 a33   p3   q3 
' ' '

Donde las a’ y las q’ son los nuevos elementos que se obtienen de las
operaciones ya mencionadas, y donde p 1 se ha eliminado en la segunda y tercera
ecuaciones. Ahora, multiplicando la segunda ecuación de 3.3 por (-a 32‘/a22’) y
sumando el resultado a la tercera ecuación de 3.3, se obtiene el sistema triangular:
a11 a12 a13   p1   q1 
 0 a' '
a 23     q2 ' 
  2=
p
 22 
 0 0 a33'  a 23
'
a32' / a 22
'
  p3  q3'  q 2' a32' / a 22' 
Que simbólicamente puede representarse más simplificadamente así:
 p1 
p 
 2 = 3.4
a11 a12 a13   q1 
 ' '   ' 
 0 a 22 a 23  q 2 

 0 0 "
a 33 
  q3

"



 p3 
Donde a33” y q3“, resultaron de las operaciones realizadas y P2 se ha eliminado
de la tercera ecuación.
El proceso de llevar el sistema de ecuaciones 3.2 a la forma de la ecuación
3.4 se conoce como triangularización.
El sistema en la forma de la ecuación 3.4 se resuelve despejando de su última
ecuación p3, sustituyendo p3 en la segunda ecuación y despejando P2 de ella, por
último, con p3 y p2 sustituidas en la primera ecuación de 3.4 se obtiene p 1. Esta parte
del proceso se llama sustitución regresiva.
Antes de ilustrar la eliminación de Gauss con un ejemplo particular, nótese
que no es necesario conservar p 1, p2 y p3 en la triangularización y que ésta puede
llevarse a cabo usando solamente la matriz coeficiente A y el vector q. Para mayor
simplicidad se empleará la matriz aumentada B.

Ing. Hermas Herrera Callejas Página: 2 de 32

Programación Aplicada Capítulo 3–Sistemas de Ecuaciones

 a11 a12 a13 q1 

 
B =  a 21 a 22 a 23 q 2  = [A | q]
 a31 a32 a33 q3 
Con esto se incorporan la notación matricial y todas sus ventajas a la solución
de sistemas de ecuaciones lineales.

Ejemplo 3.1
Resuelva por eliminación de Gauss el sistema:
4p1 - 9p2 + 2p3 = 500
2p1 - 4p2 + 6p3 = 300
p1 - p2 + 3p3 = 400

SOLUCIÓN
La matriz aumentada del sistema es
4 9 2 500 
2 4 6 300 

1 1 3 400
Triangularización.- Al sumar la primera ecuación multiplicada por (-2/4) a la
segunda, y la primera ecuación multiplicada por (-1/4) a la tercera, resulta:
4 9 2 500
0 0.5 5 50 
 

0 1.25 2.5 275

Obsérvese que en este paso la primera fila se conserva sin cambio.
Sumando la segunda fila multiplicada por (-1.25/0.5) a la tercera se obtiene la matriz
4 9 2 500
0 0 .5 5 50 
 

0 0  10 150 

Que en términos de sistemas de ecuaciones quedaría como
4p1 - 9p2 + 2p3 = 500
0.5p2 + 5p3 = 50 3.5
-10p3 = 150
Un proceso de sustitución regresiva produce el resultado buscado. La tercera
ecuación de 3.5 da el valor de p3 = -15. De la segunda ecuación se obtiene entonces
0.5P2 = 50 - 5p3 = 125. Y por tanto P2 = 250
Finalmente al sustituir p2 y p3 en la primera ecuación de la forma 3.5 resulta
4p1 = 500 + 9p2 – 2p3 = 2780. De modo que p1 = 695
Con la sustitución de estos valores en el sistema original se verifica la
exactitud de los resultados.

3.2.2 ELIMINACIÓN DE GAUSS CON PIVOTEO

En la eliminación de p1 de la segunda y tercera ecuaciones de la forma 3.2 se

tomó como base la primera, por lo cual se denomina ecuación pivote o, en términos
de la notación matricial, fila pivote. Para eliminar P2 de la tercera ecuación de la
forma 3.3, la fila pivote utilizada fue la segunda. El coeficiente de la incógnita que se
va a eliminar en la fila pivote se llama pivote. En la eliminación que dio como
resultado el sistema de ecuaciones 3.4, los pivotes fueron a 11 y a22’. Esta elección

Ing. Hermas Herrera Callejas Página: 3 de 32

Programación Aplicada Capítulo 3–Sistemas de Ecuaciones

natural de los pivotes a 11, a22’, a33”, etc., es muy conveniente tanto para trabajar con
una calculadora como con una computadora, desafortunadamente falla cuando
alguno de esos elementos es cero, puesto que los multiplicadores quedarían
indeterminados (por ejemplo si a 11 fuera cero, el multiplicador (a21/a11) no está
definido). Una manera de evitar esta posibilidad es seleccionar como pivote el
coeficiente de máximo valor absoluto en la columna relevante de la matriz reducida.
Como antes, se tomarán las columnas en orden natural de modo que se vayan
eliminando las incógnitas también en orden natural p1, p2, p3, etc. Esta técnica,
llamada pivoteo, se ilustra con la solución del siguiente sistema.

Ejemplo 3.2
Resuelva el sistema
10p1 + p2 - 5p3 = 100
-20p1 + 3p2 + 20p3 = 200 3.6
5p1 + 3p2 + 5p3 = 600
SOLUCIÓN
La matriz aumentada es
 10 1 5 100 
 20 3 20 200

 5 3 5 600 
El primer pivote debe ser (-20), ya que es el elemento de máximo valor absoluto en la
primera columna. Vamos a la forma triangular en la eliminación. Para esto es
necesario, por ejemplo en la ecuación 3.6, intercambiar la segunda fila (donde se
encuentra el elemento de máximo valor absoluto) con la primera, con lo que se
obtiene
 20 3 20 200
 10 1 5 100 

 5 3 5 600 
Se elimina entonces P1 de la segunda y tercera filas de la ecuación 3.6. Para
ello, se suma a la segunda fila la primera multiplicada por (-10/(-20)), y a la tercera
fila la primera multiplicada por (-5/(-20)). Con esto se reduce en la primera
eliminación a
 20 3 20 200
 0 2.5 5 200

 0 3.75 10 650
El siguiente pivote debe seleccionarse entre la segunda y tercera filas (segunda
columna) y en este caso es (3.75). se intercambian la segunda y la tercera filas de la
matriz anterior para obtener
 20 3 20 200
 0 3.75 10 650

 0 2.5 5 200
Sumando a la tercera fila la segunda multiplicada por (-2.5/3.75), al eliminar p 2
resulta
 20 3 20 200 
 0 3.75 10 650 
 

 0 0  1.666  233.3


Ing. Hermas Herrera Callejas Página: 4 de 32

Programación Aplicada Capítulo 3–Sistemas de Ecuaciones

Que tiene ya la forma triangular y está lista para la sustitución regresiva que
proporciona los siguientes valores:
 233.3
p3   140 . De la tercera ecuación
 1.666
650  10(140)
p2   200 . De la segunda ecuación. Y finalmente:
3.75
200  3( 200)  20(140)
p1   100 . De la primera ecuación
 20
Para terminar el tema, comparemos las técnicas de eliminación de Gauss con
pivoteo y sin éste. Por brevedad, la primera se denominará GP y la segunda G.
1. La búsqueda del coeficiente de mayor valor absoluto que se usará como pivote y
el intercambio de filas significa mayor programación en GP.
2. Encontrar en GP un pivote igual a cero significaría que se trata de una matriz
coeficiente A singular y que el sistema A p = q no tiene solución única. Encontrar en
G un pivote igual a cero detendría el proceso de triangularización.
A pesar de la programación adicional y el mayor tiempo de máquina que se
emplea en el método de Gauss con pivoteo, sus otras ventajas borran totalmente
estas desventajas en la práctica.

3.2.3 ELIMINACIÓN DE GAUSS - JORDAN

Es posible extender los métodos vistos de modo que las ecuaciones se

reduzcan a una forma en que la matriz coeficiente del sistema sea diagonal y ya no
se requiera la sustitución regresiva. Los pivotes se eligen como en el método de
Gauss con pivoteo, y una vez intercambiadas las filas se eliminan los elementos
arriba y debajo del pivote. El sistema del ejemplo 3.3 ilustra este método.

Ejemplo 3.3
Por eliminación de Jordan, resuelva el sistema
4p1 - 9p2 + 2p3 = 500
2p1 - 4p2 + 6p3 = 300
p1 - p2 + 3p3 = 400

SOLUCION
La matriz aumentada del sistema es
4 9 2 500 
2 4 6 300 

1 1 3 400
Como en la primera columna el elemento de máximo valor absoluto se
encuentra en la primera fila, ningún intercambio es necesario y el primer paso de
eliminación produce
4 9 2 500
0 0.5 5 50 
 

0 1.25 2.5 275

El elemento de máximo valor absoluto en la parte relevante de la segunda columna
(filas 2 y 3) es 1.25; por tanto, la fila 3 debe intercambiarse con la 2.

Ing. Hermas Herrera Callejas Página: 5 de 32

Programación Aplicada Capítulo 3–Sistemas de Ecuaciones

4 9 2 500
0 1.25 2.5 275
 

0 0.5 5 50 

Sumando la segunda fila multiplicada por (-(-9)/1.25), a la primera fila y la segunda
multiplicada por (-0.5/1.25) a la tercera, se obtiene el nuevo arreglo
4 0 20 2480
0 1.25 2.5 275 
 

0 0 4  60 

Donde se han eliminado los elementos de arriba y abajo del pivote (nótese que en
este paso el primer pivote no se modifica porque sólo hay ceros debajo de él). Por
último, sumando la tercera fila multiplicada por (-20/4) a la primera fila y la tercera
multiplicada por (-2.5/4) a la segunda
4 0 0 2780 
0 1.25 0 312.5
 

0 0 4  60 

Que escrita de nuevo como sistema de ecuaciones da
4p1 = 2780
1.25p2 = 312.5
4p3 = -60
De donde el resultado final se obtiene fácilmente
p1  2780 / 4  695
p 2  312.5 / 1.25  250
p 3  60 / 4  15

3.2.4 MÉTODOS DE FACTORIZACIÓN DE DOOLITLE Y CROUT.

Consiste en descomponer la matriz A en dos matrices: una matriz triangular

inferior L y una matriz triangular superior U para aplicarse al sistema A p = q sin
intercambio de filas. El resultado anterior permite resolver el sistema A p = q, ya que
sustituyendo A por LU se tiene
LUp = q
Se hace U p = g, donde g es un vector desconocido  g1 g 2 g 3 ... g n  , que
T

se puede obtener fácilmente resolviendo el sistema

Lg = q
Con sustitución progresiva o hacia adelante, ya que L es triangular inferior.
Una vez calculado g, se resuelve
Up = g
Con sustitución regresiva, ya que U es triangular superior y de esa manera se
obtiene el vector solución p.
Para encontrar las matrices triangulares se analiza la factorización de A en las
matrices generales L y U, dadas a continuación
 l1,1 0 0 u1,1 u1, 2 u1,3   a1,1 a1, 2 a1,3 
     
l 2,1 l 2, 2 0  0 u 2, 2 u 2,3  = a 2 ,1 a 2, 2 a 2 ,3  3.7
l 3,1 l 3, 2 l 3,3   0 0 u 3,3   a3,1 a 3, 2 a3,3 
Se multiplican
a)Primera fila de L por las tres columnas de U

Ing. Hermas Herrera Callejas Página: 6 de 32

Programación Aplicada Capítulo 3–Sistemas de Ecuaciones

l1,1u1,1  a1,1
l1,1u1, 2  a1, 2
l1,1u1,3  a1, 3
b)Segunda fila de L por las tres columnas de U
l 2 ,1u1,1  a 2 ,1
l 2,1u1, 2  l 2, 2 u 2, 2  a 2, 2
l 2,1u1, 3  l 2 , 2 u 2, 3  a 2 ,3
c) Tercera fila de L por las tres columnas de U
l 3,1u1,1  a 3,1
l 3,1u1, 2  l 3, 2 u 2 , 2  a3, 2
l 3,1u1, 3  l 3, 2 u 2, 3  l 3,3 u 3, 3  a 3, 3
Se llega a un sistema de
ecuaciones con 12 incógnitas nueve
l1,1 , l 2,1 , l 2, 2 , l 3,1 , l 3, 2 , l 3, 3 , u1,1 , u1, 2 , u1, 3 , u 2, 2 , u 2 , 3 , u 3, 3 ,
por lo que será necesario
establecer tres condiciones arbitrarias sobre las incógnitas para resolver dicho
sistema. La forma de seleccionar las condiciones ha dado lugar a diferentes
métodos; por ejemplo, si se toman de modo que l1,1  l 2, 2  l 3,3  1 , se obtiene el
método de Doolitle; si en cambio se selecciona u1,1  u 2, 2  u 3,3  1 , el algoritmo
resultante es llamado método de Crout.
Se continuará el desarrollo de la factorización. Tómese l1,1  l 2, 2  l3,3  1
Con estos valores se resuelven las ecuaciones directamente en el orden en que
están dadas
De (a) u1,1 = a1,1 u1,2 = a1,2 u1,3 = a1,3 3.8
De (b) y sustituyendo los resultados (Ec. 3.8)
a 2,1 a 2,1
l 2,1  
u1,1 a1,1
a 2,1
u 2, 2  a 2, 2  l 2,1u1, 2  a 2, 2  a1, 2 3.9
a1,1
a 2,1
u 2,3  a 2 ,3  l 2,1u1,3  a 2,3  a1,3
a1,1
De (c) y sustituyendo los resultados de las ecuaciones 3.8 y 3.9
a3,1 a3,1
l3,1  
u1,1 a1,1
a 3,1
a 3, 2  a1, 2
a3, 2  l 3,1u1, 2 a1,1
l 3, 2   3.10
u 2, 2 a 2,1
a 2, 2  a1, 2
a1,1
 a3,1 
 a 3, 2  a1, 2 
u 3,3  a3,3  l 3,1u1,3  l 3, 2 u 2,3  a3,3 
a3,1
a1,3  
a1,1  a  a 2,1 a 
a1,1  a 2,1   2,3 a1,1 1,3 
a
 2, 2  a1, 2 
 a1,1 
Las ecuaciones 3.8, 3.9 y 3.10, convenientemente generalizadas constituyen un
método directo para la obtención de L y U, con la ventaja sobre la triangularización

Ing. Hermas Herrera Callejas Página: 7 de 32

Programación Aplicada Capítulo 3–Sistemas de Ecuaciones

de que no se tiene que escribir repetidamente las ecuaciones o arreglos modificados

de Ap = q. A continuación se resuelve un ejemplo.

Ejemplo 3.4
Resuelva por el método de Doolitle el sistema
4p1 -9p2 + 2p3 = 5
2p1 -4p2 + 6p3 = 3
p1 -p2 + 3p3 = 4

SOLUCIÓN
Con l1,1  l 2, 2  l3,3  1 , se procede al cálculo de la primera fila de U
u1,1 = 4 u1,2 = -9 u1,3 = 2
Cálculo de la primera columna de L
l1,1 = 1 (dato) l2,1 = 2/4 = 0.5 l3,1 = ¼ = 0.25
Cálculo de la segunda fila de U
u2,1 = 0 (recuérdese que U es triangular superior)
u2,2 = -4 –(2/4)(-9) = 0.5
u2,3 = 6 –(2/4)(2) = 5
Cálculo de la segunda columna de L
l1,2 = 0 (ya que L es triangular inferior)
l2,2 = 1 (dato)
l3,2 = (-1-(1/4)(-9))/(-4-(2/4)(-9)) = 2.5
Cálculo de la tercera fila de U, o más bien sus elementos faltantes, ya que por
ser triangular superior
u31 = u32 = 0
u3,3 = 3 -(1/4)(2) -[(-l-(1/4)(-9))/(-4-(2/4)(-9))](6-(2/4)(2)) = - 10
Con esto se finaliza la factorización.
Las matrices L y U quedan como sigue
 1 0 0 4 9 2 
 0 .5 0 0 5 
L=  1  U=  0. 5 
0.25
 2 .5 1
 0
 0 
 10
Cuyo producto, como ya se comprobó, da A.
Se resuelve el sistema Lg = q, donde q es el vector de términos independientes
del sistema original
 1 0 0  g 1  5 
 0 .5 1

0 g 2 = 
  
    3

0.25 2.5 1   g 3  
 4

g1 = 5
g2 = 3 –0.5(5) = 0.5
g3 = 4 –0.25(5) –2.5 (0.5) = 1.5
Y, finalmente, al resolver el sistema Up = g se tiene la solución del sistema
original
4  9 2   p1   5 
0 0.5    0.5
 p2  = 
 5  

0 0  10
  p 3  
1.5 


Ing. Hermas Herrera Callejas Página: 8 de 32

Programación Aplicada Capítulo 3–Sistemas de Ecuaciones

p3 = -0.15
p2 = (0.5 –5(-0.15))/0.5 = 2.5
p1 =(5 + 9(2.5) -2(-0.15))/4 = 6.95
 6.95 
 2.5 
p =  

  0.15

Las ecuaciones 3.8, 3.9 y 3.10 se generalizan para factorizar la matriz coefi-
ciente del sistema Ap = q, que puede resolverse por eliminación de Gauss sin
intercambio de filas; se tiene entonces
i 1
u i , j  ai , j   li ,k u k , j j = i+1, … ,n
k 1
j 1
1
li , j  ( ai , j   u k , j l i ,k ) i = j+1, … ,n 3.11
ui, j k 1

l i ,i  1 i = i+1, … ,n
0
Con la convención en las sumatorias que  0
k 1
Puede observarse al seguir las ecuaciones 3.8, 3.9, 3.10 o bien las ecuaciones
3.11, que una vez empleada ai,j el cálculo de ui,j o li,j según sea el caso, esta
componente de A no vuelve a emplearse como tal, por lo que las componentes de L
y U generadas pueden guardarse en A y ahorrar memoria de esa manera.

3.3 MÉTODOS ITERATIVOS

Al resolver un sistema de ecuaciones lineales por eliminación, la memoria de

máquina requerida es proporcional al cuadrado del orden de A, y el trabajo com-
putacional es proporcional al cubo del orden de la matriz coeficiente A. Debido a
esto, la solución de sistemas lineales grandes (n  50), con matrices coeficiente
densas, se vuelve costoso y difícil en una computadora con los métodos de
eliminación, ya que se requiere amplia memoria. Además, como el número de
operaciones que se debe ejecutar es muy grande, se pueden producir errores de
redondeo también muy grandes. Sin embargo, se han resuelto sistemas de orden
1000, y aun mayor, con los métodos que se estudiarán en esta sección.
Estos sistemas de un número muy grande de ecuaciones se presentan en la so-
lución numérica de ecuaciones diferenciales parciales, en la solución de los modelos
resultantes en la simulación de columnas de destilación, etc. En favor de estos sis -
temas, puede decirse que tienen matrices con pocos elementos distintos de cero y
que éstas poseen ciertas propiedades (simétricas, bandeadas, diagonal dominantes,
etc.), que permiten garantizar éxito en la aplicación de los métodos de esta sección.
La solución se obtiene mediante aproximaciones sucesivas; la ventaja que tienen
sobre los métodos directos es que pueden manejar sistemas grandes de ecuaciones,
porque no es necesario que todas las ecuaciones estén en la memoria del
computador. La desventaja relativa radica en el hecho de que los métodos pueden
crear inestabilidad y perder la convergencia de las soluciones.
Considérese el siguiente sistema:
Ap = q
Asumiendo que se obtiene una solución aproximada p 0, esta aproximación
tendrá un error ε0 respecto a la solución exacta, tal que:

Ing. Hermas Herrera Callejas Página: 9 de 32

Programación Aplicada Capítulo 3–Sistemas de Ecuaciones

p = p0 + ε0 Por lo tanto:
Ap = Ap0 + Aε0 O sea:
Aε0 = q - Ap0
Este vector se llama vector residual, y naturalmente es 0 para solución exacta.

3.3.1 MÉTODOS DE JACOBI Y GAUSS-SEIDEL

Se puede aplicar esta técnica para elaborar métodos de solución de Ap = q, de

la manera siguiente.
Se parte de Ap = q para obtener la ecuación
Ap - q = 0, 3.12
Ecuación vectorial correspondiente a f(p) = 0. Se busca ahora una matriz B y
un vector c, de manera que la ecuación vectorial
p = Bp + c, 3.13
Sea sólo un arreglo de la ecuación 3.12; es decir, de manera que la solución de
una sea también la solución de la otra. La ecuación 3.13 correspondería a p = g(p).
A continuación se propone un vector inicial p (0) como primera aproximación al vector
solución p. Luego, se calcula con la ecuación 3.14 la sucesión vectorial p (1), p(2), ...,de
la siguiente manera
p(k+1) = Bp(k) +c, k = 0, 1, 2, ... Donde
(k )
 k k

k T (k)
p  p1 p 2 ... p n p 3.14
(0) (1) (n)
Para que la sucesión p , p , ..., p , ..., converja al vector solución p es
necesario que eventualmente p jm, 1  j  n (los componentes del vector p (m)), se
aproximen tanto a pj, 1  j  n (los componentes correspondientes a p) que todas las
diferencias p j  p j , 1  j  n sean menores que un valor pequeño previamente
m

fijado, y que se conserven menores para todos los vectores siguientes de la iteración;
es decir
lim p mj  p j 1 j n 3.15
m 
La forma como se llega a la ecuación 4.13 define el algoritmo y su
convergencia. Dado el sistema Ap = q, la manera más sencilla es despejar p 1 de la
primera ecuación, p2 de la segunda, etc. Para ello, es necesario que todos los
elementos de la diagonal principal de A, por razones obvias, sean distintos de cero.
Para ver esto en detalle considérese el sistema general de tres ecuaciones
(naturalmente puede extenderse a cualquier número de ecuaciones).
Sea entonces
a11p1 + a12p2 + a13p3 = q1
a21p1 + a22p2 + a23p3 = q2
a31p1 + a32p2 + a33p3 = q3
Con , a11, a22 y a33 distintos de cero.
Se despeja p1 de la primera ecuación, p2 de la segunda y p3 de la tercera con lo que
se obtiene

Ing. Hermas Herrera Callejas Página: 10 de 32

Programación Aplicada Capítulo 3–Sistemas de Ecuaciones

a12 a13 q1
p1   p2  p3 
a11 a11 a11
a 21 a q
p2   p1  23 p3  2 3.16
a 22 a 22 a 22
a a q
p3   31 p1  32 p 2  3
a33 a33 a33
Que en notación matricial queda
 a12 a13   q1 
 0     
 a11 a11   a11 
 p1   p1 
 a 21 a  q2 
p 
= 0  23  p 
+ 4.17
  2
 a a 22    2
a 
 p3   a22   p3   q22 
  31 a32
 0   3
 a33 a33   a33 
Y ésta es la ecuación 3.14 desarrollada, con
 a12 a13   q1 
 0     
 a11 a11   a11 
 a 21 a  q2 
B=  0  23  y c =
 a a 22  a 
 a22 a32   q22 
  31  0   3
 a33 a33   a33 
Una vez que se tiene la forma 3.17, se propone un vector inicial p (0) que puede
ser p(0) = 0, o algún otro que sea aproximado al vector solución p. Se presenta a
continuación un algoritmo para resolver sistemas de ecuaciones lineales con el
método iterativo, en sus dos versiones, desplazamientos simultáneos y
desplazamientos sucesivos

[Link] ITERACIÓN DE JACOBI (DESPLAZAMIENTOS SIMULTÁNEOS)

 p1k 
 
Si p (k )   p 2k  3.18
 p3k 
 
Es el vector aproximación a la solución p después de k iteraciones, entonces se
tiene para la siguiente aproximación
 1 
 (q1  a12 p 2k  a13 p3k ) 
 p1k 1   a11 
( k 1)  k 1   1 k 
p   p2   (q  a p  a 23 p3 )
k
3.19
 a 22 2 21 1 
 p3  
k 1

  1
 (q3  a31 p1k  a32 p 2k ) 
 a33 

Ing. Hermas Herrera Callejas Página: 11 de 32

Programación Aplicada Capítulo 3–Sistemas de Ecuaciones

O bien, para un sistema de n ecuaciones con n incógnitas y usando notación

más compacta y de mayor utilidad en programación, se tiene
 
1  n
k
p k 1
  qi   aij p j  , para 1  i  n 3.20
aii 
i
j 1
 
 j i 

[Link] ITERACIÓN DE GAUSS-SEIDEL (DESPLAZAMIENTOS SUCESIVOS)

En este método los valores que se van calculando en la (k+1)-ésima iteración se

emplean para calcular los valores faltantes de esa misma iteraciónes decir, con
p ( k ) se calcula p ( k 1) de acuerdo con
 1 
 (q1  a12 p 2k  a13 p3k ) 
 p1k 1   a11 
( k 1)  k 1   1 k 1 k 
p   p2   (q  a 21 p1  a 23 p3 ) 3.21
 a 22 2 
 p3  
k 1

  1
 (q3  a31 p1k 1  a32 p 2k 1 )
 a33 
O bien, para un sistema de n ecuaciones
1  i 1 n 
pik 1    qi   aij p kj 1   aij p kj  , para 1  i  n 3.22
a ii  j 1 j i 1 

Ejemplo 3.5
Resuelva el siguiente sistema por los métodos de Jacobi y Gauss-Seidel
4 p1  p2 5
 p1  4 p2  p3 9
3.23
 p2  4 p3  p4  3
 p3  4 p4  15
SOLUCIÓN
Despejando p1 de la primera ecuación, p2 de la segunda, etc., se obtiene
p2 5
p1  
4 4
p1 p3 9
p2   
4 4 4 3.24
p2 p 3
p3   4 
4 4 4
p3 15
p4   
4 4
Vector inicial.- Cuando no se tiene una aproximación al vector solución, se
emplea generalmente como vector inicial el vector cero, esto es p ( 0 )   0 0 0 0 T

a) Método de Jacobi
El cálculo de p (1) en el método de Jacobi se obtiene remplazando p ( 0) en cada
una de las ecuaciones de 3.24

Ing. Hermas Herrera Callejas Página: 12 de 32

Programación Aplicada Capítulo 3–Sistemas de Ecuaciones

0 5 5
p1     1.25
4 4 4
0 0 9 9
p2      2.25
4 4 4 4
0 0 3 3
p3       0.75
4 4 4 4
0 15 15
p4      3.75
4 4 4
T
  5 9 3 15
Y entonces p (1)   , , , 
 4 4 4 4 
( 2)
Para calcular p se sustituye (1) en cada una de las ecuaciones de 3.24. Para
p
simplificar la notación se han omitido los superíndices.
p1  2.25
4  54  1.8125
p2  1.25
4  0.75 4  94  2.375
p3  2.25
4  3.75 4  34  0.75
p4   0.75 4  15 4  3.5625
A continuación se presentan los resultados de subsecuentes iteraciones, en forma tabular

p1k p 2k p3k p 4k  1k  2k  3k  4k
K
0 0,000000 0,000000 0,000000 0,000000
1 1,250000 2,250000 -0,750000 3,750000 1,2500 2,2500 -0,7500 3,7500
2 1,812500 2,375000 0,750000 3,562500 0,5625 0,1250 1,5000 -0,1875
3 1,843750 2,890625 0,734375 3,937500 0,0313 0,5156 -0,0156 0,3750
4 1,972656 2,894531 0,957031 3,933594 0,1289 0,0039 0,2227 -0,0039
5 1,973633 2,982422 0,957031 3,989258 0,0010 0,0879 0,0000 0,0557
6 1,995605 2,982666 0,992920 3,989258 0,0220 0,0002 0,0359 0,0000
7 1,995667 2,997131 0,992981 3,998230 0,0001 0,0145 0,0001 0,0090
8 1,999283 2,997162 0,998840 3,998245 0,0036 0,0000 0,0059 0,0000
9 1,999290 2,999531 0,998852 3,999710 0,0000 0,0024 0,0000 0,0015
10 1,999883 2,999536 0,999810 3,999713 0,0006 0,0000 0,0010 0,0000
11 1,999884 2,999923 0,999812 3,999953
12 1,999981 2,999924 0,999969 3,999953
13 1,999981 2,999987 0,999969 3,999992
14 1,999997 2,999988 0,999995 3,999992
15 1,999997 2,999998 0,999995 3,999999
16 1,999999 2,999998 0,999999 3,999999
17 1,999999 3,000000 0,999999 4,000000
18 2,000000 3,000000 1,000000 4,000000
Tabla 3.1 Solución del sistema 3.23 por el método de Jacobi.

b) Método de Gauss-Seidel
Para el cálculo del primer elemento del vector p (1) , se sustituye p ( 0 ) en la
primera ecuación de 3.24, para simplificar la notación se han omitido los su-
períndices.
0 5 5
p1     1.25
4 4 4

Ing. Hermas Herrera Callejas Página: 13 de 32

Programación Aplicada Capítulo 3–Sistemas de Ecuaciones

Para el cálculo de p2 de p (1) , se emplea el valor de p1 ya obtenido (1/4) y los

valores de p2, p3 y p4 de p ( 0 ) . Así
1.25 0 9
p2     2.5625
4 4 4
Con los valores de p1 y p2 ya obtenidos y con p3 y p4 de p(0) se evalúa p3 de p (1) .
2.5625 0 3
p3     0.109375
4 4 4
Finalmente, con los valores de p1, p2 y p3 calculados previamente y con p4 de p(0),
se obtiene la última componente de p (1)
 0.109375 15
p4    3.722656
4 4
Entonces p (1)  [1.25 2.5625  0.109375 3.722656]T
Para la segunda iteración (cálculo de p 2) se procede de igual manera.
2.5625 5
p1    1.890625
4 4
1.890625  0.109375 9
p2     2.6953125
4 4 4
2.6953125 3.722656  3
p3     0.854492
4 4 4
0.854492 15
p4    3.963623
4 4
Con lo que p ( 2)  [1.890625 2.695313 0.854492 3.963623]T
En la tabla 3.2 se presentan los resultados de las iteraciones subsecuentes.

p1k p 2k p3k p 4k  1k  2k  3k  4k
K
0 0,000000 0,000000 0,000000 0,000000
1 1,250000 2,562500 -0,109375 3,722656 1,2500 2,5625 -0,1094 3,7227
2 1,890625 2,695313 0,854492 3,963623 0,6406 0,1328 0,9639 0,2410
3 1,923828 2,944580 0,977051 3,994263 0,0332 0,2493 0,1226 0,0306
4 1,986145 2,990799 0,996265 3,999066 0,0623 0,0462 0,0192 0,0048
5 1,997700 2,998491 0,999389 3,999847 0,0116 0,0077 0,0031 0,0008
6 1,999623 2,999753 0,999900 3,999975 0,0019 0,0013 0,0005 0,0001
7 1,999938 2,999960 0,999984 3,999996 0,0003 0,0002 0,0001 0,0000
8 1,999990 2,999993 0,999997 3,999999
9 1,999998 2,999999 1,000000 4,000000
10 2,000000 3,000000 1,000000 4,000000
Tabla 3.2 Solución del sistema 3.23 por el método de Gauss-Seidel.
En la aplicación de estas dos variantes son válidas las preguntas siguientes:
1. ¿La sucesión de vectores p (1) , p ( 2 ) , p (3) , ... , converge o se aleja del vector
solución p  [ p1 p 2 ... p n ] ?
T

2. ¿Cuando se detendrá el proceso iterativo?

Las respuestas correspondientes, conocidas como criterio de convergencia, se
dan a continuación:
1. Si la sucesión converge a p, cabe esperar que los elementos de p ( k ) se
vayan acercando a los elementos correspondientes de p, es decir, p1K , a p1 , p 2K , a
p 2 , etc., o que se alejen en caso contrario.

Ing. Hermas Herrera Callejas Página: 14 de 32

Programación Aplicada Capítulo 3–Sistemas de Ecuaciones

2. Cuando:
k 1 k 1
a) Los valores absolutos p1  p1 , p 2  p 2 , etc., sean todos menores de un
k k

número pequeño ε cuyo valor será dado por el programador. O bien

b) Si el número de iteraciones ha excedido un máximo predeterminado MAXIT.
Por otro lado, es natural pensar que si la sucesión p ( 0) , p (1) ,…, converge a p,
la distancia de p ( 0 ) a p, de p (1) a p, etc., se va reduciendo. También es cierto que la
distancia entre cada dos vectores consecutivos p ( 0 ) y p (1) , p (1) y p ( 2 ) , etc., se
decrementa conforme el proceso iterativo avanza; esto es, la sucesión de números
reales:
p (1)  p ( 0 ) , p ( 2 )  p (1) ,..., p ( k 1)  p ( k ) 3.25
Convergirá a cero.
Si, por el contrario, esta sucesión de números diverge, entonces puede
pensarse que el proceso diverge. Con esto, un criterio más es
c) Detener el proceso una vez que p ( k 1)  p ( k ) < ε
Al elaborar un programa de cómputo para resolver sistemas de ecuaciones
lineales, generalmente se utilizan los criterios (a), (b) y (c) o la combinación de (a) y
(b), o la de (b) y (c).
Si se observan las columnas de las tablas 3.1 y 3.2, se advertirá que todas son
sucesiones de números convergentes, por lo que ambos métodos convergen a un
vector, presumiblemente la solución del sistema 3.23.
Si se tomara el criterio (a) con ε = 10 2 y el método de Jacobi, ε se satisface
en la sexta iteración de la tabla 3.1; en cambio si ε = 10 3 , se necesitan 10
iteraciones.
Si se toma ε = 10 3 el método de Gauss-Seidel y el criterio (a), se requerirían
sólo seis iteraciones, como puede verse en la tabla 3.2.
Aunque hay ejemplos en los que Jacobi converge y Gauss-Seidel diverge y
viceversa, en general puede esperarse convergencia más rápida por Gauss-Seidel, o
una manifestación más rápida de divergencia. Esto se debe al hecho de ir usando
los valores más recientes de p k 1 que permitirán acercarse o alejarse más
rápidamente de la solución.

3.3.2 RE-ARREGLO DE ECUACIONES.

Para motivar el re-arreglo de ecuaciones, se propone resolver el siguiente

sistema con el método de Gauss-Seidel y con ε = 10 2 aplicado a p
( k 1)
 p (k )
 p1  3 p2  5 p3  2 p4  10
p1  9 p2  8 p3  4 p4  15
3.26
p2  p4  2
2 p1  p2  p3  p4  3
Al resolver para p1 de la primera ecuación, para p 2 de la segunda, p 3 de la
cuarta y p 4 de la tercera se obtiene

Ing. Hermas Herrera Callejas Página: 15 de 32

Programación Aplicada Capítulo 3–Sistemas de Ecuaciones

p1  3 p2  5 p3
 10  2 p4
p 158 4
p2   1   p3  p4
9 99 9
p3   2 p1  p 2 3  p4
p4   p2 2
Con el vector cero como vector inicial, se tiene la siguiente sucesión de vectores.
Nótese que el proceso diverge.

P1k P2k P3k P4k

K DIFERENCIAS
0 0,00000 0,00000 0,00000 0,00000
1 -10,00000 2,77778 14,22222 -0,77778 -10,00000 2,77778 14,22222 -0,77778
2 67,88889 -18,17284 -121,38272 20,17284 77,88889 -20,95062 -135,60494 20,95062
3 -631,08642 170,71742 1108,62826 -168,71742 -698,97531 188,89026 1230,01097 -188,89026
Tabla 3.3. Aplicación del método de Gauss-Seidel al sistema 3.26.
Si el proceso iterativo diverge, como es el caso, un re-arreglo de las
ecuaciones puede originar convergencia; por ejemplo, en lugar de despejar p1 de la
primera ecuación, p 2 de la segunda, etc., cabe despejar las diferentes pi de
diferentes ecuaciones, teniendo cuidado de que los coeficientes de las pi
despejadas sean distintos de cero.
Esta sugerencia presenta, para un sistema de n ecuaciones, n! distintas
formas de re-arreglar dicho sistema. A fin de simplificar este procedimiento, se
utilizará el siguiente teorema

Teorema 3.1 Los procesos de Jacobí y Gauss-Seidel convergirán si en la matriz

coeficiente cada elemento de la diagonal principal es mayor (en valor absoluto) que
la suma de los valores absolutos de todos los demás elementos de la misma fila o
columna (matriz diagonal dominante). Es decir, se asegura h convergencia si:
n
a ii  
j 1
a ij
1 i  n
j i

y 3.27
n
a ii  
j 1
a ji
1 i  n
j i

Este teorema no será de mucha utilidad si se toma al pie de la letra, ya que

contados sistemas de ecuaciones lineales poseen matrices coeficiente
diagonalmente dominantes; sin embargo, si se arreglan las ecuaciones para tener el
sistema lo más cercano posible a las condiciones del teorema, algún beneficio se
puede obtener. Ésta es la pauta para reordenar las ecuaciones y obtener o mejorar
la convergencia, en el mejor de los casos. A continuación se ilustra esto, re-
arreglando el sistema 4.26, despejando p1 de la ecuación 4, p 2 de la ecuación 2,
p 3 de la ecuación 1 y p 4 de la ecuación 3 para llegar a:

Ing. Hermas Herrera Callejas Página: 16 de 32

Programación Aplicada Capítulo 3–Sistemas de Ecuaciones

p2 p3 p4 3
p1     
2 2 2 2
p1 8 p3 4 p4 15
p2     
9 9 9 9
p1 3 p2 2 p4 10
p3    
5 5 5 5
p4   p2 2
Los resultados para las primeras 18 iteraciones con el vector cero como vector
inicial se muestran en la tabla 3.4
Antes de continuar las iteraciones, puede observarse en la tabla 3.4 que los
valores de p (18) parecen converger al vector
p  1 2
T
0 1
Con la sustitución de estos valores en el sistema 3.26, se comprueba que p1
= -1, p 2 = 0.0, p3 = 1 y p 4 = 2 es el vector solución y por razones obvias se detiene
el proceso.

P1k P2k P3k P4k

K DIFERENCIAS
0 0,00000 0,00000 0,00000 0,00000
1 -1,50000 1,83333 0,60000 0,16667 -1,50000 1,83333 0,60000 0,16667
2 -2,63333 1,35185 0,59556 0,64815 -1,13333 -0,48148 -0,00444 0,48148
3 -2,14963 1,08807 0,65798 0,91193 0,48370 -0,26379 0,06242 0,26379
4 -1,91705 0,88950 0,71811 1,11050 0,23258 -0,19856 0,06014 0,19856
5 -1,74856 0,72907 0,76865 1,27093 0,16849 -0,16043 0,05053 0,16043
6 -1,61339 0,59784 0,81025 1,40216 0,13516 -0,13124 0,04160 0,13124
7 -1,50296 0,49026 0,84439 1,50974 0,11044 -0,10758 0,03414 0,10758
8 -1,41245 0,40204 0,87239 1,59796 0,09051 -0,08822 0,02800 0,08822
9 -1,33824 0,32970 0,89535 1,67030 0,07422 -0,07234 0,02296 0,07234
10 -1,27737 0,27038 0,91418 1,72962 0,06086 -0,05932 0,01883 0,05932
11 -1,22747 0,22173 0,92962 1,77827 0,04991 -0,04865 0,01544 0,04865
12 -1,18654 0,18183 0,94229 1,81817 0,04093 -0,03990 0,01266 0,03990
13 -1,15297 0,14911 0,95267 1,85089 0,03356 -0,03272 0,01038 0,03272
14 -1,12545 0,12228 0,96119 1,87772 0,02753 -0,02683 0,00852 0,02683
15 -1,10287 0,10028 0,96817 1,89972 0,02257 -0,02200 0,00698 0,02200
Tabla 3.4. Aplicación del método de Gauss-Seidel al sistema 3.26, re-arreglando las ecuaciones para
obtener una aproximación a un sistema diagonal dominante.
3.4 SISTEMAS DE ECUACIONES NO LINEALES

Se vio cómo encontrar las raíces de una ecuación de la forma f(x) = 0.

Además se vieron técnicas iterativas de solución de sistemas de ecuaciones lineales
Ax = b. Ahora veremos sistemas de ecuaciones no lineales donde se tiene un
sistema de varias ecuaciones con varias incógnitas, cuya representación es:
f1(x1, x2, x3, …, xn) = 0
f2(x1, x2, x3, …, xn) = 0
.
. 3.28
fn(x1, x2, x3, …, xn) = 0
donde fi(x1,x2,x3 …, xn) para 1 ≤ i ≤ n es una función (lineal o no) de las variables

Ing. Hermas Herrera Callejas Página: 17 de 32

Programación Aplicada Capítulo 3–Sistemas de Ecuaciones

independientes x1, x2, x3, …, xn.

Si la ecuación 3.28 consiste sólo en una ecuación de una incógnita (n = 1), se
tiene la ecuación del tipo f(x) = 0. En cambio si n > 1 se tendrá un sistema de
ecuaciones lineales y f1, f2, … fn son todas funciones lineales de x1, x2, x3, …, xn.
Por todo esto, es fácil entender que los métodos iterativos de solución de la
ecuación 3.28 son extensiones de los métodos para ecuaciones no lineales en una
incógnita y emplean las ideas que se aplicaron al desarrollar los algoritmos iterativos
para resolver A x = b.
A continuación se dan algunos ejemplos.
a) 1 1 , x 2 )  x12  x 22  4  0
f ( x
f 2 ( x1 , x 2 )  x 2  x12  0
b) f1 ( x1 , x 2 )  10( x 2  x12 )  0
f 2 ( x1 , x 2 )  1  x1  0
c) f 1 ( x1 , x 2 , x3 )  x1 x 2 x3  10 x13  x 2  0
f 2 ( x1 , x 2 , x3 )  x1  2 x 2 x3  sen( x 2 )  15  0
f 3 ( x1 , x 2 , x3 )  x 22  5 x1 x3  3x33  3  0

3.4.1 DIFICULTAD EN SOLUCIÓN DE SISTEMAS DE ECUACIONES NO LINEALES

Antes de desarrollar los métodos iterativos para resolver sistemas de

ecuaciones no lineales con varias Incógnitas, se destacarán algunas de las
dificultades que se presentan al aplicar estos métodos.
 Es imposible graficar las superficies multidimensionales definidas por las
ecuaciones de los sistemas para n > 2.
 No es fácil encontrar ‘buenos” valores iniciales.

Para atenuar estas dificultades se darán algunas sugerencias aplicables antes

de un intento formal de solución de la ecuación 3.28.

[Link] REDUCCIÓN DE ECUACIONES

Resulta muy útil tratar de reducir analíticamente el número de ecuaciones y de

incógnitas antes de intentar una solución numérica. En particular, trátese de resolver
alguna de las ecuaciones para alguna de las incógnitas. Después, sustitúyase la
ecuación resultante para esa incógnita en todas las demás ecuaciones; con esto el
sistema se reduce en una ecuación y una incógnita. Continúese de esta manera
hasta donde sea posible.
Por ejemplo, en el sistema
f1 ( x1 , x 2 )  10( x 2  x12 )  0
f 2 ( x1 , x 2 )  1  x1  0
se despeja x1 en la segunda ecuación
x1 = 1
y se sustituye en la primera
10(x2 – 12) = 0
cuya solución, x2 = 1, conjuntamente con x1 = 1 proporciona una solución del sistema
dado, sin necesidad de resolver dos ecuaciones con dos incógnitas.

Ing. Hermas Herrera Callejas Página: 18 de 32

Programación Aplicada Capítulo 3–Sistemas de Ecuaciones

[Link] PARTICIÓN DE ECUACIONES

A veces resulta más sencillo dividir las ecuaciones en subsistemas menores y

resolverlos por separado. Considérese por ejemplo el siguiente sistema de cinco
ecuaciones con cinco incógnitas
f 1 ( x1 , x 2 , x3 , x 4 , x5 )  0
f 2 ( x1 , x 2 , x 4 , x5 )  0
f 3 ( x1 , x3 , x 4 , x5 )  0
f 4 ( x2 , x4 )  0
f 5 ( x1 , x 4 )  0
En vez de atacar las cinco ecuaciones al mismo tiempo, se resuelve el
subsistema formado por f 2 , f 4 , f 5 . Las soluciones de este subsistema se utilizan
después para resolver el subsistema compuesto por las ecuaciones f1 y f3
En general, una partición de ecuaciones es la división de un sistema de
ecuaciones en subsistemas llamados bloques. Cada bloque de la partición es el
sistema de ecuaciones más pequeño que incluye todas las variables que es preciso
resolver.

[Link] TANTEO DE ECUACIONES

Supóngase que se quiere resolver el siguiente sistema de cuatro ecuaciones

con cuatro incógnitas
f 1 ( x 2 , x3 )  0
f 2 ( x 2 , x3 , x 4 )  0
f 3 ( x1 , x 2 , x3 , x 4 )  0
f 4 ( x1 , x 2 , x3 )  0
No se pueden dividir en subsistemas, sino que es preciso resolverlas
simultáneamente. Sin embargo, es posible abordar el problema por otro camino.
Supóngase que se estima un valor de x3. Se podría obtener así x2 a partir de f1, x4 de
f2 y x1 de f3. Finalmente, se comprobaría con f 4 la estimación hecha de x 3 Sí f4 fuese
cero o menor en magnitud que un valor predeterminado o criterio de exactitud Є, la
estimación x3 y los valores de x2, x4 y x1 obtenidos con ella, serían una aproximación a
la solución del Sistema dado. En caso contrario, habría que proponer un nuevo valor
de x3 y repetir el proceso.
Nótese la íntima relación que guarda este método con el método de punto fijo,
ya que un problema multidimensional se reduce a unidimensional en x 3:
h(x3) = 0.

[Link] VALORES INICIALES

a) De consideraciones físicas
Sí el sistema de ecuaciones 3.28 tiene un significado físico, con frecuencia es
posible acotar los valores de las incógnitas a partir de consideraciones físicas. Por
ejemplo, si alguna de las variables x i, representa la velocidad de flujo de un fluido,
ésta no podrá ser negativa. Por tanto x i ≥ 0 En el caso de que xi represente una
concentración expresada como fracción peso o fracción molar de una corriente de

Ing. Hermas Herrera Callejas Página: 19 de 32

Programación Aplicada Capítulo 3–Sistemas de Ecuaciones

alimentación, se tiene que 0 ≤ xi ≤ 1.

b) De consideraciones geométricas
En caso de tener un sistema de dos ecuaciones con dos incógnitas
f 1 ( x1 , x 2 )  0
f 2 ( x1 , x 2 )  0
Cada una define, en general, una curva en el plano x 1 – x2 y el problema de
resolver el sistema puede verse como el problema de encontrar el punto o los puntos
de intersección de estas dos curvas. Graficando (puede usarse el software GC, el
Math-CAD, o un programa que grafique) pueden obtenerse buenos valores iniciales.
Sea por ejemplo el sistema
f 1 ( x1 , x 2 )  0
f 2 ( x1 , x 2 )  0

Al graficar f1 y f2 se obtiene la figura 3.1, en donde se podrán apreciar valores

iniciales muy cercanos a la solución.

Figura 3.1 Solución gráfica de un sistema de dos ecuaciones

Por último, resulta muy conveniente conocer bien las características de cada método
de solución del sistema 3.28 para efectuar la elección más adecuada del mismo.
3.4.2 MÉTODO DE PUNTO FIJO MULTIVARIABLE

Los algoritmos discutidos aquí son, en principio, aplicables a sistemas de

cualquier número de ecuaciones. Sin embargo, para ser más concisos y evitar
flotación complicada, se considerará sólo el caso de dos ecuaciones con dos incóg-
nitas. Estas generalmente se escribirán como
f1(x, y) = 0
f2(x, y) = 0 3.29
y se tratará de encontrar pares de valores (x, y) que satisfagan ambas ecuaciones.
Como en el método de punto fijo, se resolverá la primera ecuación para alguna de las
variables, x por ejemplo, y la segunda para y.
x = g1(x, y)
y = g2(x, y) 3.30
Al igual que en los métodos mencionados, se tratará de obtener la estimación

Ing. Hermas Herrera Callejas Página: 20 de 32

Programación Aplicada Capítulo 3–Sistemas de Ecuaciones

(k + 1)-ésima a partir de la estimación k-ésima con la expresión

x k 1  g1 ( x k , y k )
y k 1  g 2 ( x k , y k ) 3.31
Se comienza con valores iniciales x , y , se calculan nuevos valores x , y1 y se
0 0 1

repite el proceso, esperando que después de cada iteración los valores de x k, yk se

aproximen a la raíz buscada x, y , la cual cumple con
x  g1 ( x , y )
y  g 2 (x, y)
Por analogía con los casos discutidos, puede predecirse el comportamiento y
las características de este método de punto fijo multivariable.
Como se sabe, en el caso de una variable la manera particular de pasar de
f(x) = 0 a x = g(x), afecta la convergencia del proceso iterativo. Entonces debe
esperarse que la forma en que se resuelve para x = g 1(x, y) y y = g2(x, y) afecte la
convergencia de las iteraciones 3.31.
Por otro lado, se sabe que el reordenamiento de las ecuaciones en el caso
lineal afecta la convergencia, por lo que puede esperarse que la convergencia del
método en estudio dependa de si se despeja x de f2 o de f1.
Finalmente, la convergencia - en caso de existir - es de primer orden, cabe
esperar que el método iterativo multivariable tenga esta propiedad.

Ejemplo 3.6
Encuentre una solución del sistema de ecuaciones no lineales
f1(x, y) = x2 – 10x + y2 + 8 = 0
f2(x, y) = xy2 + x – 10y + 8 = 0.

SOLUCIÓN
Con el despeje de x del término (-10x) en la primera ecuación y de y del término (-
10y) en la segunda ecuación, resulta
x2  y2  8
x
10
xy  x  8
2
y
10
o con la notación de la ecuación 3.31
(xk )2  ( y k )2  8
x k 1 
10
x k ( y k )2  xk  8
y k 1 
10
Con los valores iniciales x0 = 0, y0 = 0, se inicia el proceso iterativo

Primera iteración
02  02  8
x1   0.8
10
0( 0) 2  0  8
y1   0.8
10

Segunda iteración

Ing. Hermas Herrera Callejas Página: 21 de 32

Programación Aplicada Capítulo 3–Sistemas de Ecuaciones

(0.8) 2  (0.8) 2  8
x2   0.928
10
0.8(0.8) 2  0.8  8
y2   0.9312
10
Al continuar el proceso iterativo, se encuentra la siguiente sucesión de vectores
k xk yk
0 0.00000 0.00000
1 0.80000 0.80000
2 0.92800 0.93120
3 0.97283 0.97327
4 0.98937 0.98944
5 0.99578 0.99579
6 0.99832 0.99832
7 0.99933 0.99933
8 0.99973 0.99973
9 0.99989 0.99989
10 0.99996 0.99996
11 0.99998 0.99998
12 0.99999 0.99999
13 1.00000 1.00000
Para observar la convergencia del proceso iterativo, se pudieron usar criterios
como distancia entre dos vectores consecutivos o bien las distancias componente a
componente de dos vectores consecutivos. También existe un criterio de
convergencia equivalente que dice: Una condición suficiente aunque no necesaria,
para asegurar la convergencia es que
g 1 g 2 g1 g 2
  M 1 ;   M 1 3.32
x x y y
para todos los puntos (x, y) de la región del plano que contiene todos los valores
(xk, yk) y la raíz buscada ( x , y ) .
Por otro lado; si M es muy pequeña en una región de interés, la iteración
converge rápidamente; si M es cercana a 1 en magnitud, entonces la iteración puede
converger lentamente. Este comportamiento es similar al del caso de una función
univariable discutido anteriormente.
Por lo general es muy difícil encontrar el sistema 3.30 a partir de la ecuación
3.29, de modo que satisfaga la condición 3.32.
De todas maneras, cualquiera que sea el sistema 3.29 a que se haya llegado y
que se vaya a resolver con este método, puede aumentarse la velocidad de
convergencia usando desplazamientos sucesivos en lugar de los desplazamientos
simultáneos del esquema 3.31. Es decir, se iteraría mediante
x k 1  g1 ( x k , y k )
y k 1  g 2 ( x k 1 , y k )
Si la iteración por desplazamientos simultáneos diverge, generalmente el
método por desplazamientos sucesivos divergiría más rápido; es decir, se detecta
más rápido la divergencia, por lo que se recomienda en general el uso de
desplazamientos sucesivos en lugar de desplazamientos simultáneos.

Ejemplo 3.7

Ing. Hermas Herrera Callejas Página: 22 de 32

Programación Aplicada Capítulo 3–Sistemas de Ecuaciones

Resuelva el sistema del ejemplo 3.6 utilizando el método de punto fijo multivariable
con desplazamientos sucesivos
f1(x, y) = x2 – 10x + y2 + 8 = 0
f2(x, y) = xy2 + x – 10y + 8 = 0.
Sugerencia: Se pueden seguir los cálculos con un pizarrón electrónico o se
programa una calculadora.

SOLUCIÓN
Al despejar x del término (-10x) y y del término (-10y) de la primera y segunda
ecuaciones, respectivamente, resulta
(x k )2  ( y k )2  8
x k 1  g1 ( x k , y k ) 
10
x ( y k ) 2  x k 1  8
k 1
y k 1  g 2 (x k , y k ) 
10
Al derivar parcialmente, se obtiene
g 1 2 x k g1 2yk
 
x 10 y 10
g 2 ( y k ) 2  1 g 2 2 x k 1 y k
 
x 10 y 10
y evaluadas en x0 = 0 y en y0 = 0
g1 g 2 1
0  para x0 y y0
x x 10
g1 g 2
0 0 para x0 y y0
y y
con lo que se puede aplicar la condición 3.32
g 1 g 2
+ = 0 + 1/10 = 1/10 < 1
x x
g 1 g 2
+ =0+0=0<1
y y
Que se satisface; si los valores sucesivos de la iteración: x 1, y1; x2, y2; x3, y3; …la
satisfacen también, se llega entonces a ( x , y ) .

Primera iteración
02  02  8
x1   0.8
10
0.8(0) 2  0.8  8
y1   0.88
10
Cálculo de la distancia entre el vector inicial y el vector [x 1, y1]T
x (1)  x ( 0 )  (0.8  0.0) 2  (0.88  0.0) 2  1.18929

Segunda iteración
(0.8) 2  (0.88) 2  8
x2   0.94144
10
0.94144(0.88) 2  0.94144  8
y2   0.96704
10
Cálculo de la distancia entre [x2, y2]T y [x1, y1]T

Ing. Hermas Herrera Callejas Página: 23 de 32

Programación Aplicada Capítulo 3–Sistemas de Ecuaciones

x ( 2 )  x (1)  (0.94144  0.8) 2  (0.96704  0.88) 2  0.16608

A continuación se muestran los resultados de las iteraciones
k xk yk |x(k+1)-xk|
0 0.00000 0.00000
1 0.80000 0.88000 1.18929
2 0.94144 0.96705 0.16608
3 0.98215 0.99006 0.04677
4 0.99448 0.99693 0.01411
5 0.99829 0.99905 0.00436
6 0.99947 0.99970 0,00135
7 0.99983 0.99991 0.00042
8 0.99995 0.99997 0.00013
9 0.99998 0.99999 0.00004
10 0.99999 1.00000 0.00001
11 1.00000 1.00000 0.00001
Nótese que se requirieron once iteraciones para llegar al vector solución (1, 1)
contra 13 del ejemplo 3.6, donde se usaron desplazamientos simultáneos.
A continuación se presenta un algoritmo para el método de punto fijo
multivariable en sus versiones de desplazamientos simultáneos y desplazamientos
sucesivos.

ALGORITMO 3.1 Método de punto fijo multivariable

Para encontrar una solución aproximada de un sistema de ecuaciones no lineales

g(x) = x, proporcionar las funciones G(I, x), I = 1, 2, …, N y los
DATOS: El número de ecuaciones N, el vector de valores iniciales x, el
criterio de convergencia EPS, el número máximo de iteraciones
MAXIT y M = 0 para desplazamientos sucesivos o M = 1 para
desplazamientos simultáneos.
RESULTADOS: Una solución aproximada x o mensaje ‘NO HUBO
CONVERGENCIA’.
PASO 1. Hacer K = 1
PASO 2. Mientras K ≤ MAXIT, repetir los pasos 3 a 14.
PASO 3. Si M = 0, hacer xaux = x. De otro modo continuar
PASO 4. Hacer I = 1
PASO 5. Mientras I ≤ N, repetir los pasos 6 y 7.
PASO 6. Si M = 0, hacer X(I) = G(I, x). De otro modo hacer
XAUX(I) = G(I, x)
PASO7. Hacer I = I + 1
PASO 8. Hacer I = 1
PASO 9. Mientras I ≤ N, repetir los pasos 10 y 11.
PASO 10. Si ABS(XAUX(I) - X(I)) > EPS ir al paso 13. De
otro modo continuar.
PASO 11. Hacer I = I + 1
PASO 12. IMPRIMIR x Y TERMINAR.
PASO 13. Si M =1 hacer x = xaux. De otro modo continuar
PASO 14. Hacer K = K + 1
PASO 15. IMPRIMIR mensaje ‘NO HUBO CONVERGENCIA’ y TERMINAR.

Ing. Hermas Herrera Callejas Página: 24 de 32

Programación Aplicada Capítulo 3–Sistemas de Ecuaciones

Sugerencia: Desarrolle este algoritmo con Math-CAD o un software equivalente.

3.4.3 MÉTODO DE NEWTON-RAPHSON MULTIVARIABLE

El método iterativo para sistemas de ecuaciones converge linealmente. Como

en el método de una incógnita, puede crearse un método de convergencia cuadrática
es decir, el método de Newton-Raphson multivariable. A continuación se obtendrá
este procedimiento para dos variables; la extensión a tres o más variables es viable
generalizando los resultados.
Supóngase que se está resolviendo el sistema
f1(x, y) = 0
f2(x, y) = 0
donde ambas funciones son continuas y diferenciables, de modo que puedan
expandirse en serie de Taylor. Esto es
f f 1  2 f 2 f 2 f 
f ( x, y )  f ( a , b )  ( x  a)  ( y  b)   ( x  a)2  2 ( x  a )( y  b)  ( y  b) 2   ...
x y 2!  xx xy yy 
donde f(x,y) se ha expandido alrededor del punto (a,b) y todas las derivadas
parciales están evaluadas en (a ,b)
Expandiendo f1 alrededor de (xk, yk)
f 1 k 1 f
f 1 ( x k 1 , y k 1 )  f 1 ( x k , y k )  ( x  x k )  1 ( y k 1  y k ) 
x y
1   2 f 1 k 1  2 f 1 k 1 k 1  2 f1 k 1 
 ( x  x k 2
)  2 ( x  x k
)( y  y k
)  ( y  y k ) 2   ...
2!  xx xy yy 
3.33
donde todas las derivadas parciales están evaluadas en (x k,yk). De la misma forma
puede expandirse f2 como sigue
f 2 k 1 f
f 2 ( x k 1 , y k 1 )  f 2 ( x k , y k )  ( x  x k )  2 ( y k 1  y k ) 
x y
1   2 f 2 k 1  2 f 2 k 1 k 1  2 f 2 k 1 
 ( x  x k 2
)  2 ( x  x k
)( y  y k
)  ( y  y k ) 2   ...
2!  xx xy yy 
3.34
Al igual que en la ecuación 3.33, todas las derivadas parciales de 3.34 están
evaluadas en (xk, yk)
Ahora supóngase que xk+1 y yk+1 están tan cerca de la raíz buscada ( x , y ) que
los lados izquierdos de las dos últimas ecuaciones son casi cero; además, asúmase
que xk y yk están tan próximos de xk+1 y yk+1 que pueden omitirse los términos a partir
de los que se encuentran agrupados en paréntesis rectangulares. Con esto las ecua-
ciones 3.33 y 3.34 se simplifican a
f 1 k 1 f
0  f1 ( x k , y k )  ( x  x k )  1 ( y k 1  y k )
x y
f f
0  f 2 ( x k , y k )  2 ( x k 1  x k )  2 ( y k 1  y k ) 3.35
x y
Para simplificar aún más, se cambia la notación con
xk+1 – xk = h
yk+1 – yk = j 3.36
y así queda la (k+1)-ésima iteración en términos de la k-ésimas, como se ve a con-

Ing. Hermas Herrera Callejas Página: 25 de 32

Programación Aplicada Capítulo 3–Sistemas de Ecuaciones

tinuación
xk+1 = xk + h
yk+1 = yk + j 3.37
La sustitución de la ecuación 3.36 en la 3.35 y el rearreglo dan como resultado
f 1 f
h  1 j   f1 ( x k , y k )
x y
f 2 f 2
h j   f 2 (x k , y k ) 3.38
x y
El cual es un sistema de ecuaciones lineales en las incógnitas h y j (recuérdese que
las derivadas parciales de la ecuación 3.38, así como f 1 y f2 están evaluadas en (xk,yk)
y, por tanto, son números reales).
Este sistema de ecuaciones lineales resultante tiene solución única, siempre
que el determinante de la matriz de coeficientes o matriz jacobiana J no sea cero; es
decir si
f 1 f 1
x y
J  0
f 2 f 2
x y
Precisando: El método de Newton-Raphson consiste fundamentalmente en
formar y resolver el sistema 3.38, esto último por alguno de los métodos vistos. Con
la solución y la ecuación 3.37 se obtiene la siguiente aproximación.
Este procedimiento se repite hasta satisfacer algún criterio de convergencia
establecido.
Es interesante notar que como en el caso unidimensional, este método puede
obtenerse encontrando un plano tangente a cada f de la ecuación 3.29 en (xk,yk) y
luego encontrar el cero común de estos planos; es decir, hallar un plano tangente en
(xk, yk) tanto a la superficie f 1 como a la superficie f 2, y luego la intersección de cada
plano tangente con el plano x-y, con lo cual se obtienen dos líneas rectas en el plano
x-y y, por último, la intersección de estas dos líneas rectas, que da el cero común de
los planos tangentes.
Cuando converge este método, lo hace con orden dos, y requiere que el
vector inicial (x0, y0) esté muy cerca de la raíz buscada ( x , y )

Ejemplo 3.8
Use el método de Ncwton-Raphson para encontrar una solución aproximada del
sistema
f1(x, y) = x2 – 10x + y2 + 8 = 0
f2(x, y) = xy2 + x – 10y +8 = 0
con el vector inicial: [x0, y0]T = [0, 0]T

SOLUCIÓN
Primero se forma la matriz coeficiente del sistema 3.38, también conocida como
matriz de derivadas parciales
f 1 f 1
 2 x  10  2y
x y
f 2 f 2
 y2 1  2 xy  10
x y

Ing. Hermas Herrera Callejas Página: 26 de 32

Programación Aplicada Capítulo 3–Sistemas de Ecuaciones

que aumentada en el vector de funciones resulta en

2 x  10 2y  x 2  10 x  y 2  8
y2 1 2 xy  10  xy 2  x  10 y  8
Primera iteración
Al evaluar la matriz en [x0, y0]T se obtiene
 10 0 8
1  10  8
que al resolverse por eliminación de Gauss da
h = 0.8, j = 0.88
al sustituir en la ecuación 3.37 se obtiene
x1 = x0 + h = 0 + 0.8 = 0.8
y1 = y0 + j = 0 + 0.88 = 0.88
Cálculo de la distancia entre x(0) y x(1)
x (1)  x ( 0 )  (0.8  0) 2  (0.88  0) 2  1.18929
Segunda iteración
Al evaluar la matriz en [x1, y1]T resulta
 8.4 1.76  1.41440
1.7744  8.592  0.61952
que por eliminación gaussiana da como nuevos resultados de h y j
h = 0.36497, j = 0.1117
de donde
x2 = x1 + h = 0.8 + 0.36497 = 1.16497
y2 = y1 + j = 0.88 + 0.1117 = 0.9917
Cálculo de la distancia entre x(1) y x(2)
x ( 2 )  x (1)  (1.16497  0.8) 2  (0.9917  0.88) 2  0.38168
Con la continuación de este proceso iterativo se obtienen los resultados siguientes
k xk yk |xk+1-xk|
0 0.00000 0.00000
1 0.80000 0.88000 1.18929
2 1.16497 0.99170 0.38168
3 1.31255 1.08099 0.17250
4 0.98255 0.98599 0.00020
Se requirieron cuatro iteraciones para llegar al vector solución (1,1) contra once del
ejemplo 3.7, donde se usó el método de punto fijo con desplazamientos sucesivos.
Sin embargo, esta convergencia cuadrática implica mayor número de cálculos, ya
que -como se puede observar- en cada iteración se requiere
a) La evaluación de 2 x 2 derivadas parciales
b) La evaluación de 2 funciones
c) La solución de un sistema de ecuaciones lineales de orden 2.
Sugerencia: Los cálculos, incluidas las derivadas parciales y la inversa de la matriz,
se pueden ejecutar en Math-CAD o con otro software

[Link] GENERALIZACIÓN

Para un sistema de n ecuaciones no lineales con n incógnitas (véase Ec. 3.28)

y retomando la flotación vectorial y matricial, las ecuaciones 3.38 quedan

Ing. Hermas Herrera Callejas Página: 27 de 32

Programación Aplicada Capítulo 3–Sistemas de Ecuaciones

f 1 f 1 f 1
h1  h2  ...  hn   f1
x1 x 2 x n
f 2 f 2 f 2
h1  h2  ...  hn   f2
x1 x 2 x n
. . .
. . .
. . .
f n f n f n
h1  h2  ...  hn   f n
x1 x 2 x n
3.39
o J h = -f
donde las funciones fi y las derivadas parciales ∂f i/∂xj, i = 1, 2, …, n; j = 1, 2, …, n
están evaluadas en el vector x(k) y
hi  xik 1  xik 1≤i≤n 3.40
De donde
xik 1  xik  hi 1≤i≤n 3.41
( k 1)
o x  x h
(k ) (k )

y la matriz de derivadas parciales (matriz jacobiana), ampliada en el vector de

funciones queda
 f 1 f 1 f 1 
 x ...  f1 
x 2 x n
 1 
 f 2 f 2 f 2  f 2 
...
 x1 x 2 x n 
 . . . 
 . 
 . . . . 
 . . . 
 . 
 f n f n
...
f n
 fn 
 x1 x 2 x n 
3.42
o bien
[ J | f ]
Se presenta a continuación un algoritmo para este método.

ALGORITMO 3.2 Método de Newton-Raphson Multivariable

Para encontrar una solución aproximada de un sistema de ecuaciones no lineales

f(x) = 0, proporcionar la matriz jacobiana ampliada con el vector de funciones (véase
Ec. 3.42) y los
DATOS: El número de ecuaciones N, el vector de valores iniciales x, el número
máximo de iteraciones MAXIT y el criterio de convergencia FF5.
RESULTADOS: El vector solución xn o mensaje ‘NO CONVERGE”.
PASO 1. Hacer K = 1
PASO 2. Mientras K = MAXIT, repetir los pasos 3 a 9.
PASO 3. Evaluar la matriz jacobiana aumentada (3.42).

Ing. Hermas Herrera Callejas Página: 28 de 32

Programación Aplicada Capítulo 3–Sistemas de Ecuaciones

PASO 4. Resolver el sistema lineal (3.39).

PASO 5. Hacer xn = x + h (operación vectorial)
PASO 6. Si |xn – x| > EPS ir al paso 8. De otro modo continuar.
PASO 7. IMPRIMIR xn y TERMINAR.
PASO 8. Hacer x = xn
PASO9. Hacer K = K ± 1
PASO 10. IMPRIMIR “NO CONVERGE” Y TERMINAR.

Ejemplo 3.9
Con el algoritmo 3.2, elabore un programa de propósito general para resolver
sistemas de ecuaciones no lineales. Luego úselo para resolver el Sistema
f1(x1, x2, x3) = 3x1 - cos(x2x3) - 0.5 = 0
f2(x1, x2, x3) = x12 - 625x22 = 0
f3(x1, x2, x3) = e-x1x2 + 20x3 + (10π -3)/3 = 0

SOLUCIÓN
La matriz jacobiana ampliada para el sistema es
 
 3 x3 sen ( x 2 x3 ) x 2 sen( x 2 x3 )  3 x1  cos( x 2 x3 )  0.5 
 2x  1250 x 2 0  x12  625 x 22 
 1 
 x 2 e  x1 x2  x1e  x1 x2 20  x1 x 2 10  3 
e  20 x3 
 3 

Al ejecutar el programa con el vector inicial [1 1 1] T debe producir los siguientes

resultados
k x1 x2 x3 Distancia
0 1.00000 1.00000 1.00000
1 0.90837 0.50065 -0.50286
1.5863
2 0.49927 0.25046 -0.51904
0.47982
3 0.49996 0.12603 -0.52045
0.12444
4 0.49998 0.06460 -0.52199
0.61446E-01
5 0.49998 0.03540 -0.52272
0.29214E-01
6 0.49998 0.02335 -0.52302
0.12052E-01
7 0.49998 0.02024 -0.52309
0.31095E-02
8 0.49998 0.02000 -0.52310
0.23879E-03
9 0.49998 0.02000 -0.52310
0.14280E-05

La solución del sistema es

X1 = 0.49998176
X2 = 0.19999269E-01

Ing. Hermas Herrera Callejas Página: 29 de 32

Programación Aplicada Capítulo 3–Sistemas de Ecuaciones

X3 =.-0.52310085

Nótese que en cada iteración se requiere

a) La evaluación de n2 derivadas parciales
b) La evaluación de n funciones
c) La solución de un sistema de ecuaciones lineales de orden n,
lo que representa una inmensa cantidad de cálculo. Debido a esto, se han elaborado
métodos donde los cálculos no son tan numerosos y cuya convergencia es en
general, superior a la del método de punto fijo (superlineal). A continuación se
presenta el método de Newton-Raphson modificado.

3.4.4 MÉTODO DE NEWTON-RAPHSON MODIFICADO

El método de Newton-Raphson modificado que se describe a continuación consiste

en aplicar el método de Newton-Raphson univariable dos veces (para el caso de un
sistema de n ecuaciones no lineales en n incógnitas, se aplicará n veces), una para
cada variable. Cada que se hace esto, se consideran las otras variables fijas.
Considérese de nuevo el sistema
f1(x, y) = 0
f2(x, y) = 0
Tomando los valores iniciales x 0, y0. se calcula a partir del método de Newton-
Raphson univariable un nuevo valor x1 así
f (x0 , y 0 )
x1  x 0  1
f1 x
∂f1/∂x evaluada en x0, y0
Nótese que se ha obtenido x1 a partir de f1 y los valores más recientes de x,y y: x0,y0
Ahora se usa f2 y los valores más recientes de x,y y (x1, y0) para calcular y1
f ( x1 , y 0 )
y1  y 0  2
f 2 y
donde ∂f2/∂y se evalúa en x1, y0. Se tiene ahora x1 y y1. Con estos valores se calcula
x2, después y2, y así sucesivamente.
Este método converge a menudo si x 0, y0 está muy cerca de ( x , y ) , y requiere
la evaluación de sólo 2n funciones por paso (cuatro para el caso de dos ecuaciones
que se está manejando).
Nótese que se han empleado desplazamientos sucesivos, pero los
desplazamientos simultáneos también son aplicables.

Ejemplo 3.10
Resuelva el sistema
f1(x, y) = x2 – 10x + y2 + 8 = 0
f2(x, y) = xy2 + x – 10y +8 = 0
con el método Newton-Raphson modificado, usando los valores iniciales x 0=0, y0=0.

SOLUCIÓN
Primero se obtiene
f 1 f 2
 2 x  10 y  2 xy  10
x y

Ing. Hermas Herrera Callejas Página: 30 de 32

Programación Aplicada Capítulo 3–Sistemas de Ecuaciones

Primera iteración
Se evalúan f1 y ∂f1/∂x en [0,0]T
f1(0,0) = 8
y
f 1 0
x  10
x
y0
se sustituye
8
x1  0   0.8
 10
Para el cálculo de y1 se necesita evaluar f2 y ∂f2/∂y en x1, y0
F2(0.8, 0) = 0.8(0)2 + 0.8 - 10(0) + 8 = 8.8
f 2 1
x  2(0.8)(0)  10  10
y
y0
se sustituye
8.8
y1  0   0.88
 10

Segunda iteración
f 1 1
f1(0.8, 0.88) = 1.4144 y x  8.4
x
y1
1.4144
x 2  0.8   0.96838
 8.4
Ahora se evalúan f2 y ∂f2/∂y en (x2, y1):
f 2 2
F2(0.96838, 0.88) = 0.91829 y x  8.29565
y
y1
de donde
0.91829
y 2  0.88   0.99070
 8.29565
Continuar las iteraciones y calcular las distancias entre cada dos vectores
consecutivos. Continuar hasta que xk ≈ 1 y yk ≈ 1. Comparar además la velocidad de
convergencia de este método con la velocidad de convergencia del método de
Newton-Raphson y el de punto fijo para este sistema particular.
En la aplicación de este método se pudo tomar f 2 para evaluar x1 y f1 a fin de
evaluar y1, asÍ
f2 (x0 , y 0 )
x x 
1 0

f 2 x
f1 ( x1 , y 0 )
y1  y 0 
f 1 y
Esto puede producir convergencia en alguno de los arreglos y divergencia en
el otro. Es posible saber de antemano si la primera o la segunda forma convergirán
para el caso de sistemas de dos ecuaciones, pero cuando n ≥ 3 las posibilidades son

Ing. Hermas Herrera Callejas Página: 31 de 32

Programación Aplicada Capítulo 3–Sistemas de Ecuaciones

varias (n!) y es imposible conocer cuál de estos arreglos tiene viabilidad de

convergencia, por lo cual la elección se convierte en un proceso aleatorio. Esta
aleatoriedad es la mayor desventaja de este método.
En general, para un sistema de n ecuaciones en n incógnitas: x 1, x2, …, xn, el
algoritmo toma la forma:
k 1 f i ( x1k 1 , x 2k 1 ,..., xik11 , xik ,..., x nk )
xi  xi k
f i 1≤i≤n 3.43
( x1k 1 , x 2k 1 ,..., xik11 , xik ,..., x nk )
xi

ALGORITMO 3.3 Método de Newton-Raphson modificado

Para encontrar una solución aproximada de un sistema de ecuaciones no lineales

f(x) = 0, proporcionar las funciones F(I, x) y las derivadas parciales D(I, x) y los

DATOS: El número de ecuaciones N, el vector de valores iniciales x, el

número máximo de iteraciones MAXIT, el criterio de convergencia
EPS y M = 0 para desplazamientos sucesivos o M = 1 para
desplazamientos simultáneos.
RESULTADOS: El vector solución xn o mensaje “NO CONVERGE”.
PASO1. Hacer K= 1
PASO2. Mientras K ≤ MAXIT, repetir los pasos 3 a 11.
PASO3. Si M = 0 hacer xaux = x (operaciones vectoriales)
PASO4. Hacer I = 1
PASO5. Mientras 1 ≤ N, repetir los pasos 6 y 7.
PASO6. Si M = 0 hacer X(I) = X(I)-F(I,x)/D(I, x), de otro mo-
do hacer XAUX(l) = X(I) - F(I, x)/D(I, x)
PASO7. Hacer I = I + 1
PASO8. Si | xaux – x | > EPS ir al paso 10. De otro modo continuar.
PASO9. IMPRIMIR x y TERMINAR.
PASO10. Si M = 1 hacer x = xaux
PASO11. Hacer K = K + 1
PASO12. IMPRIMIR “NO CONVERGE” Y TERMINAR

Ing. Hermas Herrera Callejas Página: 32 de 32