UNIVERSIDAD NACIONAL DE INGENIERIA

FACULTAD DE INGENIERIA ELECTRICA Y ELECTRONICA

Laboratorio N°5
Métodos Numéricos

Sistemas Lineales

Profesor: Dr. Israel Dı́az Acha

Alumno Codigo Especialidad

Huaraca Lapa Sebastián Alonso 20200075I Ing. Eléctrica
Malaca Huarac Luis Figo 20204178G Ing. Eléctrica
Galván Pérez Walter Antony 20200343C Ing. Eléctrica
López Cervantes Rogger Francis 20203014K Ing. Eléctrica
Mas Garcı́a Jhosselin Saray 20202244B Ing. Eléctronica

Fecha de ejecución: Sábado 27 de Noviembre

Fecha de entrega: Miércoles 22 de Diciembre

1
Laboratorio N°4
Sistemas Lineales

1. Métodos Directos

a) Sea T (n) = 15 + 25 + 35 + · · · + n5 = a1 · n + a2 · n2 + · · · + a6 · n6 la formula de la suma

de las primeras n quintas potencias de numeros naturales.

Encuentre un sisrema lineal Ax = B, para resolver x = (a1 , a2 , · · · , a6 ) a partir de los

valores de T (1) = 1, T (2) = 33, · · · , T (6).

Por el metodo de eliminacion Gaussiana, implemente un programa en MatLab que permita

triangularizar superiormente la matriz A y encuentre los valores a1 , a2 , · · · , aτ .

Solucion:

T (1) = 1 = a1 + a2 + a3 + a4 + a5 + a6

T (2) = 33 = 2a1 + 4a2 + 8a3 + 16a4 + 32a5 + 64a6

T (3) = 276 = 3a1 + 9a2 + 27a3 + 81a4 + 243a5 + 729a6

T (4) = 1300 = 4a1 + 16a2 + 64a3 + 256a4 + 1024a5 + 4096a6

T (5) = 4425 = 5a1 + 25a2 + 125a3 + 625a4 + 3125a5 + 15625a6

T (6) = 12201 = 6a1 + 36a2 + 216a3 + 1296a4 + 7776a5 + 46656a6

entonces el sistema de la forma Ax = B, es:

    
1 1 1 1 1 1 a
  1  1 
    
2 4 8 16 32 64  a   33 
 2
   
    
3 9 27 81 243 729  a3   276 
    
   =  
    
4 16 64 256 1024 4096  a4   1300 
    
    
5 25 125 625 3125 15625 a   4425 
   5  
    
6 36 216 1296 7776 46656 a6 12201

UNI - Facultad de Ingenieria Electrica y Electronica 2

Laboratorio N°4
Sistemas Lineales

Metodo de eliminacion Gaussiana

Primero escribiremos la matriz ampliada:

 
1 1 1 1 1 1 1 
 
2 4 8 16 32 64 33 
 
 
3 9 27 81 243 729 276 
 
 
 
4 16 64 256 1024 4096 1300 
 
 
5 25 125 625 3125 15625 4425 
 
 
6 36 216 1296 7776 46656 12201

Formaremos la matriz triangular superior mediante operaciones:

 
f2 : f1 × − aa21 + f2 = f1 (−2) + f2
 11 
f3 : f1 × − aa11
31
+ f3 = f1 (−3) + f3
 
f4 : f1 × − aa11
41
+ f4 = f1 (−4) + f4
 
f5 : f1 × − aa11
51
+ f5 = f1 (−5) + f5
 
f6 : f1 × − aa11
61
+ f6 = f1 (−6) + f6

Quedaria:  
1 1 1 1 1 1 1 
 
0 2 6 14 30 62 31 
 
 
0 6 24 78 240 726 273 
 
 
 
0 12 60 252 1020 4092 1296 
 
 
0 20 120 620 3120 15620 4420 
 
 
0 30 210 1290 7770 46650 12195
 
f3 : f2 × − aa22
32
+ f3 = f2 (−3) + f3
 
f4 : f2 × − aa22
42
+ f4 = f2 (−6) + f4
 
f5 : f2 × − aa22
52
+ f5 = f2 (−10) + f5
 
f6 : f2 × − aa62
22
+ f6 = f2 (−15) + f6

UNI - Facultad de Ingenieria Electrica y Electronica 3

Laboratorio N°4
Sistemas Lineales

Quedaria:  
1 1 1 1 1 1 1
 
0 2 6 14 30 62 31 
 
 
0 0 6 36 150 540 180 
 
 
 
