Transcript Métodos directos para sistemas lineales

Sistemas de lineales
Programación Numérica
Objetivos
• Resolver un sistema de ecuaciones lineales
acopladas
• Obtener la inversa de una matriz
• Calcular el determinante de una matriz
• Factorizar matrices
Solución de sistemas pequeños
La solución de un sistema lineal de 2 ecuaciones puede
obtenerse fácilmente de forma gráfica
a 
b
x2   11  x1  1
a12
 a12 
a11x1 + a12x2 = b1
a21x1 + a22x2 = b2
Despejando x2 y graficando
 a11 
b1
x2    x1 
a12
 a12 
 a21 
b2
 x1 
x2  
a22
 a22 
solución
a 
b
x2   21  x1  2
a22
 a22 
Ejemplo
Sea el sistema
3x1 + 2x2 = 18
–x1 + 2x2 = 2
Despejando
x2 = 9 – 1.5x1
x2 = 1 + 0.5x1
(4, 3)
Tipos de soluciones
No hay solución
Ejemplo
Infinidad de soluciones
Ejemplo
Sistema mal condicionado
Ejemplo
½x1 + x2 = 1
½x1 + x2 = 1
½x1 + x2 = 1
½x1 + x2 = ½
x1 + 2x2 = 1
0.501x1 + x2 = 0.9
Determinantes
Solución de un sistema de 3x3
La matriz del sistema es:
 a11
A  a21
a31
a12
a22
a32
a13 
a23 
a33 
El determinante de A es:
a11
a12
a13
D  a21
a22
a23
a31
a32
a33
a22
 a11
a32
a23
a21
 a12
a33
a31
a23
a21
 a13
a33
a31
a22
a32
 a11 a22 a33  a32 a23   a12 a21a33  a31a23   a13 a21a32  a31a22 
Regla de Cramer
Cada incógnita de un sistema de ecuaciones lineales algebraicas
puede expresarse como la razón de dos determinantes con
denominados D y con numerador obtenido a partir de D, al
reemplazar la columna de coeficientes de la incógnita en cuestión
por las constantes independientes b1, b2, …, bn. Por ejemplo x1 se
calcula con:
x1 
b1
a12
a13
b2
a22
a23
b3 a32 a33
a11 a12 a13
a21
a22
a23
a31
a32
a33
Eliminación de Gauss simple
Definición
Un sistema lineal es un sistema de la forma:
E1: a11x1 + a12x2 + a13x3 +...+ a1nxn = b1
E2: a21x1 + a22x2 + a23x3 +...+ a2nxn = b2
.
En: an1x1 + an2x2 + an3x3 +...+ annxn = bn
Donde x1, x2, x3.., xn son las incógnitas.
Operaciones básicas
Se utilizarán tres operaciones para simplificar un sistema
lineal:
1. La ecuación Ei puede multiplicarse por una constante l
distinta de cero. (lEi) (Ei)
2. La ecuación Ej puede multiplicarse por una constante l
distinta de cero y sumarse a otra ecuación Ei. (Ei +lEj) 
(Ei)
3. El orden de las ecuaciones Ei y Ei puede intercambiarse.
(Ei)  (Ei)
Ejemplo
positivo
E1: x1 + x2
+ 3x4 = 4
E2: 2x1+ x2 – x3 +
x4 = 1
E3: 3x1- x2 – x3 + 2x4 = -3
E4: -x1+2x2 +3x3 -
x4 = 4
Las operaciones:
(E2 –2E1)E2, (E3 –3E1) E3, (E4 + E1)E4
Transforman el sistema en:
Simplificar
E1: x1 + x2
+ 3x4 = 4
E2:
- x2 – x3 - 5x4 = -7
E3:
- 4x2 – x3 -7x4 =-15
E4:
+ 3x2 +3x3 + 2x4 = 8
Las operaciones:
(E3 –4E2)E3, (E4 +3E2) E4
Transforman el sistema en:
Simplificar
E1: x1 + x2
+ 3x4 = 4
E2:
- x2 – x3 - 5x4 = -7
E3:
3x3 +13x4 = 13
E4:
-13x4 = -13
El sistema se puede resolver hacia atrás.
x4 = 1
x3 = (13-13x4)/3 = 0
x2 = -(x3 + 5x4 -7) = 2
x1 = -x2 - 3x4 + 4 = -1
Comprobación
E1 : x 1 + x 2
+ 3x4 = 4
E2: 2x1+ x2 – x3 +
x4 = 1
x1 = -1
x2 = 2
E3: 3x1- x2 – x3 + 2x4 = -3
x3 = 0
E4: -x1+2x2 +3x3 -
x4 = 1
x4 = 4
E1: (-1)+(2)+3(1)= 4
E2: 2(-1)+(2)–(0)+(1) = 1
E3: 3(-1)-(2)–(0)+2(1)= -3
E4: -(-1)+2(2)+3(0)-(1) = 4
Solución en Matlab
%sistema de ecuaciones
A = [1,1,0,3,4;2,1,-1,1,1;3,-1,-1,2,-3;-1,2,3,-1,4]
x = [0 0 0 0]’;
%operaciones
%(E2 - 2E1) -> E2
A(2,:) = A(2,:) - 2*A(1,:);
%(E3 - 3E1) -> E3
A(3,:) = A(3,:) - 3*A(1,:);
%(E4 + E1) -> E4
A(4,:) = A(4,:) + A(1,:);
%(E3 - 4E2) -> E3
A(3,:) = A(3,:) - 4*A(2,:);
%(E4 + 3E2) -> E4
A(4,:) = A(4,:) + 3*A(2,:);
x(4) = A(4,5)/A(4,4);
x(3) = (A(3,5)-A(3,4)*x(4))/A(3,3);
x(2) = (A(2,5)-A(2,3:4)*x(3:4))/A(2,2);
x(1) = (A(1,5)-A(1,2:4)*x(2:4))/A(1,1);
Actividad
Resolver
x1 + x2
+ x4 = 2
2x1 + x2 - x3 + x4= 1
4 x1 - x2 - 2x3 + 2 x4 = 0
3x1 - x2 - x3 + 2 x4 = -3
Matrices
Una matriz es un arreglo rectangular de elementos de n
renglones y m columnas.
 a11
