Apuntes de C´ Calculo a´ lculo Num´ Numerico e´ rico (M´etodos (M e´ todos Matem´ Matematicos) a´ ticos)
A NTONIO RODRIGUEZ L U I S S AINZ DE LOS T ERREROS
Calculo a´ lculo Num´ Numerico e´ rico
Contenido
Prim
Ult
Ant
Sig
Volver
Cerrar
Abandonar
Contenido 2. Calculo a´ lculo matricial 3 2.1. Sistem Sistemas as de ecua ecuacione cioness lineal lineales es . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 2.1.1. Metodos e´ todos exactos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 2.1.2. M´etodos etodos aproximados . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 2.2.. Aut 2.2 Autov ovalo alores res . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 2.2.1. Metodo e´ todo de la potencia . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
Metodos e´ todos Matem´ Matematicos a´ ticos
Contenido
Prim
Ultim
Ant
Sig
Volver
Cerrar
Abandonar
Contenido 2. Calculo a´ lculo matricial 3 2.1. Sistem Sistemas as de ecua ecuacione cioness lineal lineales es . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 2.1.1. Metodos e´ todos exactos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 2.1.2. M´etodos etodos aproximados . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 2.2.. Aut 2.2 Autov ovalo alores res . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 2.2.1. Metodo e´ todo de la potencia . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
Metodos e´ todos Matem´ Matematicos a´ ticos
Contenido
Prim
Ultim
Ant
Sig
Volver
Cerrar
Abandonar
C A P´I TULO 2
Calculo a´ lculo matricial
Calculo a´ lculo Num´ Numerico e´ rico
Contenido
Prim
Ult
Ant
Sig
Volver
Cerrar
Abandonar
2.1.
Sistemas de ecuaciones lineales
Los u´ nicos algoritmos exactos que emplearemos durante este curso ser a´ n algunos de los que se usan para resolver sistemas de ecuaciones lineales y que estudiaremos en la secci´on 2.1.1. Por ser exactos estar´an exentos de error de truncamiento, aunque sujetos al error de redondeo al ser implementados en un ordenador. Tambi´en existen m´etodos aproximados de resolucio´ n de sistemas lineales, que pueden ser m´as eficaces que los primeros en algunos casos, aunque conlleven errores tanto de redondeo como de truncamiento. Los estudiaremos en la secci´on 2.1.2. En cualquiera de los casos, nos restringiremos al estudio de sistemas de ecuaciones lineales compatibles determinados, de mayor inter´e s en las aplicaciones, es decir, sistemas de la forma (2.1) Ax = b con A de tama˜no n
2.1.1.
× n y det A = 0.
M´etodos exactos
Algoritmo de eliminacio´ n de Gauss El algoritmo de eliminaci´on de Gauss es de sobra ´ conocido del curso de Algebra Lineal. Aunque los c´alculos podr´ıan realizarse manualM´etodos Matem´aticos
Contenido
Prim
Ultim
Ant
Sig
Volver
Cerrar
Abandonar
mente para sistemas de tama˜no grande se hace imprescindible el uso del ordenador. Como veremos, el error de redondeo que ello conlleva nos obligar´a a introducir modificaciones en el algoritmo. En su versi´o n m´as simple, el algoritmo de eliminaci´on de Gauss considera la matriz ampliada (A b) del sistema (2.1)
|
(1) 11 (1) 21
(1) a12 (1) a22
(1)
(1)
a a
.. .
.. . .. .
(1) a1n (1) a2n
(1) b1 (1) b2
(1)
(1)
.. .
.. .
an1 an2 . . . ann bn
(2.2)
(donde el super´ındice (1) indica que nos encontramos en la primera fase del algoritmo) y, (1) suponiendo que sea a11 = 0, tomamos ese elemento como pivote y “hacemos ceros”por debajo de e´ l. Para ello, modificaremos las filas de la matriz en la forma
(1)
F i
→ F − i
ai1
(1) a11
F 1 ,
i = 2, . . . , n
donde llamaremos multiplicadores a los factores mi1 C´alculo Num´erico
Contenido
Prim
Ult
(2.3)
(1) (1) i1 /a11
≡a Ant
Sig
(1)
En el caso en que a11 =
Volver
Cerrar
Abandonar
(1)
0, al no estar definidos los multiplicadores, se busca la primera fila para la cual a p1 = 0 y se hace el cambio F 1 F p (a este intercambio de filas se le conoce como pivoteo), para
↔
despu´es realizar las operaciones (2.3). Con todo esto, la matriz del sistema se transforma en (1) 11
a 0 .. .
0
(1) a12 (2) a22 (2)
.. . .. .
(1) a1n (2) a2n
(1) b1 (2) b2
(2)
(2)
.. .
.. .
an2 . . . ann bn
(2.4)
donde todos los elementos a partir de la segunda fila han cambiado (de ah´ı el super´ındice (2) (2) ). En la siguiente etapa, se toma a22 como pivote en la segunda columna, o se pivota (2) si a22 = 0 (la condici´on det A = 0 garantiza que podemos encontrar pivotes en todas las columnas de la matriz), para despu´es realizar las operaciones F i F i mi2 F 2 , con (2) (2) mi2 = ai2 /a22 , para i = 3, . . . , n, consiguiendo de esta forma ceros por debajo del segundo pivote.
El proceso se repite n M´etodos Matem´aticos
→ −
− 1 veces hasta llegar a la matriz triangular superior Contenido
Prim
Ultim
Ant
Sig
Volver
Cerrar
Abandonar
(1) 11
(1) a12 (2) a22
a 0 .. .
0
0
... ... ..
(1) a1n (2) a2n
(1) b1 (2) b2
(n )
(n)
.. .
.
. . . ann bn
(2.5)
que se corresponde con el sistema de ecuaciones (1) a11 x1
+
(1) a12 x2 (2) a22 x2
+ +
··· + ··· +
(1) a1 n−1 xn−1 (2) a2 n−1 xn−1
(n−1) xn−1 −1 n−1
an
+ +
(1) a1n xn (2) a2n xn
(n−1)
= = .. .
(1) b1 (2) b2 (n−1)
+ an 1n xn = bn (n ) (n) ann xn = bn −
(2.6)
equivalente al sistema inicial (de igual soluci´on) y que puede ser resuelto por sustituci´ on regresiva (segunda etapa del m´etodo de eliminaci´on de Gauss despu´es de la triangula(n) (n) ci´o n): de la u´ ltima ecuaci´o n de (2.6) se despeja xn = bn /ann , que se introduce en la pen´ultima ecuaci´on para obtener xn 1 y as´ı sucesivamente. −
C´alculo Num´erico
Contenido
Prim
Ult
Ant
Sig
Volver
Cerrar
Abandonar
Se puede comprobar que el n u´ mero de operaciones aritm´eticas elementales que se realizan siguiendo el m´etodo de eliminaci´on de Gauss en un sistema n n es O(n3 ).1 ´ Este m´etodo para resolver sistemas de ecuaciones es mucho m a´ s efectivo que otros como el m´e todo de Kramer o el c´alculo de la matriz inversa (para hacer x0 = A 1 b) que involucran un nu´ mero mucho mayor de operaciones.
×
−
Pivoteo parcial Sin embargo, el algoritmo de eliminaci´on de Gauss expuesto anteriormente no tiene en cuenta los errores de redondeo que se producen al realizar ciertas operaciones como la divisio´ n entre n´umeros peque˜nos. En particular, al calcular los multipli( j ) ( j ) cadores mij = aij /a jj , se produce un error de redondeo considerable, que se arrastra en ( j )
las operaciones siguientes, si el pivote a jj
0. Ve´amoslo en el siguiente ejemplo.
Ejemplo. Consid´erese el siguiente sistema de ecuaciones:
0 003x + 59 14x = 59 17 . 1 5.291x1
.
−
. 6.130x2 = 46.78 2
(2.7)
´ cuya u´ nica soluci´on es x1 = 10, x2 = 1. Esta puede obtenerse realizando las operaciones del m´etodo de eliminaci´on de Gauss considerando el elemento a11 = 0.003 como pivote 1
Se dice que una funci´on f (n) es del orden de n p cuando n C y n0 tales que f (n) Cn p n n0 .
→ ∞ y se escribe f (n) = O(n ) si existen
M´etodos Matem´aticos
Ant
|
|
∀
Contenido
Prim
Ultim
p
Sig
Volver
Cerrar
Abandonar
y haciendo el cambio F 2 F 2 ampliada del sistema resulta
→ −m
21 F 1 ,
0 003 .
−
donde m21 = 5.291/0.003 = 1763 .¯ 6, la matriz
59.14 104309 .37¯6
59 17 .
(2.8)
−104309.37¯6
de donde por sustitucio´ n regresiva se obtiene la soluci o´ n exacta x1 = 10, x2 = 1. Supongamos ahora que el ordenador trabajara s o´ lo con 4 cifras significativas. De esta forma, el multiplicador se aproximar´ıa como m21 1763. Tras multiplicarlo por la fila primera y restar el resultado a la segunda en aritm´etica de 4 d´ıgitos la nueva matriz ampliada resulta
0 003 .
−
59.14 104300
59 17 .
(2.9)
−104400
de donde, por sustitucio´ n regresiva y en aritm´e tica de 4 d´ıgitos se obtiene la soluci´on err´onea x1 = 10, x2 = 1.001
−
completamente alejada de la soluci´on exacta.
Para evitar este problema se utiliza el llamado m´etodo de pivoteo parcial, cuya finalidad es escoger pivotes lo m´as grandes posible. Para ello, en la fase k -´esima del proceso de C´alculo Num´erico
Contenido
Prim
Ult
Ant
Sig
Volver
Cerrar
Abandonar
(k )
eliminaci´on de Gauss, en lugar de escoger como pivote el elemento akk , se compara e´ ste previamente con el resto de los elementos situados por debajo de e´ l en su misma columna y se escoge el de mayor magnitud, es decir, se halla p k tal que (k )
(k )
|a | = m´ax |a | pk
y despu´es se hace el cambio F k
k in
ik
↔ F . p
Ejemplo. Consideremos de nuevo el sistema del ejemplo 2.1.1 y apliqu´emosle el m´etodo de eliminaci´on de Gauss con pivoteo parcial. Al escoger el pivote para la primera columna previamente comparamos a11 = 0.003 < 5.291 = a21 por lo que hacemos el cambio F 1
↔ F obteniendo por tanto la matriz ampliada 5.291 −6.130 46.78 2
0.003
59.14 59.17
(2.10)
que, utilizando el nuevo multiplicador m21 = 0.000567 , y en aritm´etica de 4 d´ıgitos, se transforma en 5.291 6.130 46.78 (2.11) 59.14 59.14
M´etodos Matem´aticos
Contenido
−
Prim
Ultim
Ant
Sig
Volver
Cerrar
Abandonar
de donde se obtiene la solucio´ n exacta x1 = 10, x2 = 1.
Adem´a s de este m´etodo existe el llamado pivoteo total en el que no so´ lo se intercambian filas sino tambi´en columnas (y por tanto inc´ognitas) para encontrar el pivote de mayor tama˜no, aunque ello obliga a guardar memoria de los intercambios realizados entre las inc´ognitas para despu´es deshacerlos.
Pivoteo parcial escalado Sin embargo, tambi´en el pivoteo parcial nos puede llevar a errores dr´asticos, como veremos en el siguiente ejemplo.
Ejemplo. Consideremos nuevamente el sistema (2.7) y multipliquemos su primera ecuaci´on por 104 . El nuevo sistema
30 00x + 591400x = 591700 .
1
2
5.291x1
− 6.130x
2
(2.12)
= 46.78
tendr´a la misma soluci´on que el antiguo. Sin embargo, al aplicar el pivoteo parcial, al ser a11 = 30 > 5.291 = a21 se utilizar´a el elemento a11 como pivote junto con el multi0.1764, lo que, en aritm´etica de 4 d´ıgitos nos conduce plicador m21 = 5.291/30.00 finalmente a la matriz 30.00 591400 591700 (2.13)
C´alculo Num´erico
Contenido
−104300 −104400 Prim
Ult
Ant
Sig
Volver
Cerrar
Abandonar
que nos lleva nuevamente al resultado err´oneo x1 =
−10, x
2
= 1.001.
Para resolver este problema el m´etodo de pivoteo parcial escalado escala previamente todas las filas de la matriz asociada dividiendo todos sus elementos por el de mayor magnitud antes de realizar la comparaci´on —que ahora se efectu´ a entre magnitudes relativas— para escoger el pivote. Es decir, para cada fila se define previamente
si = m´ax aij j =1,...,n
| |
de modo que, para determinar por ejemplo el pivote correspondiente a la primera columna se tomar´a
|a | = p1
s p
m´ax
i=1,...,n
|a | i1
si
F p . El escalado se utiliza so´ lo para comparar. El pivote a utilizar y se har´a el cambio F 1 ser´a a p1 y no a p1 /s p para evitar el error de redondeo producido en la divisio´ n.
↔
Ejemplo. Si finalmente utilizamos pivoteo parcial escalado para resolver el sistema (2.12), primero calculamos
|a | = 11
s1
M´etodos Matem´aticos
30.00 = 0.5073 591400 Contenido
× 10
Prim
4
−
<
Ultim
|a | = 5.291 = 0.8631 21
s2
Ant
6.130
Sig
Volver
Cerrar
Abandonar
lo que nos lleva a realizar el cambio F 1
↔ F y a considerar el sistema 5.291x − 6.130x = 46.78 2
1
2
(2.14)
30.00x1 + 591400x2 = 591700
que, tras resolver en aritm´etica de 4 d´ıgitos nos conduce al resultado correcto x1 = 10, x2 = 1.
Sistemas mal condicionados Por u´ ltimo, conviene llamar la atenci´on sobre un tipo de problema que se da con frecuencia en el C a´ lculo Num´erico y en particular en la resoluci o´ n de sistemas de ecuaciones. Se dice que un problema est´a mal condicionado cuando una peque˜na modificaci´on en los datos produce variaciones dr´asticas en los resultados. Esto es as´ı para ciertos sistemas de ecuaciones, como veremos en el siguiente ejemplo.
Ejemplo. Consid´erense los sistemas de ecuaciones
S 1
≡
x1 + 2x2 = 10 , 1.1x1 + 2x2 = 10.4
S 2
≡
x1 + 2x2 = 10 1.05x1 + 2x2 = 10.4
(2.15)
cuyas soluciones son x1 = (4, 3)t para S 1 y x2 = (8, 1)t para S 2 . Se observa que, aunque la matriz del sistema se ha modificado s´olo ligeramente, la soluci´on ha cambiado mucho. C´alculo Num´erico
Contenido
Prim
Ult
Ant
Sig
Volver
Cerrar
Abandonar
Por otro lado, si introducimos la solucio´ n de S 2 en S 1 obtenemos
8 + 2 1 = 10 1.1 8 + 2 1 = 10.8 10.4
·
·
·
(2.16)
con lo que valores muy alejados de la solucio´ n satisfacen de manera aproximada el siste ma S 1 . Algo similar ocurre al introducir la soluci´on de S 1 en S 2 Podemos, por tanto, interpretar los sistemas mal condicionados como aquellos que son satisfechos de forma aproximada por un rango amplio de soluciones. Aunque no entraremos en detalle en este tipo de sistemas, es claro que los sistemas del ejemplo anterior estaban muy pr´oximos a ser compatibles indeterminados (el determinante de la matriz asociada en ambos sistemas es pr´oximo a 0). Ya que al representar gr´aficamente cada una de las ecuaciones en el plano se obtienen rectas casi paralelas, es dif´ıcil encontrar su punto de interseccio´ n (que se corresponde con la soluci´on del sistema).
Descomposici´on LU La descomposici´on LU de una matriz es una estrategia para aprovechar las operaciones realizadas sobre la matriz asociada de un sistema Ax = b por el m´etodo de eliminaci´on de Gauss y utilizarla para resolver distintos sistemas con igual matriz asociada y diferentes entradas b. M´etodos Matem´aticos
Contenido
Prim
Ultim
Ant
Sig
Volver
Cerrar
Abandonar
´ Como se recordar´a del curso de Algebra Lineal las operaciones realizadas en la matriz A del sistema durante el m´etodo de eliminaci´on de Gauss pueden implementarse mediante ´ el producto matricial por las llamadas matrices elementales. Estas provienen de realizar sobre la matriz identidad las mismas operaciones que se desean realizar sobre A. En particular, para conseguir ceros del debajo del pivote en el paso k -´esimo del m´etodo de eliminaci´on de Gauss, bastar´a con premultiplicar la matriz A(k) del sistema en dicho paso por la matriz
M (k)
1 0 =
0
··· ··· · ·· ···
..
.
..
.
.. . .. . .. . .. .
..
.
..
.
0
· · · −m
−m
k +1,k
..
.
..
.
..
.
0
..
.
..
.
.. .
..
.
..
.
0
·· ·
0
.. . .. .
n,k
0
0 .. . .. . .. . .. .
(2.17)
1
que continene en la columna k por debajo del pivote los multiplicadores usados para hacer ceros por debajo del mismo en las ecuaciones posteriores. Despu´es de multiplicar la matriz A por la sucesio´ n de matrices M (1) , M (2) , . . . , M ( n 1) , obtenemos la matriz −
C´alculo Num´erico
Contenido
Prim
Ult
Ant
Sig
Volver
Cerrar
Abandonar
asociada del sistema en forma triangular (ver (2.5)), a la que llamaremos U : (1) 11
M (n
1)
(2)
··· M
−
M (1)
a 0 A= .. .
(1)
(1)
≡ U
a12 . . . a1n (2) (2) a22 . . . a2n ..
0
0
. (n)
. . . ann
(2.18)
Podemos despejar A en la ecuaci´on (2.18) multiplicando por la izquierda a ambos lados de la ecuaci´on por (M (n 1) M (2) M (1) ) 1 = (M (1) ) 1 (M (2) ) 1 (M (n 1) ) 1 . Es f a´ cil comprobar que (M (k) ) 1 es como la matriz dada en (2.17) pero cambiando el signo a los multiplicadores, por lo que finalmente obtenemos −
−
·· ·
−
−
1 m =
21
(M (1) ) 1 (M (2) ) −
1
−
(n−1)
··· (M
)
1
−
.. . .. .
0 1 .. . .. .
···
−
··· ·· · ·· · 0 ..
.
..
mn1 mn2
de modo que finalmente llegamos a la factorizaci´on
···
mn
. 1n
−
−
0 0 .. . .. .
1
≡ L
A = LU M´etodos Matem´aticos
Contenido
Prim
Ultim
−
(2.19)
(2.20) Ant
Sig
Volver
Cerrar
Abandonar
donde L es triangular inferior con unos en la diagonal y U es triangular superior y vienen dadas en (2.19) y (2.18) respectivamente. De esta forma, podemos reexpresar el sistema de ecuaciones segu´ n:
Ax = LU x = b
(2.21)
U x = y
(2.22)
Ly = b
(2.23)
de modo que, haciendo el cambio obtenemos el sistema Una vez resuelto el sistema (2.23) para y, introducimos el resultado en el sistema (2.22), cuya solucio´ n es x. De esta forma, hemos trasladado la dificultad de la resoluci´on del sistema (2.21) a la resolucio´ n de los sistemas (2.23) y (2.22), m´as sencillos, ya que sus matrices asociadas son triangulares. El sistema (2.23) se resolver´ıa por sustitucio´ n “hacia delante” y el sistema (2.22) por sustitucio´ n “hacia atr´as”(regresiva). El nu´ mero de operaciones a realizar en cada uno de ellos es O(n2 ), por lo que este m´etodo es venta joso frente al de eliminaci´on de Gauss (O(n3 )). Ahora bien, previamente hay que haber obtenido la descomposici´on A = LU , la cual es equivalente a la eliminaci´on gaussiana. ¿D´onde est´a entonces la ventaja? C´alculo Num´erico
Contenido
Prim
Ult
Ant
Sig
Volver
Cerrar
Abandonar
Durante el m´etodo de eliminaci´on de Gauss el ordenador puede aprovechar las posiciones de memoria guardadas para los elementos de A por debajo de la diagonal —que se van a convertir en ceros durante el proceso de eliminaci o´ n— y almacenar en ellas los respectivos multiplicadores. Una vez calculadas L y U una sola vez, pueden ser utilizadas para resolver sistemas Ax = b para diferentes entradas b.
Ejemplo. Resolver el siguiente sistema de ecuaciones realizando previamente la factorizaci´on LU de su matriz asociada.
x 2x 3x
x2 + 3x4 x2 x3 + x4 x2 x3 + 2x4 + 2x2 + 3x3 x4 + +
1 1
− −
−
1
−x
1
−
= = = =
−
4 1 3 4
(1)
Considerando el elemento a11 = 1 de A(1) como pivote y haciendo los cambios F 2 F 2 2F 1 , F 3 F 3 3F 1 y F 4 F 4 + F 1 obtenemos A(2) :
−
→ −
A(1)
=
M´etodos Matem´aticos
−
→
1 2 3 1
1 1 1 2
0 1 1 3
− − −
Contenido
3 1 ; 2
A(2)
−1
Prim
1 0 = 0 0
Ultim
Ant
Sig
1 1 4 3
0 1 1 3
− − − −
Volver
→
3 −5 −7 2
Cerrar
Abandonar
(2)
Si despu´es tomamos a22 = 1 como pivote y hacemos los cambios F 3 F 3 4F 2 , (3) F 4 + 3F 1 obtenemos A(3) . El tercer pivote ser´ıa a33 = 3, pero no es preciso hacer ceros (3) con e´ l ya que directamente hemos obtenido a43 = 0.
−
A(3) = A(4)
→
1 0 = 0
1 1 0 0
0 1 3 0
3 −5 ≡ U 13
− −
0
−13
Los multiplicadores han sido m21 = 2, m31 = 3, m41 = m43 = 0, por lo que obtenemos
L= C´alculo Num´erico
Contenido
1 2 3 1
− − Prim
0 1 4 3 Ult
−
−1, m
32
= 4, m42 =
−3 y
0 0 0
0 0 1 0 1
Ant
Sig
Volver
Cerrar
Abandonar
El sistema Ly = b resulta
y 2y + y 3y + 4y + y 1
1
2
1
2
3
= = = =
8 7 14 7
−y − 3y + y − cuya soluci´on es (y , y , y , y ) = (8, −9, 26, −26), que introducimos como t´ermino in1
1
2
3
2
4
4
dependiente del sistema U x = y
x + x + 3x −x − x − 5x 3x + 13x 1
2
4
2
3
4
3
4
−13x cuya soluci´on es (x , x , x , x ) = (2, 0, −1, 3).
4
1
2
3
= = = =
8 9 26 26
− −
4
Ahora bien, no toda matriz A admite una descomposici´on LU . En el c´alculo de las matrices L y U expuesto arriba hemos supuesto que las u´ nicas operaciones elementales realizadas sobre la matriz de coeficientes del sistema han sido manipulaciones con las filas del tipo F i F i mik F k . Ahora bien, cuando durante el proceso de eliminaci o´ n
→ −
M´etodos Matem´aticos
Contenido
Prim
Ultim
Ant
Sig
Volver
Cerrar
Abandonar
(k )
de Gauss se anula alg´un elemento de la diagonal, akk = 0, es preciso pivotar, es decir, intercambiar filas, operaci´on que no hemos contemplado al obtener la expresi´on para L. S´olamente aquellas matrices regulares para las que no es necesario pivotar durante el m´etodo de eliminaci´on de Gauss admiten una factorizaci´on LU . Tendremos que buscar condiciones bajo las cuales una matriz admita dicha descomposici´on. Para ello veamos previamente la siguiente definici´on. Definici´ on. matriz diagonal dominante. Se dice que una matriz cuadrada A = (aij ) de orden n es estrictamente diagonal dominante si:
|a | > ii
|
aij ;
j =i
|
i = 1, . . . , n
Es decir, una matriz estrictamente diagonal dominante es aquella en la que en cada fila el valor absoluto del elemento correspondiente a la diagonal supera a la suma de los valores absolutos del resto de los elementos de la fila. Estas matrices poseen propiedades muy importantes, entre ellas, las que se enuncian en el siguiente teorema.
Teorema 2.1.1. Si una matriz A es estrictamente diagonal dominante, entonces C´alculo Num´erico
Contenido
Prim
Ult
Ant
Sig
Volver
Cerrar
Abandonar
A es no singular. Se pueden aplicar en A las operaciones del m etodo de eliminacion ´ ´ de Gauss sin intercambio de filas.
Hemos encontrado, pues, una condicio´ n suficiente para que una matriz A admita una factorizaci´on LU . Basta con que sea estrictamente diagonal dominante.
2.1.2.
M´etodos aproximados
Veremos ahora un par de m´etodos iterativos aproximados de resoluci´on de sistemas de ecuaciones lineales que, por tanto, s´ı est´an sujetos a errores de truncamiento. Ambos son eficientes para sistemas grandes con matrices con un elevado n´umero de ceros. Dichos sistemas aparecen con frecuencia en problemas de resoluci´on de ecuaciones diferenciales en derivadas parciales. M´etodos Matem´aticos
Contenido
Prim
Ultim
Ant
Sig
Volver
Cerrar
Abandonar
M´etodo iterativo de Jacobi
Consid´erese el sistema de ecuaciones
a a
11 x1
+ a12 x2 + 21 x1 + a22 x2 +
·· · + a ·· · + a
1n xn
= b1 2n xn = b2
an1 x1 + an2 x2 +
·· · + a
(2.24)
.. .
x = bn
nn n
donde, supuesto que aii = 0 i podemos despejar de cada ecuaci´on la inco´ gnita correspondiente
∀
x = (b − a x = (b − a 1
1
12 x2
2
2
22 x1
.. .
xn = (bn
−a
x
n1 1
1n xn )/a11
−···−a −···−a
2n xn )/a22
−···−a
(2.25)
x )/ann
nn n
De esta forma, al estar las inc o´ gnitas despejadas, el problema parecer´ıa estar resuelto, aunque no es as´ı, ya que las inc´ognitas tambi´en est´an presentes a la derecha del igual en cada ecuaci´on en (2.25). Sin embargo, podr´ıamos introducir una primera soluci´on apro(0) (0) (0) ximada x(0) (x1 , x2 , . . . , xn )t a la derecha del igual en (2.25) que producir´ıa una (1) (1) (1) soluci´on x(1) (x1 , x2 , . . . , xn )t a la izquierda. Posteriormente, podr´ıamos introducir
≡ ≡
C´alculo Num´erico
Contenido
Prim
Ult
Ant
Sig
Volver
Cerrar
Abandonar
x(1) como entrada para obtener x(2) , y as´ı sucesivamente, produciendo el proceso iterativo (k +1) 1 (k +1) 2
= (b1 = (b2
−a −a
(k +1)
= (bn
−a
x x
.. .
xn
(k ) 12 x2 (k ) 21 x1
−···−a −···−a
(k ) n1 1
−···−a
x
(k ) 1n xn )/a11 (k ) 2n xn )/a22 (k ) n,n−1 n−1
x
(2.26)
)/ann
Podemos reexpresar matricialmente el proceso (2.26) reescribiendo la matriz del sistema como A = T I + D + T S , donde D = diag(a11 , a22 , . . . , ann ), y
0 a T = I
..
.
.. .
..
.
..
an1
···
an n
21
0
··· ·· ·
.. . .. .
. 1
−
0
,
0 T = .. . .. .
S
0
a12 ..
.
· ·· ..
.
..
.
··· ···
a1n .. .
an
1n
−
0
De esta forma, la ecuaci´on (2.26) podr´ıa expresarse como x(k+1) = D T S )x(k) ) que, haciendo T J D 1 (T I + T S ) y cJ D 1 b resulta
≡−
≡
−
1
−
(b
−
x(k+1) = T J x(k) + cJ M´etodos Matem´aticos
Contenido
Prim
Ultim
Ant
− (T
I
+
(2.27) Sig
Volver
Cerrar
Abandonar
El proceso iterativo formalmente definido en (2.27) es el que utilizaremos para obtener sucesivas soluciones aproximadas del sistema (2.24). Estamos dando por supuesto que la sucesi´on [x(n) ]n=0 es convergente y que converge a la soluci´on x de (2.24).2 M´as adelante veremos condiciones suficientes que garantizan cu´ando esto es as´ı. Una vez supuesta la convergencia, hemos de introducir un criterio de paro para el ´ algoritmo (2.27). Este puede consistir en fijar una toleracia TOL al comienzo del problema y parar las iteraciones cuando el error relativo est´e por debajo de ella, es decir, parar en un N tal que x(N ) x(N 1) < TOL ( ) N x ∞
||
||
−
||
−
||
Ejemplo. Aplicar el m´etodo de Jacobi para resolver el sistema:
10x −x 2x
1 1 1
− + −
x2 + 2x3 11x2 x3 + 3x4 x2 + 10x3 x4 3x2 x3 + 8x4
− −
Por definici´on, decimos que l´ımn→∞ x(n) = n N . 2
C´alculo Num´erico
Contenido
−
x
si ε > 0
Prim
∀
Ult
= = = =
6 25 11 15
(2.28)
−
∃ N (ε) tal que ||
(n)
x
Ant
Sig
Volver
− || < ε para todo x
Cerrar
Abandonar
cuya soluci´on u´ nica es x = (1, 2, 1, 1)t , a partir de la aproximaci´on inicial x(0) = (0, 0, 0, 0)t y con una tolerancia TO L = 5 10 4 . Las ecuaciones del proceso iterativo son
−
x x x
(k +1) 1
=
(k +1) 2
=
(k +1) 3
=
(k +1)
=
x4
×
1 (k ) x2 10
−
1 (k ) x3 5 1 (k ) + x 11 3
1 (k ) x 11 1 1 (k ) 1 (k ) x + x 5 1 10 2 3 (k ) x + 8 2
−
−
3 5 25 + 11 11 10 15 + 8 +
−
3 (k) x 11 4 1 (k) + x 10 4
−
1 (k ) x 8 3
−
y producen el resultado: k (k ) x1 (k ) x2 (k ) x3 (k ) x4
0
1
2
3
4
5
6
7
8
9
10
0
0.6000
1.0473
0.9326
1.0152
0.9890
1.0032
0.9981
1.0006
0.9997
1.0001
0
2.2727
1.7159
2.0530
1.9537
2.0114
1.9922
2.0023
1.9987
2.0004
1.9998
0
-1.1000
-0.8052
-1.0493
-0.9681
-1.0103
-0.9945
-1.0020
-0.9990
-1.0004
-0.9998
0
1.8750
0.8852
1.1309
0.9739
1.0214
0.9944
1.0036
0.9989
1.0006
0.9998
M´etodos Matem´aticos
Contenido
Prim
Ultim
Ant
Sig
Volver
Cerrar
Abandonar
donde hemos parado la iteracio´ n en k = 10 al ser (10)
(9)
||x − x || = 0.327 × 10 ||x ||
4
−
(10)
< TO L
M´etodo iterativo de Gauss-Seidel El m´etodo de Gauss-Seidel supone una aparente mejora del m´etodo de Jacobi ya que no espera al final de cada iteraci o´ n para actualizar el valor de las variables ya calculadas. En este m´etodo, una vez introducidos los valores (0) (0) (0) (1) (x1 , x2 , . . . , xn ) en la primera de las ecuaciones (2.26), y obtenido el valor de x1 , (0) (0) introduciremos dicho dicho valor junto con x2 , . . . , xn en la segunda ecuaci´on para ob(1) (1) (1) (0) (0) tener x2 . A su vez, en la tercera ecuacio´ n introducimos x1 y x2 junto con x3 , . . . , xn (1) para obtener x3 y as´ı sucesivamente. El proceso puede ser reexpresado ahora a partir de C´alculo Num´erico
Contenido
Prim
Ult
Ant
Sig
Volver
Cerrar
Abandonar
(2.26) como
(k +1) 11 x1 (k +1) 21 x1
a a
(k+1)
an1 x1
= b1 = b2
(k+1)
+ a22 x2 +
·· · + a
.. .
(k +1) n,n−1 n−1
x
(k+1)
+ ann xn
−
≡
(k ) 1n xn (k ) 2n xn
−···−a −···−a
= bn (2.29)
o, equivalentemente: (T I + D)x(k+1) = b D) 1 T S y cG S (T I + D) 1 b resulta −
(k ) 12 x2 (k ) 23 x3
−a −a
−
− T x S
x(k+1) = T G
S x
(k )
(k )
−
, de donde haciendo T G
S
−
+ cG
≡ −(T + I
(2.30)
S
−
Ejemplo. Ahora resolveremos el sistema (2.28) por el m´etodo de Gauss-Seidel. En este caso, las ecuaciones del proceso iterativo son M´etodos Matem´aticos
Contenido
Prim
Ultim
Ant
Sig
Volver
Cerrar
Abandonar
x x x
(k +1) 1 (k +1) 2
=
(k +1) 3
=
(k +1)
=
x4
1 (k ) x 10 2
=
1 (k ) x 5 3 1 (k ) x 11 3
−
3 5 25 + 11 11 10 15 + 8 +
1 (k+1) 3 (k ) x + x 11 1 11 4 1 (k+1) 1 (k+1) 1 (k ) x1 + x2 + x4 5 10 10 3 (k+1) 1 (k+1) + x2 x 8 8 3
−
−
−
−
y producen el resultado:
k (k ) x1 (k ) x2 (k ) x3 (k ) x4 C´alculo Num´erico
0
1
2
3
4
5
0
0.6000
1.0302
1.0066
1.0009
1.0001
0
2.3273
2.0369
2.0036
2.0003
2.000
0
-0.9873
-1.0145
-1.0025
-1.0003
-1.000
0
0.8789
0.9843
0.9984
0.9998
1.000
Contenido
Prim
Ult
Ant
Sig
Volver
Cerrar
Abandonar
donde hemos parado la iteracio´ n en k = 5 al ser (5)
(4)
||x − x || = 2.09 × 10 ||x ||
4
−
(5)
< TOL
´ Convergencia de los metodos iterativos La comparaci´on entre los ejemplos de Jacobi (subsecci´on:2.1.2) y Gauss-Seidel (subsecci´on: 2.1.2), podr´ıa hacernos pensar que el m´etodo de iterativo de Gauss-Seidel es mejor que el de Jacobi para resolver sistemas de ecuaciones, ya que aqu´el convergi´o mucho m´as r´apidamente que e´ ste. Sin embargo, esto no es siempre as´ı. Hay casos en que el m´etodo de Jacobi converge m´as r´apidamente, o tambi´en en que uno de los dos m e´ todos converge y el otro no. Veamos primeramente una condici´on suficiente para que los m´etodos iterativos converjan.
Teorema 2.1.2. Si A es estrictamente diagonal dominante, entonces con cualquier elecci´ on de x(0) tanto el m´ etodo de Jacobi como el de Gauss-Seidel generan sucesiones ´ unica ´ del sistema Ax = b. [x(k) ]k=0 que convergen a la soluci on ∞
Obs´ervese que la matriz asociada del sistema de los ejemplos 2.1.2 y 2.1.2 es estrictamente diagonal dominante, aunque esta condici´on es suficiente. Puede haber sistemas M´etodos Matem´aticos
Contenido
Prim
Ultim
Ant
Sig
Volver
Cerrar
Abandonar
con matrices que no sean estrictamente diagonal dominantes y que den lugar a esquemas iterativos convergentes. Para intentar distinguir entre las propiedades de convergencia de ambos m´etodos, introduciremos previamente la siguiente definici´on. Definici´ on. Radio espectral . El radio espectral ρ(A) de una matriz A de autovalores λ1 , . . . , λ p se define como ρ(A) = m´ax λi
| |
1i p
Pues bien, aunque no entraremos en detalle en este punto, puede comprobarse que una estimaci´on para el error absoluto en el paso k del proceso iterativo x(k+1) = T x(k) + c viene dada por (k )
k
(0)
||x − x|| ρ(T ) ||x − x||
(2.31)
es decir, la propagaci´on del error inicial depende del radio espectral. Si queremos que el proceso converja, e´ ste tendr´a que ser ρ(T ) < 1. El m´etodo —Jacobi o Gauss-Seidel— que m´as r´apidamente converger´a en cada caso particular ser´a aquel cuya matriz asociada tenga un menor radio espectral. C´alculo Num´erico
Contenido
Prim
Ult
Ant
Sig
Volver
Cerrar
Abandonar
Refinamiento iterativo Independientemente del m´etodo —exacto o iterativo— que uti˜ licemos para resolver num´ericamente un sistema de ecuaciones Ax = b, la soluci´on x que obtendremos ser´a aproximada, de tal forma que no resolver a´ exactamente el sistema x b = 0. Llamaremos vector resino que al introducirla en el mismo, obtendremos A˜ x (este vector tendr´ıa que calcularse en precisi´on doble sidual a la diferencia r = b A˜ ˜=x x ˜ entre la ya que de otro modo obtendr´ıamos r = 0). Es claro que la diferencia y soluci´on exacta y la aproximada satisface
−
−
−
y=r A˜
(2.32)
La soluci´on del sistema (2.32) —que tendr´ıa que obtenerse en precisio´ n doble— nos ˜ para sustituirla por x ˜+y ˜ . Este es el llamado m´etodo de permite mejorar la soluci´on x refinamiento iterativo.
M´etodos Matem´aticos
Contenido
Prim
Ultim
Ant
Sig
Volver
Cerrar
Abandonar
2.2.
Autovalores
El c´alculo de autovalores y autovectores de matrices est´a presente en m´ultiples aplicaciones a la ingernier´ıa, desde el estudio de las vibraciones de una viga al de la estabilidad de un aeroplano. Son muchos los m e´ todos num´ericos —exactos y aproximados— desarrollados con el fin de diagonalizar matrices. Nosotros nos limitaremos al estudio de uno de ellos, dirigido a obtener un autovalor (y su autovector asociado) de particular importancia en las aplicaciones, el que tiene mayor magnitud, al que llamaremos valor propio dominante.
2.2.1.
M´etodo de la potencia
Consideraremos una matriz A diagonalizable con n autovalores distintos λ1 , λ2 , . . . , λn tales que λ1 > λ2 > . . . > λn . Cada valor propio λi tiene asociado un vector propio vi tal que Avi = λi vi para i = 1, . . . , n. El sistema [v1 , v2 , . . . , vn ] formar´a una base de n n R de modo que cualquier vector x R podr´ a escribirse como una combinacio´ n lineal de los autovectores
| | | |
| |
∈
x = α1 v1 + α2 v2 + C´alculo Num´erico
Contenido
Prim
··· + α v
Ult
n
Ant
(2.33)
n
Sig
Volver
Cerrar
Abandonar
Multiplicando por A a ambos lados de la igualdad (2.33) obtenemos
Ax = α1 λ1 v1 + α2 λ2 v2 +
· ·· + α λ v n n
(2.34)
n
y si repetimos la operaci´on k veces resulta:
Ak x = α1 λk1 v1 + α2 λk2 v2 +
k n n
··· + α λ v
(2.35)
n
Podemos reescribir la ecuaci´on (2.35) en la forma
λk2 A x = λ1 (α1 v1 + α2 k v2 + λ1 k
k
···
λkn + αn k vn ) λ1
(2.36)
donde es claro que, al ser λ1 el autovalor dominante, siempre y cuando α1 = 0 el resto de los t´erminos ser´an despreciables frente al primero en la ecuaci´on (2.36) para valores altos de k . Lo anterior nos da la idea para definir la siguiente iteraci´on. Escogemos un vector x(0) aleatorio con x(0) = 1 y formamos la sucesio´ n:
|| ||
y M´etodos Matem´aticos
(k +1)
(k )
= Ax ,
Contenido
Prim
(k +1)
x
Ultim
=
y(k+1)
||y
Ant
(k+1)
Sig
(2.37)
|| Volver
Cerrar
Abandonar
Por lo expuesto arriba es claro que [x(k) ]k=0 forma una sucesio´ n de vectores unitarios tal que l´ım x(k) = v1 ∞
k→∞
Por tanto, consideraremos x(k) como la aproximaci´on al autovector dominante v1 en el paso k. En dicho paso, tomaremos como aproximaci o´ n al autovalor dominante
λ(k) = x(k) Ax(k) ya que x(k) Ax(k)
·
(k )
(k )
≈λ x ·x 1
= λ1 x(k)
|| ||
·
2
= λ1 .
Ejemplo. Utilizar el m´etodo de la potencia para calcular el autovalor dominante y el correspondiente autovector para la matriz
−4 A = −4 −
(0)
utilizando la aproximaci´on inicial x 10 3 . El resultado exacto es: λ1 = 6 y
0 0
14 13 1 0 2
√ √ √ = (1/ 3, 1/ 3, 1/ 3)
t
y una tolerancia TO L =
−
v1 C´alculo Num´erico
1 √ = (28, 20, −7) ≈ (0.797400, 0.569572, −0.199350) 1233 Contenido
Prim
Ult
Ant
Sig
Volver
Cerrar
Abandonar
El resultado del proceso iterativo (2.37) para este caso es:
k
λ(k)
1 2 3 4 5 6 7 8 9
6.933334 6.475747 6.230470 6.112702 6.055632 6.027620 6.013756 6.006863 6.003428
(k )
(k )
x1
0.778499 0.796858 0.798244 0.797990 0.797721 0.797564 0.797482 0.797441 0.797420
(k )
x2
x3
0.622799 0.597644 0.583332 0.576326 0.572909 0.571228 0.570396 0.569983 0.569777
0.077850 -0.088540 -0.150097 -0.176237 -0.188194 -0.193884 -0.196649 -0.198009 -0.198683
donde hemos parado el proceso en k = 9 ya que (9)
(8)
||x − x || = 7.04 × 10
4
−
< TO L
´ Metodo de la potencia inversa Podemos modificar el m´etodo de la potencia expuesto anteriormente para encontrar el autovalor de A m´as cercano a un cierto nu´ mero q . Este M´etodos Matem´aticos
Contenido
Prim
Ultim
Ant
Sig
Volver
Cerrar
Abandonar
m´etodo se basa en que si q = λi para i = 1, . . . , n, los autovalores de la matriz (A son
1
λ1
1
(A
− qI )
k
x=
1 (λ1
−
1
−
1
− q , λ − q , ··· , λ − q 2
n
Si suponemos que λ1 es el autovalor m´as cercano a q y hacemos actuar (A la expresi´on (2.33) obtenemos −
− qI )
(λ1 (α1 v1 + α2 q )k (λ2
− q ) − q )
k k
v2 +
·· ·
(λ1 + αn (λn
− qI )
1
−
sobre
k
− q ) − q )
k
vn )
donde, nuevamente, los u´ ltimos t´erminos son despreciables frente al primero a condici o´ n de que α1 = 0. En el proceso iterativo ahora habr´a que plantear la ecuaci´on
y(k+1) = (A
− qI )
1 (k )
−
x
(2.38)
que nos obligar´a en cada paso a resolver el sistema de ecuaciones x(k) = (A
− qI )y
(k +1)
(2.39)
por alguno de los m´etodos discutidos en la secci´on 2.1. C´alculo Num´erico
Contenido
Prim
Ult
Ant
Sig
Volver
Cerrar
Abandonar