0 0 24 168 840 3720 1110 
 
 
0 0 60 480 2820 15000 4110 
 
 
0 0 120 1080 7320 45720 11730
 
f4 : f3 × − aa43 + f4 = f3 (−4) + f4
 33 
f5 : f3 × − aa33
53
+ f5 = f3 (−10) + f5
 
f6 : f3 × − aa63
33
+ f6 = f3 (−20) + f6

Quedaria:  
1 1 1 1 1 11 
 
0 2 6 14 30 62 31 
 
 
0 0 6 36 150 540 180 
 
 
 
0 0 0 24 240 1560 390 
 
 
0 0 0 120 1320 9600 2310
 
 
0 0 0 360 4320 34920 8130
 
f5 : f4 × − aa54 + f5 = f4 (−5) + f5
 44 
f6 : f4 × − aa44
64
+ f6 = f4 (−15) + f6

Quedaria:  
1 1 1 1 1 1 1 
 
0 2 6 14 30 62 31 
 
 
0 0 6 36 150 540 180 
 
 
 
0 0 0 24 240 1560 390 
 
 
0 0 0 0 120 1800 360 
 
 
0 0 0 0 720 11520 2280

UNI - Facultad de Ingenieria Electrica y Electronica 4

Laboratorio N°4
Sistemas Lineales
 
f6 : f5 × − aa65
55
+ f6 = f5 (−6) + f6

 
1 1 1 1 1 1 1 
 
0 2 6 14 30 62 31 
 
 
0 0 6 36 150 540 180
 
 
 
0 0 0 24 240 1560 390
 
 
0 0 0 0 120 1800 360
 
 
0 0 0 0 0 720 120

Ahora llamanos los valores de a1 , a2 , a3 , a4 , a5 y a6 :

1
a6 = 6
360−1800( 16 )
a5 = 120 = 12
390−1560( 16 )−240( 12 ) 5
a4 = 24 = 12
180−540( 16 )−150( 12 )−36( 125
)
a3 = 6 =0
31−61( 16 )−30( 12 )−14( 12
5
)−6(0) 1
a2 = 2 = − 12
1 1 5 1


a1 = 1 − 6 − 2 − 12 − 0 − − 12 =0

 
1 5 1 1
∴ x = 0, − , 0, , ,
12 12 2 6

UNI - Facultad de Ingenieria Electrica y Electronica 5

Laboratorio N°4
Sistemas Lineales

Algoritmo en Matlab

Listing 1: Algoritmo que permite triangularizar la matriz y dar los valores de la solucion

1 A = zeros (6 ,7) ;

2 a =0;

3 for i =1:6

4 for j =1:7

5 if ( j ==7)

6 A (i , j ) = a +( i .^5) ;

7 a = A (i , j ) ;

8 else

9 A (i , j ) = i .^ j ;

10 end

11 end

12 end

13 fprintf ( ' La matriz aumentada : ') ;

14 A

15

16 for i =2:6

17 b = A (i ,1) ;

18 for j =1:7

19 A (i , j ) = A (i , j ) + A (1 , j ) *( - b / A (1 ,1) ) ;

20 end

21 end

22

23 for i =3:6

24 b = A (i ,2) ;

25 for j =1:7

26 A (i , j ) = A (i , j ) + A (2 , j ) *( - b / A (2 ,2) ) ;

27 end

28 end

UNI - Facultad de Ingenieria Electrica y Electronica 6

Laboratorio N°4
Sistemas Lineales

29

30 for i =4:6

31 b = A (i ,3) ;

32 for j =1:7

33 A (i , j ) = A (i , j ) + A (3 , j ) *( - b / A (3 ,3) ) ;

34 end

35 end

36

37 for i =5:6

38 b = A (i ,4) ;

39 for j =1:7

40 A (i , j ) = A (i , j ) + A (4 , j ) *( - b / A (4 ,4) ) ;

41 end

42 end

43

44 for i =6:6

45 b = A (i ,5) ;

46 for j =1:7

47 A (i , j ) = A (i , j ) + A (5 , j ) *( - b / A (5 ,5) ) ;

48 end

49 end

50

51 fprintf ( ' Matriz triangular superior : ') ;

52 A

53 a6 = A (6 ,7) / A (6 ,6) ;

54 a5 =( A (5 ,7) -A (5 ,6) * a6 ) / A (5 ,5) ;

55 a4 =( A (4 ,7) -A (4 ,6) * a6 - A (4 ,5) * a5 ) / A (4 ,4) ;

56 a3 =( A (3 ,7) -A (3 ,6) * a6 - A (3 ,5) * a5 - A (3 ,4) * a4 ) / A (3 ,3) ;