a
A  aij    21
 

an1
a12
a22

an 2
 a1m 
 a2 m 
 

 anm 
Representación matricial
El sistema de ecuaciones
a11 x1 a12 x2  a1m xn  b1
a21 x1 a22 x2  a2 m xn  b2



an1 x1 an 2 x2  anm xn  bn
Se puede representar por una matriz ampliada
 a11
a
A, b   21
 

an1
a12
a22

an 2
 a1n . b1 
 a2 n . b2 
 . .

 ann . bn 
Ejemplo
(E2 –2E1)  E2, (E3 –3E1)  E3,
(E4 + E1)  E4
1 1 0 3. 4
 2 1 1 1 . 1 


 3  1  1 2 .  3


 1 2 3  1. 4 
1 1 0 3 . 4 
0  1  1 5 .  7 


0  4  1  7 .  15


0 3 3  1 . 8 
(E3 –4E2)  E3, (E3 +3E2)  E4
3 . 4 
1 1 0
0  1  1 5 .  7 


0 0 3 13 .  13


0
0
0

13
.

13


Eliminación Gaussiana con sustitución hacia
atrás
Se forma la matriz aumentada
 a11

a
A  A, b    21
 

an1
a12
a22

an 2
 a1n . a1,n 1 

 a2 n . a2,n 1 
. 
 .

.
a
 ann
n , n 1 
Siempre que aii  0, realizamos la operación
(Ej–(aji/aii) Ei)  (Ej)
Para i = 1, 2, 3, ..., n-1 y j = i + 1, i + 2, ..., n.
a11

0
A


0
a12  a1n . a1,n 1 

a22  a2 n . a2,n 1 
. 

 .

.
a
0  ann
n , n 1 
Sustitución hacia atrás
an,n 1
xn 
an , n
an1,n1  an1,n xn
xn1 
an1,n1
En general
xi 
ai ,n 1  ai ,n xn  ai ,n1 xn 1  ...  ai ,i 1 xi 1