57 a2 =( A (2 ,7) -A (2 ,6) * a6 - A (2 ,5) * a5 - A (2 ,4) * a4 - A (2 ,3) * a3 ) / A (2 ,2) ;

58 a1 =( A (1 ,7) -A (1 ,6) * a6 - A (1 ,5) * a5 - A (1 ,4) * a4 - A (1 ,3) * a3 - A (1 ,2) * a2 )

UNI - Facultad de Ingenieria Electrica y Electronica 7

Laboratorio N°4
Sistemas Lineales

/ A (1 ,1) ;

59

60 fprintf ( ' Los valores de a1 , a2 , a3 , a4 , a5 y a6 son :\ n ')

61 fprintf ( ' a1 = %.04 f \ n ' , a1 ) ;

62 fprintf ( ' a2 = %.04 f \ n ' , a2 ) ;

63 fprintf ( ' a3 = %.04 f \ n ' , a3 ) ;

64 fprintf ( ' a4 = %.04 f \ n ' , a4 ) ;

65 fprintf ( ' a5 = %.04 f \ n ' , a5 ) ;

66 fprintf ( ' a6 = %.04 f \ n ' , a6 ) ;

UNI - Facultad de Ingenieria Electrica y Electronica 8

Laboratorio N°4
Sistemas Lineales

Resultado de la ejecucion

Figura 1: Resultados de la ejecución

Podemos ver que nos sale los mismos valores obtenidos de forma teórica. Entonces, el

algoritmo funciona de manera correcta.

(i+j)!
b) Sea A = (aij ) es una matriz cuadrada de 6x6, definida por aij = i!j! .

Implemente un algoritmo en Matlab que permita encontrar la descomposicion LU por el

metodo de Doo-little.

Solucion:

UNI - Facultad de Ingenieria Electrica y Electronica 9

Laboratorio N°4
Sistemas Lineales

Hallando la matriz A

 
a11 a12 a13 a14 a15 a16 
 
a a22 a23 a24 a25 a26 
 21 
 
a31 a32 a33 a34 a35 a36 
 
A=



a41 a42 a43 a44 a45 a46 
 
 
a a52 a53 a54 a55 a56 
 51 
 
a61 a62 a63 a64 a65 a66

 
2! 3! 4! 5! 6! 7!
 1!1! 1!2! 1!3! 1!4! 1!5! 1!6! 
 
 3! 4! 5! 6! 7! 8! 
 2!1! 2!2! 2!3! 2!4! 2!5! 2!6! 
 
 4! 5! 6! 7! 8! 9! 
 3!1! 3!2! 3!3! 3!4! 3!5! 3!6! 
A=
 5!

6! 7! 8! 9! 10! 


 4!1! 4!2! 4!3! 4!4! 4!5! 4!6! 
 
 6! 7! 8! 9! 10! 11! 
 5!1! 5!2! 5!3! 5!4! 5!5! 5!6! 
 
7! 8! 9! 10! 11! 12!
6!1! 6!2! 6!3! 6!4! 6!5! 6!6!

 
2 3 4 5 7 6
 
3 6 10 15 21 28 
 
 
4 10 20 35 56 084
 
A= =L·U
 
5 15 35 70 126 210
 
 
6 21 56 126 252 462
 
 
7 28 84 210 462 924

Ahora que ya tenemos la matriz A, haremos el algoritmo de matlab que nos permita

encontrar los valores de L y U.

UNI - Facultad de Ingenieria Electrica y Electronica 10

Laboratorio N°4
Sistemas Lineales

Algoritmo en Matlab

Listing 2: Algoritmo para la descomposición LU por el método de Doolitle

1 A = zeros (6) ; L = zeros (6) ; U = zeros (6) ;

2 for i =1:6

3 for j =1:6

4 A (i , j ) = factorial ( i + j ) /( factorial ( i ) * factorial ( j ) ) ;

5 end

6 end

8 for j =1:6

9 L (j , j ) =1;

10 end

11 for j =1:6

12 U (1 , j ) = A (1 , j ) ;

13 end

14 for i =2:6

15 for j =1:6

16 for k =1: i -1

17 v =0;

18 if k ==1

19 v =0;

20 else

21 for p =1: k -1

22 v = v + L (i , p ) * U (p , k ) ;

23 end

24 end

25 L (i , k ) =( A (i , k ) -v ) / U (k , k ) ;

26 end

27 for k = i :6

28 s2 =0;

UNI - Facultad de Ingenieria Electrica y Electronica 11

Laboratorio N°4
Sistemas Lineales