ai ,i
ai ,n 1   j i 1 ai , j x j
n
ai ,i
Algoritmo Gauss
Entrada: número de incógnitas y de ecuaciones n, matriz aumentada (aij)
donde 1<= i <= n y 1<= j <= n+1.
Salida: solución x1, x2 ... xn o mensaje de que no hay solución
1. Para i = 1 hasta n-1
2.
Sea p el entero más pequeño con
i <= p <= n y api <> 0
3.
Si no se encuentra p
4.
Salida(“no existe solución”)
5.
terminar
continuación
6.
Si p<>i realice (Ei)<->(Ep)
7.
Para j = i+1 hasta n
8.
mij = aij/aii
9.
Realice (Ej-mijEi) -> (Ej)
10. Si ann=0 entonces
11.
Salida(“no existe solución”)
12. xn = an,n+1/ann
13. Para i = n – 1 hasta 1
14.
xi = [ai,n+1-Sj=i+1 aijxj]/aii
Código en C
void gauss(double A[][10],int n, double x[]){
//resuelve un sistema de ecuaciones Mx = b
// A es la matriz aumentada (n x (n+1)) A = [M b]
int i,j,k,p; //p es el pivote
double temp,s,m;
for(i = 0;i<n-1; i++){
p = i;
while(p < n-1 && A[p][i] == 0)
p = p + 1;
if(p == n-1){
std::cout<<"No hay solución por pivoteo\n";
return;
}
if(A[n-1][n-1]==0){
if(p != i)
std::cout<<"No solución\n";
for(j=0; j<n+1; j++){
return;
temp = A[i][j];
}
A[i][j] = A[p][j];
x[n-1]=A[n-1][n]/A[n-1][n-1];
A[p][j] = temp;
for(i = n-2;i>=0; i--){
}
s = A[i][n];
for(j = i+1; j<n; j++){
for (j = i+1; j<n;j++)
m = A[j][i]/A[i][i];
s = s - A[i][j]*x[j];
for(k = 0;k<n+1;k++)
x[i] = s/A[i][i];
A[j][k] = A[j][k] - m*A[i][k];
}
}
}
}
Método de Gauss
function x = ecuaciones(A)
% resuelve un sistema de ecuaciones Mx = b
% A es la matriz aumentada (n x (n+1)) A = [M b]
[n dumy] = size(A);%n = número de renglones
for i = 1 : n-1
p = i;
while p < n && A(p,i) == 0; p = p + 1; end
if p == n
fprintf('No hay solución...');
return;
end
if p ~= i
temp = A(i,:);
A(i,:) = A(p,:);
A(p,:) = temp;
end
for j = i+1 : n
m = A(j,i)/A(i,i);
for k = 1 : n+1
A(j,k) = A(j,k) - m*A(i,k);
end
end
end
Método de Gauss
if A(n,n)==0
fprintf('No hay solución...');
return;
end
x(n) = A(n,n+1)/A(n,n);
for i = n-1 : -1: 1
s = A(i,n+1);
for j = i+1 : n
s = s - A(i,j)*x(j);
end
x(i) = s/A(i,i);
end
Inversa
La inversa de una matriz A de n x n es una matriz A -1 tal que A x
A-1 = A-1 x A = In
Una matriz que tiene inversa se le llama matriz no singular.
Una matriz que no tiene inversa es una matriz singular.
Se cumple que:
a. La inversa es única.
b. La inversa de una matriz no singular es no singular y (A-1)-1 = A
c. Si A y B son no singulares, (AB)-1 = B-1 A-1
Obtención de la inversa
Sea Bj la j-ésima columna de la matriz B de n x n.
 b1 j 
b 
2j
Bj   
  
 
 bnj 
Si AB = C, entonces la j-esima columna de C está dada por:
n
 c1 j 
 a11 a12  a1n   b1 j   k 1 a1k bkj 

c 
a
 b   n
a

a
k 1 a2 k bkj 
22
2n   2 j 
 2 j   C  AB   21

j
j

  
 

    


 

   n
 cnj 
an1 an 2  ann   bnj  k 1 ank bkj 
Cont.
Supongamos que A-1 existe y que A-1 = B = (bij). Entonces AB =
Iy
0 

 
0 
 
AB j  1 
0 
 

0 
 
Para encontrar B podemos resolver n ecuaciones lineales donde la
j-ésima columna de la inversa es la solución del sistema lineal
que en el lado derecho tiene la j-ésima columna de I.
Ejemplo
 1 2  1
A   2 1 0 
 1 1 2 
 b11  2b21  b31
AB   2b11  b21
 b11  b21  2b31
 1 2  1 b11 b12
AB   2 1 0  b21 b22
 1 1 2  b31 b32
b12  2b22  b32
2b12  b22
 b12  b22  2b32
b13 
b23 
b33 
b13  2b23  b33 
2b13  b23 
 b13  b23  2b33 
Si B = A1, entonces AB = I se tiene
b11  2b21  b31  1
b12  2b22  b32  0
b13  2b23  b33  0
2b11  b21  0
2b12  b22  1
2b13  b23  0
 b11  b21  2b31  0  b12  b22  2b32  0  b13  b23  2b33  1