29 for p =1: i -1

30 s2 = s2 + L (i , p ) * U (p , k ) ;

31 end

32 U (i , k ) = A (i , k ) - s2 ;

33 end

34 end

35 end ; L , U

Resultado de la ejecucion

Figura 2: resultados de MatLab

UNI - Facultad de Ingenieria Electrica y Electronica 12

Laboratorio N°4
Sistemas Lineales

c) Sea P (x1 , x2 , x3 , ...., x6 ) = (x1 + x2 + x3 + .... + x6 )2 + (x1 + x2 + x3 + .... + x5 )2 + ... +

(x1 + x2 )2 + (x1 )2

Se define la matriz simetrica A = (aij ) a partir del polinomio cuadratico P de modo que

se cumple:
X
P (x1 , x2 , x3 , ...., x6 ) = ai,j xi xj
i,j

Encuentre los coeficientes de la Matriz A y pruebe que A es una matriz definida positiva.

(Puede usar el Matlab si desea calcular los autovalores de la Matriz A).

Implemente un programa en Matlab que permita encontrar la factorizacion de Cholesky

de la Matriz A.

Solución:

Primero hallamos A, de la forma de P se forma la matriz B y B T

   
1 1 1 1 1 1 1 1 1 1 1 1
   
1 1 1 1 1 0 1 1 1 1 1 0
   
   
1 1 1 1 0 0 1 1 1 1 0 0
   
B= T
yB = 
   
1 1 1 0 0 0 1 1 1 0 0 0
   
   
1 1 0 0 0 0 1 1 0 0 0 0
  
   
1 0 0 0 0 0 1 0 0 0 0 0

Entonces A = B.B T

 
6 5 4 3 2 1
 
5 5 4 3 2 1
 
 
4 4 4 3 2 1
 
A=



3 3 3 3 2 1
 
 
2 2 2 2 2 1
 
 
1 1 1 1 1 1

Calculamos los autovalores de A para determinar si es una matriz positiva.

Algoritmo en Matlab

UNI - Facultad de Ingenieria Electrica y Electronica 13

Laboratorio N°4
Sistemas Lineales

Listing 3: Autovalores de la matriz A

1 %Hallamos los autovalores de la matriz A

2 A =[6 5 4 3 2 1; 5 5 4 3 2 1; 4 4 4 3 2 1; 3 3 3 3 2 1; 2 2 2

2 2 1; 1 1 1 1 1 1];

3 lambda = eig ( A ) ;

4 disp ( lambda ) ;

Resultado de la ejecución

Figura 3: Autovalores obtenidos del codigo en matlab

Como se muestra en los resultados de ejecucion, todos los autovalores son positivos, por

lo tanto podemos afirmar que A es una matriz definida positiva.

Ahora hallamos la matriz L a partir del metodo de Cholesky, para ello utilizaremos Matlab

y ademas comprobamos que se cumple A = L × LT .

Algoritmo en Matlab

Listing 4: Metodo de Cholesky

1 %Calculo de L por el metodo de Cholesky

2 clear

3 clc

4 L = zeros (6) ;

UNI - Facultad de Ingenieria Electrica y Electronica 14

Laboratorio N°4
Sistemas Lineales

5 A =[6 5 4 3 2 1; 5 5 4 3 2 1; 4 4 4 3 2 1; 3 3 3 3 2 1; 2 2 2

2 2 1; 1 1 1 1 1 1];

6 for i =1:6

7 sum =0;

8 for j =1: i -1

9 sum = sum + ( L (i , j ) ) ^2;

10 end

11 L (i , i ) = sqrt ( A (i , i ) - sum ) ;

12 for k = i +1:6

13 sum =0;

14 for j =1: i -1

15 sum = sum + L (k , j ) * L (i , j ) ;

16 end

17 L (k , i ) =( A (k , i ) - sum ) / L (i , i ) ;

18 end

19 end

20 disp (" Matriz L ")

21 disp ( L ) ;

22 disp (" Matriz L ^ T ")

23 disp ( transpose ( L ) ) ;

24 B = L * transpose ( L ) ;

25 disp (" A = L * L ^ T ") ;

26 disp ( B ) ;

UNI - Facultad de Ingenieria Electrica y Electronica 15

Laboratorio N°4
Sistemas Lineales

Resultado de la ejecución

Figura 4: Matrices L y LT

Como se muestra en la figura 4, la matriz L se obtiene a partir del metodo de Cholesky y

ademas se demuestra que A = L × LT .

2. Métodos iterativos

12x − 2y − z = 3

−x + 12y − 2z = 8