2
5
1
b12 
b13  
9
9
9
4
1
2
b21 
b22  
b23 
9
9
9
1
1
1
b31  
b32 
b33 
3
3
3
b11  
 2
 9
 4
1
A 
 9
 1
 3
5
9
1

9
1
3
1
 
9
2 

9 
1 
3 
Inversa en C
void inversa(double B[][20],int n,double x[][20]){
double I[20][20];
double temp,m,s;
int i,j,k,p;
for(i = 0; i<20; i++)
for(j = 0; j<20; j++)
if(i==j) I[i][j] = 1; else I[i][j] = 0;
double A[20][20];
for(i = 0; i<n; i++)
for(j = 0; j<2*n; j++)
if(j<n) A[i][j] = B[i][j]; else A[i][j] = I[i][j-n];
for(i = 0; i<20; i++)
for(j = 0; j<20; j++)
x[i][j] = 0;
for(i = 0; i<n-1;i++){
p = i;
while(p < n && A[p][i] == 0) p = p + 1;
if( p == n){
std::cout << "Matriz singular...\n";
return;
}
if(
p != i)
for( k = 0; k<2*n; k++){
temp = A[i][k];
A[i][k] = A[p][k];
A[p][k] = temp;
}
for(j = i+1;j<n; j++){
m = A[j][i]/A[i][i];
for( k = 0; k<2*n; k++)
A[j][k] = A[j][k] - m*A[i][k];
}
}
if( A[n-1][n-1]==0){
std::cout << "Matriz singular...\n";
return;
}
for( k = 0; k<n ; k++){
x[n][k] = A[n-1][n+k]/A[n-1][n-1];
for( i = n-1; i>=0; i--){
s = A[i][n+k];
for( j = i+1 ; j<n; j++)
s = s - A[i][j]*x[j][k];
x[i][k] = s/A[i][i];
}
}
}
Inversa en Matlab
function x = inversa(B)
[n dumy] = size(B);
I = eye(n);
A = [B,I];
x = zeros(n);
for i = 1 : n-1
p = i;
while p < n && A(p,i) == 0; p = p + 1; end
if p == n
fprintf('Matriz singular...');
return;
end
if p ~= i
temp = A(i,:);
A(i,:) = A(p,:);
A(p,:) = temp;
end
for j = i+1 : n
m = A(j,i)/A(i,i);
for k = 1 : 2*n
A(j,k) = A(j,k) - m*A(i,k);
end
end
end
Inversa en Matlab
if A(n,n)==0
fprintf('Matriz singular...');
return;
end
for k = 1: n
x(n,k) = A(n,n+k)/A(n,n);
for i = n-1 : -1: 1
s = A(i,n+k);
for j = i+1 : n
s = s - A(i,j)*x(j,k);
end
x(i,k) = s/A(i,i);
end
end
Comandos de Matlab para matrices
inv(m) – invierte una matriz cuadrada m.También m^-1.
/ – divide, A/B es equivalente a B*inv(A)
\ – divide, A\B es equivalente a inv(A)*B
det(m) – calcula el determinante de m.
trace(m) – suma los elementos de la diagonal de m.
Tarea 30
1. Determine cual de las siguientes
matrices son no singulares y calcule su
inversa.
a.| 4 2 6| b.|1 2 0| c.|4 0 0|
| 3 0 7|
|2 1 -1|
|0 0 0|
|-2 -1 -3|
|3 1 1|
|0 0 3|
d.|1
|1
|2
|-1
1 -1 1| e.|4 0 0
2 -4 -2|
|6 7 0
1 1 5|
|9 11 1
0 -2 -4|
|5 4 1
0| f.|2 0 1
0|
|1 1 0
0|
|2 -1 3
1|
|3 -1 4
2|
2|
1|
3|
3
2

 1

4
4
4
5
6
2
0
7
1
6

8
0

1
1
7

0

4
3
4
5
6
0
5
1
2
6

3
5

4
Álgebra lineal
Definición: dos matrices son iguales si tienen el mismo tamaño
n x m y si aij = bij para todo i = 1, 2, .., n y j = 1, 2, ..., m.
Definición: la suma de dos matrices A y B del mismo tamaño n
x m es igual a una matriz C donde cij = aij + bij para todo i = 1,
2, .., n y j = 1, 2, ..., m.
Definición: la multiplicación de una matriz A de tamaño n x m
por una escalar l es una matriz de tamaño n x m cuyos
elementos son laij para todo i = 1, 2, .., n y j = 1, 2, ..., m.
Propiedades de las matrices
Si A, B y C son matrices de tamaño n x m, y l y m son números
reales, se cumplen las siguientes propiedades:
a. A + B = B + A
e. l(A + B) = lB + lA
b. (A + B) + C = A + (B + C)
f. (l + m)A = lA + mA
c. A + 0 = 0 + A = A
g. l (mA) = (lmA
d. A + (-A) = -A + A = 0
h. 1A = A
Multiplicación de matrices
Sea A una matriz de n x m y B una matriz de m x p. El producto
de la matriz A por la matriz B se define como la matriz C de n x p
donde los elementos de C se calculan por:
m
cij   aik bkj  ai1b1 j  ai 2b2 j  ...  aimbmj
k 1
Definiciones
Una matriz cuadrada es la que tiene igual número de renglones
que de columnas.
Una matriz diagonal es una matriz cuadrada con D =(dij) con dij
= 0 simpre que i<>j.
La matriz identidad de orden n, In = (dij) es una matriz diagonal
con dij = 1 si i = j y dij =0 si i<>j.
Una matriz triangular superior de n x n U = (uij) tiene, para
toda j = 1, 2, ..., n, los elementos
uij = 0, para cada i = j +1, j + 2, ..., n
Una matriz triangular inferior de n x n L = (lij) tiene, para toda
j = 1, 2, ..., n, los elementos
lij = 0, para cada i = 1, 2, ..., j – 1
Traspuesta
La traspuesta de una matriz A de n x m es una matriz At tal que
la i_ésima columna de At es la misma que el i-ésimo renglón de
A.
Se cumple que
a. (At)t = A
b. (A + B)t = At + Bt
c. (AB)t = Bt At
d. Si A-1 existe, (A-1)t = (At)-1
Determinante
Si A es una matriz de 1 x 1, entonces det(A) = x.
Si a es de orden mayor que 1, calcule el determinante de
A como sigue:
Escoja cualquier fila o columna. para cada elemento A[i,
j] en esa fila o columna, forme el producto
(-1)(i+ j)*A[i, j]*det(menor(A[i, j]))
Donde det(menor(A[i, j])) es el determinante del menor
de A[i, j].
det(A) = suma de todos los productos para la columna o
fila seleccionada.
Propiedades
a. Si un renglón o columna tiene solo ceros, el determinante es cero.
b. Si se intercambian 2 renglones o columnas, el signo del determinante cambia
c. Si dos columnas o renglones son iguales, el determinante es cero.
d. Si se multiplica un renglón o columna por un numero real el determinante se
multiplica por ese número real.
e. Si se suma un múltiplo de un renglón o columna a otro renglón o columna, el
determinante no se altera.
f. El determinante de un producto de matrices es igual al producto de los
determinantes de cada una.
g. El determinante de la inversa es el inverso del determinante de la matriz
original.
Cálculo eficiente del
determinante
Convertir la matriz en una matriz triangular superior y
calcular el determinante de ésta.
El determinante de una matriz triangular superior es igual al
producto de los elementos de la diagonal principal.
La matriz triangular se obtiene mediante la eliminación
gaussiana.
Determinante en Matlab
function x = determinante(A)
[n dumy] = size(A);
for i = 1 : n-1
p = i;
while p < n && A(p,i) == 0;
p = p + 1;
end
if p == n
x = 0;
return;
end
if p ~= i
temp = A(i,:);
A(i,:) = A(p,:);
A(p,:) = temp;
end
for j = i+1 : n
m = A(j,i)/A(i,i);
for k = 1 : n
A(j,k) = A(j,k) - m*A(i,k);
end
end
end
if A(n,n)==0
x = 0;
return;
end
p = 1;
for i = 1:n
p = p*A(i,i);
end
x = p;
Descomposición LU
La descomposición LU transforma una matriz a en el producto de dos matrices
triangulares
A = LU
(1)
Donde L es triangular inferior y U es triangular superior.
Un sistema Ax = b se puede escribir como
LUx = b
Si hacemos
Ux = d
Obtenemos
Ld = b
(2)
(3)
(4)
Con (4) se obtiene z y después con (3) se obtiene x.
En los casos en que se tengan varios sistemas con la misma matriz de
coeficientes es más eficiente aplicar este método que resolver cada uno
mediante Gauss.
Supongamos un sistema de 3x3
 a11
a
 21
 a31
a12
a22
a32
a13   x1   b1 
a23   x2   b2 
a33   x3  b3 
En el primer paso de eliminación gaussiana se multiplica el renglón 1 por
f 21 = a21/a11 y se le resta al segundo para eliminar a21.
Luego se multiplica el renglón 1 por f 31 = a31/a11 y se le resta al tercer para
eliminar a31.
El paso final es multiplicar el segundo renglón por f 32 = a’32/a’22 y restar el
resultado al tercer renglón.
Los valores de f 21, f 31 y f 32 se pueden retener en a21, a31 y a32 .
Después de la eliminación
 a11
f
 21
 f 31
a12
a '22
f 32
a13 
a '23 
a ' '33 
Esta matriz representa un almacenamiento eficiente de la
descomposición LU de A.
A=LU
donde
a11
U   0
 0
a12
a '22
0
1
L   f 21
 f 31
a13 
a '23 
a ' '33 
0
1
f 32
0
0
1
Ejemplo
 3  0. 1  0 . 2 
A   0.1
7
 0.3
0.3  0.2 10 
Los valores de f 21, f 31 son
f 21 = 0.1/3 = 0.033333 y
f 31 = 0.3/3 = 0.1
 0.1
 0.2 
 3
0.03333 7.00333  0.29333


 0.1
 0.19
10.020 
El valore de f 32 es
f 32 = 0.19/7.003333 = 0.027130
 0.1
 0.2 
 3
0.03333 7.00333  0.29333


 0.1
 0.027130 10.0120 
 0.1
 0.2 
3
U  0 7.00333  0.29333
0
0
10.0120 
0
0
 1
L  0.03333
1
0
 0.1
 0.027130 1
0
0 3
 0.1
 0.2 
 1
A  0.03333
1
0 0 7.00333  0.29333 
 0.1
 0.027130 1 0
0
10.0120 
3
 0.1
 0.2 

0.099999 7


0
.
3


 0.3
 0.2 9.99996
Algoritmo
Sub Decompose(a, n)
DoFor k=1, n – 1
DoFor i=k+1, n
factor = ai,k/ak,k
ai,k = factor
DoFor j=k+1, n
ai,j = ai,j – factor * ak,j
EndDo
EndDo
EndDo
End Decompose
Solución del sistema
Para resolver un sistema de ecuaciones procedemos
multiplicando L por un vector d y sustituyendo hacia adelante.
1 0
l
 21 1
l31 l32
0  d1   b1 
0 d 2   b2 
1  d 3  b3 
d1  b1
l21d1  d 2  b2
l31d1  l32 d 2  d 3  b3
Con estos valores evaluamos U por X y sustituimos hacia atrás
a11
0

 0
a12
a '22
0
a13   x1   d1 
a '23   x2   d 2 
a '33   x3   d 3 
x3  d 3
a '22 x2  a '23 x1  d 2
a11 x1  a12 x2  a13 x3  d1
Ejemplo
 3  0.1  0.2  x1   7.85 
 0. 1
  x    19.3
7

0
.
3

 2  

0.3  0.2 10   x3   71.4 
0
0  d1   7.85 
d1  7.85
d1  7.85
 1
Ld  0.03333
1
0 d 2    19.3
0.3333d1  d 2  19.3
d 2  19.5617
 0.1
 0.027130 1  d 3   71.4  0.1d1  0.02713d 2  d 3  71.4 d 3  70.0843
 0.1
 0.2   x1   7.58 
x1  3
3
Ux  0 7.00333  0.29333  x2    19.5617 x2  2.5
0
0
10.0120   x3   70.0843  x3  7.00003
Sea el sistema:
-6x1 + x2
4x1 + 3x2
5x1 - 3x2
Ejemplo
+ x3 = -37
+ x3 = -25
+ x3 = -34
La descomposición LU nos da
 6 1 1
A   4
3 1
 5  3 1
0
0
 1
L   .6667
1
0
  .8333  .5909 1
0
0  z1   37
 1
Ly   .6667
1
0  z 2    25
  .8333  .5909 1  z3   34
  37
b    25
  34
1
1 
 6
U   0 3.6667 1.6667
 0
0
2.8182
z1 = -37
-.6667z1 + z2 = -25
-.8333z1 - .5909z2 + z3 = -34
1
1   x1    37  2.8182x3 = -94.180
 6
Ux   0 3.6667 1.6667  x2    49.667 3.6667x1 + 1.6667x2 = -49.667
 0
0
2.8182  x3    94.180 -6x1 + x2 + x3 = -37
X1 = 0.87097
X2 = 1.64516
X3 = -33.41935
Pivoteo
La descomposición LU se implementa con el algoritmo de Gauss. Si se modifica
el orden de las filas, se deberá tomar en cuenta este hecho.
El algoritmo siguiente genera una matriz P de permutación para tomar en cuenta
el intercambio de filas.
La matriz de permutación se obtiene intercambiando renglones en una matriz
identidad. Por ejemplo, la siguiente matriz intercambia el primer y segundo de
una matriz A al premultiplicarla por A.
El guión regresa
0 1 0 
P  1 0 0
0 0 1
PA = LU
La ecuación original se escribe
PAx = Py
O
LUx = y
Comandos de Octave
[l, u, p] = lu (a) – regresa la matriz triangular inferior (l), la matriz
triangular superior (u) y la matriz de permutaciones (p)
[l, u] = lu (a) – regresa la matriz triangular inferior (l) y la matriz
triangular superior (u), l se regresa permutada (pl)
Algoritmo en Matlab
function [l u p]= descomposicionLU(a)
%l – matriz triangular inferior
%u – matriz triangular superior
%p – matriz de permutaciones
lu = a;
[n muda] = size(a);
p = eye(n);
psigno = 1;
for k=1:n
j = k;
for i = k+1:n
% encontrar pivote
if abs(lu(i,k))>abs(lu(j,k))
j = i;
end
end
if j!= k
% intercambia con el pivote
t = lu(j,:); lu(j,:) = lu(k,:); lu(k,:) = t;
t = p(j,:); p(j,:) = p(k,:); p(k,:) = t;
psigno = -psigno;
end
if lu (k,k) != 0
for i = k+1:n
lu(i,k) = lu(i,k)/lu(k,k);
for m=k+1:n
lu(i,l) = lu(i,m)-lu(i,k)*lu(k,m);
end
end
end
end %fin lazo k
for i = 1:n
for j = 1:n
if i>j
l(i,j) = lu(i,j);
elseif i==j, l(i,j) = 1; elsel(i,j) = 0;
end
end
end
for i = 1:n
for j = 1:n
if i<=j, u(i,j) = lu(i,j); else u(i,j) = 0;
end
end
end
Substitución hacia adelante y
hacia atrás
function y=sub_adelante(L,w)
[nRow,nCol]=size(L);
y=zeros(nRow,1);
y(1)=w(1)/L(1,1);
for n=2:nRow
y(n)=( w(n) - L(n,1:n-1)*y(1:n-1) ) / L(n,n);
end
function x=sub_atras(U,v)
[nRow,nCol]=size(U);
x=zeros(nRow,1);
x(nRow)=v(nRow)/U(nRow,nRow);
for n=(nRow-1):-1:1
x(n)=(v(n)-(U(n,n+1:end)*x(n+1:end))) / U(n,n);
end
Guión para resolver un sistema
function x=resolver(A,b)
% Resuelve Ax=b
%ENTRADA:
% A: matriz cuadrada
% b: vector de tamaño n
%
%SALIDA:
% x: la solución a A*x=b
[L,U,P] = descomposicionLU(A);
y = sub_adelante(L,P*b);
x = sub_atras(U,y);
Tarea 31
1. Factorice las siguientes matrices en la
descomposición LU con Octave y con los guiones para LU.
Haga a mano el producto PLU y confirme que es igual a
la matriz original. Utilice 3 cifras después del punto
decimal.
a.| 4 2 6| b.|1 2 6| c.|1 3 8|
| 3 5 7|
|2 1 -1|
|0 7 0|
|-2 -1 3|
|3 1 1|
|1 7 3|
d.| 1 -1 -1 1| e.|4 0 4 0| f.|1 3 1 2|
| 1 2 -4 -2|
|6 7 1 5|
|8 1 -6 2|
| 2 1 1 5|
|7 1 1 0|
|2 -4 3 -1|
|-1 6 -2 -4|
|2 4 1 1|
|3 -2 4 3|
2. Encuentre los determinantes de las matrices
anteriores
Aplicación
Una red eléctrica puede representarse por medio de un sistema de ecuaciones.
Supongamos una red que consista solo de resistencias y fuentes de voltaje
independientes.
El sistema a resolver es
RI = V
Donde R es la matriz de resistencias de la red, V es un vector donde cada
componente es el voltaje de la malla e I es el vector de corrientes a determinar.
La matriz de resistencias tiene la siguiente forma
 R11
 R
 21
 .

  Rn1
 R12
R22
.
.
...  R1n 
...
. 
.
. 

.
Rnn 
Donde Rii es la resistencia total de la malla i, y –Rij es la resistencia total entre la
malla i y la malla j.
Análisis de redes
R2
R1
R3
R5
I1
R4
V3
 R1  R2  R4

 R5


 R4

0

V1
I2
R6
V2
I3
I4
R7
R8
 R5
 R4
R3  R5  R6
0
0
R4  R7
 R6
0
R9
  I1   0 
 I    V 
 R6
1 
 2   
  I 3  V3  V2 
0
  

R6  R8  R9   I 4   V2 
0
Solución en Matlab
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
Calculo de las corrientes de malla en una red eléctrica
a - vector utilizado para leer los datos desde un archivo
nm = numero de mallas
nr = numero de resistencias
m1 = vector con la malla uno de cada resistencia
m2 = vector con la malla dos de cada resistencia
res = vector con los valores de las resistencias leídas
r = matriz de resistencias
v = vector con los voltajes totales de malla
amp = vector de corrientes de malla
la descripción de la red se lee desde un archivo, el formato es
nm
nr
res(1)
m1(1)
m2(1)
...
res(nr)
m1(nr)
m2(nr)
v(1)
...
v(nm)
Solución en Matlab
nn = 1;
%almacena todos los datos en el vector a
a = textread('malla.txt');
nm = a(nn);nn = nn+1;
nr = a(nn);nn = nn+1;
%lee cada elemento
for n = 1:nr
res(n) = a(nn);nn = nn+1;
m1(n) = a(nn);nn = nn+1;
m2(n) = a(nn);nn = nn+1;
end
%descripcion de los datos leidos
for n = 1:nr
fprintf('R%d=%6.3f entre malla %d y %d\n',n,res(n),m1(n),m2(n));
end
%Leer voltajes de malla y los despliega
for n = 1:nm
v(n) = a(nn);nn = nn+1;
fprintf('Voltaje de malla %d = %12f8:\n',n,v(n))
end
Solución en Matlab
%crea matriz de resistencias
r = zeros(nm);
for n = 1:nr
j = m1(n);
k = m2(n);
if k~=0
r(k,k) = r(k,k) + res(n);
r(j,k) = r(j,k) - res(n);
r(k,j) = r(j,k);
end
r(j,j) = r(j,j) + res(n);
end
%imprime matriz de resistencias
fprintf('Matriz de resistencias');
r
%calcula corrientes
amp = inv(r)*v';
%imprimir resultados
for n = 1:nm
fprintf('Corrientes de malla %d = %12f8:\n',n,amp(n))
end
Ejemplo numérico
Número de mallas
Número de resistencias
38.9
I5
I1
3492
1017
949
Valor de resistencia
837
4821
42.1
I2
I4
986
Malla 1
Malla 2
842
59.76
5732
I3
5670
5
9
1017
1
3
3492
1
2
837
2
4
949
2
0
842
3
5
4821
3
4
986
4
0
5732
4
5
5670
5
0
38.9
0
42.1
0
59.76
Resultados
R1=1017.000 entre malla 1 y 3
R2=3492.000 entre malla 1 y 2
R3=837.000 entre malla 2 y 4
R4=949.000 entre malla 2 y 0
R5=842.000 entre malla 3 y 5
R6=4821.000 entre malla 3 y 4
R7=986.000 entre malla 4 y 0
R8=5732.000 entre malla 4 y 5
R9=5670.000 entre malla 5 y 0
Voltaje de malla 1 = 38.900000
Voltaje de malla 2 = 0.000000
Voltaje de malla 3 = 42.100000
Voltaje de malla 4 = 0.000000
Voltaje de malla 5 = 59.760000
Matriz de resistencias
r=
4509 -3492 -1017
0
0
-3492 5278
0 -837
0
-1017
0 6680 -4821 -842
0 -837 -4821 12376 -5732
0
0 -842 -5732 12244
Corrientes de malla 1 =
Corrientes de malla 2 =
Corrientes de malla 3 =
Corrientes de malla 4 =
Corrientes de malla 5 =
0.036338
0.027346
0.028966
0.020833
0.016626
2  1 1  4 2 3 1 2  1
1 1  1
a) 3 3 9 b)  3 1 4 c) 0 1  1 d ) 0 2 3 
3 3 5  2 4 5 2 4 0 
0  1 1 
2
 1

 0
2
 3

 1
1
0   x1  1
2  1  x2   2
 1 2   x3  3
 1 1   x1  4
4  1  x2   5
 1 1   x3  6