−2x − y + 12z = −12

a) Encuentre la matriz Mj de Jacobi para el proceso iterativo Xn+1 = Mj ∗ xn + C. Pruebe

que el método de Jacobi es convergente por el método del radio espectral de los autovalores

(use calculadora). Calcule el número de iteraciones necesarias para conseguir un error de

T < 10−4 con la norma suma ∥.∥1 .

Solución:

UNI - Facultad de Ingenieria Electrica y Electronica 16

Laboratorio N°4
Sistemas Lineales

Cálculo de la matriz de Jacobi

Para el enunciado dado transformaremos el sistema de ecuaciones a una ecuación matricial

de la siguiente manera.
    
 12 −2 −1 x
   3 
    
−1 12 −2 y  =  8 
    
    
−2 −1 12 z −12

Teniendo en cuenta que se quiere calcular de forma general esta expresión por el método

de Jacobi definiremos la matriz pedida a continuación.

AX = b

(D + L + U )X = b

DX = −(L + U )X + b

X = −D−1 (L + U )X + D−1 b

X = −D−1 RX + D−1 b

A la expresión −D−1 R la llamaremos la matriz de Jacobi y en el ejercicio dado tendrı́a la

siguiente forma:

   
12 0 0   0 −2 +1
   
D=
 0 12 0  , R = −1 0 −2
  
   
0 0 12 −2 −1 0
 
1 1
0 6 12 

MJ = −D−1 R = 
 
1 1
 12 0 6


 
1 1
6 12 0

UNI - Facultad de Ingenieria Electrica y Electronica 17

Laboratorio N°4
Sistemas Lineales

Convergencia del método de jacobi

Para ello calcularemos el radio espectral (ρ) de la matriz de Jacobi (D−1 R). Teniendo en

cuenta esto la matriz por analizar serı́a la siguiente.

 
1 1
0 6 12 

MJ = −D−1 R = 
 
1 1
 12 0 6


 
1 1
6 12 0

Cuyos autovalores (λ) serán calculados a partir de la siguiente relación:

P (λ) = det(MJ − λI) = 0

Reduciéndose a la expresión:

192λ3 − 8λ − 1 = 0

Dándonos como resultados los siguientes autovalores:

1 1 195i 1 195i
λ1 = , λ2 = − + , λ3 = − −
4 8 2702 8 2702

El radio espectral se define: ρ(MJ ) =max(|λ|,λ es el máximo autovalor de MJ ).

Conociendo ya estos autovalores se puede concluir que el radio espectral de la matriz ana-
1
lizada es ρ = 4 y a su vez se observa que este valor es menor que 1 comprobando ası́ la

fiabilidad del método de Jacobi para la solución de este problema.

Cálculo del número de iteraciones

Conociendo la fórmula del número de iteraciones procederemos a utilizarla.

T (1−λ )
log( ∥x1 −x0p∥1 )
n≥
log λp

Para encontrar estos valores calcularemos nuestra fórmula de iteración para el método de

UNI - Facultad de Ingenieria Electrica y Electronica 18

Laboratorio N°4
Sistemas Lineales

jacobi.       
1 1 1
Xn+1   0 6 12  Xn    4
      
 1 1  Y  +  2 
 n+1  =  12
Y 0 6  n  3 
      
1 1
Zn+1 6 12 0 Zn −1

1 2 23
De esta expresión podemos calcular ∥x1 − x0 ∥1 = 4 + 3 +1= 12 y reemplazaremos todo

en la fórmula.
10−4 (1− 14 )
log( 23 )
12
n≥
log 14

n ≥ 7, 32067

Por lo tanto el número de iteraciones será de 8.

b) Implemente un programa en Matlab que permita encontrar la solución del sistema por

el método iterativo de Jacobi con un error T < 10−4 con la norma suma ∥.∥1 , elija la

condición inicial (x0 , y0 , z0 ) = (0, 0, 0) como punto inicial.

Solución:

Algoritmo en Matlab

Listing 5: Método de Jacobi

1 clear

2 clc

3 format long

4 A = input ( ' ingrese la matriz estrictamente diagonal dominante A

: ') ;

5 b = input ( ' ingrese el vector b como vector fila : ') ;

6 x0 = input ( ' ingrese el vector x0 como vector fila : ') ;

7 x = x0 ';

8 Tol = input ( ' ingrese la tolerancia : ') ;

9 max = input ( ' ingrese el numero maximo de iteraciones : ') ;

10 disp ( ' ')

UNI - Facultad de Ingenieria Electrica y Electronica 19

Laboratorio N°4
Sistemas Lineales

11 disp ( ' ')

12 n = length ( b ) ;

13 D = zeros ( n ) ;

14 T = zeros ( n ) ;

15 c = zeros (n ,1) ;

16 for i =1 : n ; D (i , i ) = A (i , i ) ; end

17 T =A - D ;

18 T = - D \ T ;

19 c = D \b ';

20 P = zeros ( max , n ) ;

21 k =1;

22 while k < max +1

23 y=T*x+c;

24 dx =y - x ;

25 x = y ; P (k ,1: n ) =x ';

26 if norm ( dx ,1) < Tol

27 last = k ;

28 max = k ;

29 end

30 k = k +1;

31 end

32 k = k +1;

33 disp ( ' ')

34 disp ( ' ')

35 disp ( ' la solucion al sistema [ A b ] : ')

36 disp ( ' ')

37 disp ([ A b '])

38 disp ( ' Usando el metodo de Jacobi con el vector inicial x0 = ')

39 disp ( ' ')

40 disp ( x0 )

UNI - Facultad de Ingenieria Electrica y Electronica 20

Laboratorio N°4
Sistemas Lineales

41 disp ([ ' con tolerancia = ' , num2str ( Tol ) , ' es ' ])

42 disp ( ' ')

43 disp (x ')

44 disp ( ' la lista de los x ' 's ')

45 disp ( ' ')

46 disp ( P (1: last , 1: n ) )

47 disp ( ' ')

48 disp ( ' el error dx = ')

49 disp ( dx ')

50 disp ([ ' para obtener la aproximacion con una tolerancia = ' ,

num2str ( Tol ) , ' y ' ])

51 disp ([ int2str ( max ) , ' fueron las iteraciones necesarias ' ])

Resultado de la ejecución

ingrese la matriz estrictamente diagonal dominante A: [12 -2 -1; -1 12 -2; -2 -1 12]

ingrese el vector b como vector fila: [3 8 -12]

ingrese el vector x0 como vector fila: [0 0 0]

ingrese la tolerancia : 0.0001

ingrese el numero maximo de iteraciones : 20

Resultado:

X = 0,263511877491641 Y = 0,536734624164095 Z = −0,911350831082819

Fueron necesarias 7 iteraciones.

c) Encuentre la matriz MGS de Gauss-Seidel para el proceso iterativo Xn+1 = MGS .xn + C.

Pruebe que el metodo de Gauss-Seidel es convergente por el metodo de radio espectral de

los autovalores (use calculadora). Calcule el numero de iteraciones necesarias para conse-

guir un error de T < 104 con la norma maximo ∥.∥∞ :

Solucion:

UNI - Facultad de Ingenieria Electrica y Electronica 21

Laboratorio N°4
Sistemas Lineales

Hallando la matriz MGS

donde: A = D + L + U

Sabemos: MGS = −(D + L)−1 U

La matriz A es:  
 12 −2 −1 
 
−1 12 −2
 
 
−2 −1 12
     
12 0 0  0 0 0 0 −2 −1
     
→D=
 0 12 0  , L = −1 0 0 y U = 0 0 −2
    
     
0 0 12 −2 −1 0 0 0 0
 
 12 0 0
 
D+L=
−1 12 0 

 
−2 −1 12
   
144 12 25  144 0 0 
T
   
→ adj(D + L) = 
 0 144 12  → adj (D + L) =  12 144 0 
  
   
0 0 144 25 12 144

→ det(D + L) = 12 × 12 × 12 = 123
 
 144 0 0 
−1 1  
→ (D + L) = 3  12 144 0 

12  

25 12 144

Reemplazamos en MGS :
  
144 0 0  0 −2 −1
1   
MGS =− 3 12 144 0  0 0 −2
12 
 

 

25 12 144 0 0 0

UNI - Facultad de Ingenieria Electrica y Electronica 22

Laboratorio N°4
Sistemas Lineales
 
0 −288 −144
1  
MGS =− 3 0 −24 −300
12 



0 −50 −49
 
1 1
0 6 12 
 
MGS = 1 25
0

72 144 
 
50 49
0 1728 1728

Convergencia del metodo de Gauss-Seidel

Para que el metodo converga,|λ| debe ser menor a 1, donde:

P (λ) = det(MGS − λI) = 0

 
1 1
−λ 6 12 
 
→ (MGS − λI) = 
 0
1
12 −λ 25
144


 
50 49
0 1728 1728 −λ

     
1 49 25 50
P (λ) = −λ −λ −λ +λ =0
12 1728 144 1728

Usando geogebra para las soluciones:

Figura 5: hallando los valores de λ en geogebra

UNI - Facultad de Ingenieria Electrica y Electronica 23

Laboratorio N°4
Sistemas Lineales

→ λ1 = −0,02 , λ2 = 0 , λ3 = 0,1319

El maximo |λ3 | = 0,1319 y este es menor que 1 |λ3 | < 1

∴ El metodo de Gauss-Seidel es convergente.

Numero de iteraciones

Sabemos:  
T (1−λ∞ )
log ∥x1 −x0 ∥ ∞
n≥
log λ∞

→ x1 = B = (D + L)−1 b
    
144 0 0  3   432 
1    1  
=B= 3 12 144 0   8 =
3
 1188 
12 
 


 12 
  

25 12 144 −12 −1557
 
1
 4 
 
x1 =  11 
 16 
 
−173
192

−173 173
→ ∥x1 − x0 ∥∞ = =
192 192

λ∞ = ∥MGS ∥ ⇒ maxima suma absoluta de filas

3 27 99 3


λ∞ = M AX 12 , 144 , 1728 = 12

Reemplazamos en la formula:

 3

10−4 (1− 12 )
log 173
192
n≥ 3


log 12

n ≥ 6,766

n=7

∴ 7 iteraciones

UNI - Facultad de Ingenieria Electrica y Electronica 24

Laboratorio N°4
Sistemas Lineales

d) Implemente un programa en Matlab que permita encontrar la solución del sistema por el

método Gauss-Seidel con un error T < 10−4 con la norma maximo ||.||∞ elija la condicion

inicial (x0 , y0 , z0 ) = (0, 0, 0) como punto inicial

solución:

Algoritmo en Matlab

Listing 6: Método Gauss-Seidel

1 % METODO ITERATIVO DE GAUSS SEIDEL

3 clc %permite borrar el area de trabajo

4 clear %permite borrar las variables almacenadas

5 format long %permite utilizar la maxima capacidad de la

maquina

7 fprintf ( ' METODO ITERATIVO DE GAUSS

SEIDEL \ n \ n \ n ')

8 %fprintf me permite ingresar comentarios de manera textual

que pueden

9 %orientar al usuario en el uso del programa

10

11 %input es un comando de solicitud de entrada de datos del

usuario .

12 a = input ( ' Ingrese la matriz de coeficientes :\ n ') ;

13 b = input ( '\ nIngrese los t r m i n o s independientes :\ n ') ;

14 x = input ( '\ nIngrese el vector con las ap ro xi ma ci ma ci on es

Iniciales :\ n ') ;

15 iter = input ( '\ nIngrese el n m e r o m x i m o de iteraciones :\ n ')

16 tol = input ( '\ nIngrese la tolerancia :\ n ') ;

17

UNI - Facultad de Ingenieria Electrica y Electronica 25

Laboratorio N°4
Sistemas Lineales

18 k = norm ( a ) * norm ( a ^ -1) ; %Se calcula el condicional de la matriz

de coeficientes

19 disp ( ' condicional = ')

20 disp ( k )

21 % la funcion disp nos permite imprimir una variable en el

espacio de trabajo

22 determinante = det ( a ) ; %se calcula el determinante de la matriz

de coeficiente

23 if determinante ==0

24 disp ( ' El determinante es cero , el problema no tiene s o l u c i n

nica ')

25 end

26

27 n = length ( b ) ; %numero de elementos del vector b

28 d = diag ( diag ( a ) ) ; %obtencion de la matriz diagonal

29 l =d - tril ( a ) ; %obtencion de la matriz diagonal superior L

30 u =d - triu ( a ) ; %obtencion de la matriz diagonal inferior u

31

32 fprintf ( '\ n SOLUCION :\ n ')

33 fprintf ( '\ nLa matriz de transicion de gauss seidel :\ n ')

34 T =(( d - l ) ^ -1) * u ; % matriz de transicion de gauss

35 disp ( T )

36 re = max ( abs ( eig ( T ) ) ) %calculo del radio espectral

37

38 if re >1

39 disp ( ' Radio Espectral mayor que 1 ')

40 disp ( ' el m t o d o no converge ')

41

42 return

43 end

UNI - Facultad de Ingenieria Electrica y Electronica 26

Laboratorio N°4
Sistemas Lineales

44 fprintf ( '\ nEl vector constante es ::\ n ')

45 C =(( d - l ) ^ -1) * b ; % vector constante C , para el metodo

46 disp ( C )

47 i =0;

48

49 err = tol +1;

50 z =[ i , x (1) ,x (2) ,x (3) , err ]; %vector que me permite graficar la

tabla

51

52 while err > tol & i < iter

53 xi = T * x + C ;

54 %disp ( xi )

55 i = i +1;

56 err = norm ( xi - x ) ; %norma 2

57 %err = max ( abs ( xi - x ) ) ; %norma 1

58 %err = norm ( xi - x ) / norm ( xi ) ; %norma relativa

59

60 x = xi ;

61 z (i ,1) = i ;

62 z (i ,2) = x (1) ;

63 z (i ,3) = x (2) ;

64 z (i ,4) = x (3) ;

65 z (i ,5) = err ;

66 end

67 fprintf ( '\ nTABLA :\ n \ n n x1

x2 x3

Error \ n \ n ')

68 disp ( z ) %impresion de la tabla .

UNI - Facultad de Ingenieria Electrica y Electronica 27

Laboratorio N°4
Sistemas Lineales

Resultado de la ejecución

ingrese la matriz estrictamente diagonal dominante A: [12 -2 -1; -1 12 -2; -2 -1 12]

ingrese el vector b como vector fila: [3;8;-12];

ingrese el vector x0 como vector fila: [0;0;0];

ingrese la tolerancia : 0.0001

ingrese el numero maximo de iteraciones : 30

Resultado:

X = 0,263510665780098 Y = 0,536733972325276 Z = −0,911353724676211

Fueron necesarias 6 iteraciones.

UNI - Facultad de Ingenieria Electrica y Electronica 28

Laboratorio N°4
Sistemas Lineales

Conclusiones y observaciones

1. Conclusiones

El método de eliminación gaussiana se complica cada vez que el orden de la matriz au-

menta, porque cada vez se tiene que hacer mas operación elementales. Pero este método

se harı́a fácil con un algoritmo donde se pueda ingresar cualquier matriz cuadrada

El método de Gauss-Seidel converge el doble de rápido, aproximadamente con respecto

a el método de Jacobi;Ademas, si la matriz proporcionada es diagonal estrictamente

dominante, el método siempre va a converger; Sin embargo, aunque no sea ası́, es posible

que converja. Por ello, ambos métodos son útiles cuando se tiene un sistema de ecuaciones

con diagonal estrictamente dominante y cuando no se conoce otro método alternativo

para hacerlo. por ejemplo, La regla de Cramer, Gauss-Jordan, etc.

El método de Gauss es mucho más eficiente que el de Jacobi, esto debido a la considerable

reducción de memoria necesaria para ejecutarlo. Esto se ve reflejado en el número de

iteraciones.

2. Observaciones

El método de Jacobi realiza muchas más operaciones que el método de Gauss-Seidel esto

debido a que con Gauss el reemplazo de la variable anterior es casi inmediato; por el

contrario, al utilizar Jacobi las variables se reemplazan luego de haber calculado todas

las enteriores.

las fórmulas para hallar el número de iteacciones en el método de Jacobi o de gauss-siedel

solo sirve como una cota, ya que, en realidad se necesita menos o igual iteracciones a la

calculada.

En el metodo de Cholesky si bien se hallo la matriz a partir de la multiplicacion de una

matriz con su transpuesta, esta matriz hallada no es la matriz que se debe de utilizar

en el metodo de Cholesky ya que no es una matriz triangular inferior como lo indica la

teoria, esta primera matriz fue hallada con el metodo directo que se indico en clase.

UNI - Facultad de Ingenieria Electrica y Electronica 29

 Laboratorio N°4
Sistemas Lineales

Algoritmo en Matlab

Listing 7: Algoritmo para graficar la respuesta en frecuencia

1 % Valores de frecuencia que se v a r i

2 f = [10 50 100 200 500 1000 2000 5000 10000 20000 50000 100000 200000

2000000 8000000];

4 % Voltajes de entrada y salida agrupados en vectores

5 Vo = [270 340 344 348 348 300 350 356 352 348 344 340 344 150 92];

6 Vi = [50 56 56 58 56 46 58 58 58 58 58 58 58 56 40];

8 % Hallamos la ganancia en t e n s i n valor a valor

9 Av = Vo ./ Vi ;

10 H =20* log10 ( Av ) ;

11

12 % Graficamos 20 log ( Av ) vs log ( f )

13 semilogx (f ,H , ' LineWidth ' ,1.5 , ' Color ' , 'b ') ;

14 grid on ;

15 xlabel ( ' log ( f ) ') ;

16 ylabel ( ' 20* log ( Av ) ') ;

17 title ( ' Respuesta en frecuencia ') ;

UNI - Facultad de Ingenieria Electrica y Electronica 30