Análisis Numérico para Ingeniería - fi.mdp.edu.ar · Factorización de Doolittle ... de la...

Análisis Numérico para Ingeniería

Clase Nro. 5

Mg. Ing. Francisco A. Lizarralde Facultad de Ingeniería - UNMDP - 2017 2

Sistemas de Ecuaciones Lineales

Tipos de Matrices

Método de Triangulación de Gauss

Método de Diagonalización de Gauss-Jordan

Tácticas de Pivoteo

Inversión de una Matriz

Refinamiento Iterativo de la Solución

Método de Thomas

Reducción de Crout

Temas a tratar:


Tipos de Matrices

Matriz Densa

Matriz Rala

Matriz Triangular Inferior

Matriz Triangular Superior

Matriz Diagonal

Matriz Tri-Diagonal

Matriz de Banda


Matrices Densas

Las matrices densas tienen pocos elementos nulos.


Matrices poco densas (Ralas)

En las matrices ralas o poco densas, la mayoría de sus valores son nulos.


Sistema Triangular Inferior

[a11 0 0 0a21 a22 0 0a31 a32 a33 0a41 a42 a43 a44

]⋅[x1

x2

x3

x4]=[

b1

b2

b3

b4]

x1=b1

a11

; xi=

bi−∑k=1

i−1

aik⋅xk

aii

; i=2,3 ,, n

Su solución se obtiene por sustitución progresiva:

Dado el siguiente Sistema Triangular Inferior


Sistema Triangular Superior

[a11 a12 a13 a14

0 a22 a23 a24

0 0 a33 a34

0 0 0 a44]⋅[

x1

x2

x3

x4]=[

b1

b2

b3

b4]

xn=bn

ann

; xi=

bi− ∑k=i1

n

aik⋅xk

aii

; i=n−1,n−2, ,1

Su solución se obtiene por sustitución regresiva:

Dado el siguiente Sistema Triangular Superior


Sistema Diagonal

[a11 0 0 00 a22 0 00 0 a33 00 0 0 a44

]⋅[x1

x2

x3

x4]=[

b1

b2

b3

b4]

x1=b1

a11

; x2=b2

a22

; x3=b3

a33

; x4=b4

a44

Su solución se obtiene como:

Dado el siguiente Sistema Diagonal


Ejemplo de aplicación

Las plantas con la presencia del sol convierten durante su proceso de fotosíntesis, el dióxido de carbono y el agua en glucosa y oxígeno.

Simbólicamente el proceso puede ser representado así:

Dióxido de carbono( CO2 ) + Agua( H

2O ) → Glucosa( C

6H

12O

6 ) + Oxígeno( O

2 )

Balancear la ecuación química de la reacción, determinando los valores enteros más pequeños para las variables, si denominamos:

x = la proporción dióxido de carbono presente

y = la proporción de agua presente,

w = la proporción de glucosa producida, y

z = la proporción de oxígeno producido.


Planteo del ejemplo de aplicación

De esta forma, a ecuación de la reacción química puede ser representada así:

x CO2 + y H

2O → w C

6H

12O

6 + z O

2

El número de átomos presente en un elemento químico está indicado mediante un subíndice. Por lo tanto, comparando la cantidad de átomos de cada elemento que hay en los reactivos y los productos, obtenemos el siguiente sistema;

x = 6w Para balancear el Carbono

2x + y = 6w + 2z Para balancear el Oxígeno

2y = 12w Para balancear el Hidrógeno

Por lo tanto, reordenando obtenemos el siguiente sistema de ecuaciones.

x - 6w = 0

2x + y - 6w - 2z = 0

2y - 12 w = 0


Notación explícita

a11 x1a12 x2a13 x3a1n xn=b1



an1 x1an2 x2an3 x3ann xn=bn

Notación explícita para un Sistema de Ecuaciones Lineales de la forma: Ax=b


Notación Matricial

[a11 a12 a13 a1n

a21 a22 a23 a2n

a31 a32 a33 a3n

⋮ ⋮ ⋮ ⋱ ⋮an1 an2 an3 ann

]⋅[x1

x2

x3

⋮xn

]=[b1

b2

b3

⋮bn

]

Notación matricial para un Sistema de Ecuaciones Lineales de la forma: Ax=b


Métodos Directos

Métodos de Resolución

Triangulación de Gauss

Diagonalización de Gauss-Jordan

Método de Thomas

Métodos de Factorización

Factorización de Crout

Factorización de Doolittle

Factorización de Choleski


Matriz Ampliada

a11 a12 a13 a1n

a21 a22 a23 a2n

a31 a32 a33 a3n

⋮ ⋮ ⋮ ⋱ ⋮an1 an2 an3 ann

∣b1

b2

b3

⋮bn

La matriz ampliada es una matriz que contiene los elementos de la matriz original junto con los elementos de los términos independientes, colocados como columnas adicionales.


Método de Gauss (I)

a11 a12 a13

a21 a22 a23

a31 a32 a33∣

b1

b2

b3

1 a '12 a '13

0 a '22 a '23

0 a '32 a '33∣

b '1

b '2

b '3

F '1=F1

a11

F '2=F2−F'1⋅a21

F '3=F3−F'1⋅a31

Matriz ampliada del sistema original

Tablero 1: Resultado de las operaciones sobre el sistema original


Método de Gauss (II)

1 a '12 a '13

0 a '22 a '23

0 a '32 a '33∣

b '1

b '2

b '3

Matriz ampliada correspondiente al Tablero 1

Tablero 2: Resultado de las operaciones sobre el Tablero 1

1 a '12 a '13

0 1 a ' '23

0 0 a ' '33∣

b '1b ' '2

b ' '3

→ F ' '2=F'2a '22

→ F ' '3=F'3−F ' '2⋅a '32


Método de Gauss (III)

1 a '12 a '13

0 1 a ' '23

0 0 a ' '33∣

b '1

b ' '2

b ' '3

1 a '12 a '13

0 1 a ' '23

0 0 1 ∣b '1

b ' '2

b ' ' '3→ F ' ' '3=

F' '3a ' '33

Tablero 2



Método de Gauss-Jordan (I)

a11 a12 a13

a21 a22 a23

a31 a32 a33∣

b1

b2

b3

1 a '12 a '13

0 a '22 a '23

0 a '32 a '33|

b '1

b '2

b '3

→ F '1=F1

a11

→ F '2=F2−F '1⋅a21

→ F '3=F3−F '1⋅a31

Matriz ampliada del sistema original

Tablero 1: Resultado de las operaciones sobre el sistema original


Método de Gauss-Jordan (II)

1 a '12 a '13

0 a '22 a '23

0 a '32 a '33∣

b '1

b '2

b '3

Matriz ampliada correspondiente al Tablero 1


1 0 a ' '13

0 1 a ' '23

0 0 a ' '33∣

b ' '1

b ' '2

b ' '3

→ F ' '1=F '1−F' '2⋅a '12

→ F ' '2=F '2a '22

→ F ' '3=F'3−F' '2⋅a '32


Método de Gauss-Jordan (III)

1 0 a ' '13

0 1 a ' '23

0 0 a ' '33∣

b ' '1

b ' '2

b ' '3

1 0 00 1 00 0 1 ∣

b' ' '1b' ' '2

b' ' '3

F' ' '1=F ' '1−F ' ' '3⋅a' '13

F' ' '2=F' '2−F ' ' '3⋅a' '23

F' ' '3=F ' '3a ' '33

Tablero 2



Resolución Simultánea de Sistemas

a11 a12 a13 … a1n

a21 a22 a23 … a2n

a31 a32 a33 … a3n

⋮ ⋮ ⋮ ⋱ ⋮an1 an2 an3 … ann

∣b1 c1 d1 … w1

b2 c2 d2 … w2

b3 c3 d3 … w3

⋮ ⋮ ⋮ ⋱ ⋮bn cn dn … wn

Es posible utilizar los métodos antes mencionados, para resolver n sistemas simultáneamente. Lo único que hay que hacer es agregar a la matriz ampliada los n vectores de términos independientes.


Tácticas de pivoteo

Las tácticas de pivoteo se emplean para evitar la amplificación de los errores producidos por trabajar con aritmética finita.

Por lo tanto se busca que los elementos de mayor valor absoluto, se encuentren en los pivotes.

Las estrategias de pivoteo más utilizadas son el pivoteo parcial por filas, o por columnas y el pivoteo total.


a11 a12 a13 … a1n

a21 a22 a23 … a2n

a31 a32 a33 … a3n

⋮ ⋮ ⋮ ⋱ ⋮an1 an2 an3 … ann

∣100⋮0

La solución del sistema A.x = v1, donde v

1 es el 1er. Versor,

x=A−1⋅v1equivale a :

Es decir, dicha solución es la primer columna de la inversa de la matriz A



a11 a12 a13 … a1n

a21 a22 a23 … a2n

a31 a32 a33 … a3n

⋮ ⋮ ⋮ ⋱ ⋮an1 an2 an3 … ann

∣010⋮0

La solución del sistema A.x = v2, donde v

2 es el 2do. Versor,

x=A−1⋅v2equivale a :

Es decir, dicha solución es la segunda columna de la inversa de la matriz A




a11 a12 a13 a1n

a21 a22 a23 a2n

a31 a32 a33 a3n

⋮ ⋮ ⋮ ⋱ ⋮an1 an2 an3 ann

∣1 0 0 00 1 0 00 0 1 0⋮ ⋮ ⋮ ⋱ ⋮0 0 0 1

Para hallar la inversa de una matriz por medio del método de Gauss-Jordan. Simplemente se ingresan tantos versores como sea el orden de la matriz y se aplica el método a todos los sistemas simultáneamente. Como la solución de cada sistema corresponde a cada una de las columnas de la matriz inversa, al finalizar el proceso de diagonalización, se obtendrá la matriz inversa completa.



SUBROUTINE Inversa(matriz)! Calcula la inversa de una matriz, por Gauss-JordanREAL(8), DIMENSION(:,:), INTENT(IN) :: matrizREAL(8), ALLOCATABLE, DIMENSION(:,:) :: matrizIdentidadINTEGER orden

orden = SIZE(matriz, DIM=1)

CALL creaMatrizIdentidad(orden, matrizIdentidad) CALL GaussJordan(matriz, matrizIdentidad) DEALLOCATE(matrizIdentidad) END SUBROUTINE

Implementación con Gauss-Jordan


Refinamiento Iterativo (I)

La solución exacta de un sistema de ecuaciones lineales es:

A⋅x−b=0

A⋅x*−b=r≠0

A⋅x*−b=r

-A⋅x−b=0

----------------

A⋅(x*−x)=r ⇒ A⋅(Δ x)=r

Como solo podemos calcular una solución aproximada x* , nos queda un residuo r.

Si restamos ambas expresiones, obtenemos :


Refinamiento Iterativo (II)

Con este nuevo valor de x, podemos calcular un nuevo residuo y compararlo con una cota establecida. Si dicha cota no se ha alcanzado, se repite el proceso.

xi+1=xi−Δ xi

A⋅Δ x=r

Por lo tanto podemos hallar la solución del siguiente sistema:

Esta solución a su vez, nos permite calcular un nuevo valor xi+1

Sin embargo, como no conocemos el valor exacto, tenemos:

A⋅(x0−x1)=r0A⋅Δ x0=r0


SI

NO

|| r || > ε

HACER x* = x* - Δx

CALCULAR Δx resolviendo A.Δx = r

CALCULAR r comor = Ax* - b

CALCULAR x* resolviendo A.x* = b

Refinamiento Iterativo (III)


Método de Thomas

El método de Thomas es una optimización del método de Gauss, para resolver sistemas cuyas matrices son tri-diagonales.

[d1 u1 0 0 0l1 d2 u2 0 00 l2 d3 u3 00 0 l3 d4 0⋮ ⋮ ⋮ ⋮ ⋱ ⋮0 0 0 dn−1 un−1

0 0 0 ln−1 dn

]⋅[x1

x2

x3

x4

⋮xn−1

xn

]=[b1

b2

b3

b4

⋮bn−1

bn

]


u = . . . 0u(1) u(2) u(3) u(4) . . . u(n)

d = . . .d(1) d(2) d(3) d(4) . . . d(n)

l = 0 . . .l(1) l(2) l(3) l(4) . . . l(n)

u1

u2

u3

u4

d1

d2

d3

d4

dn

l1

l2

l3

ln-1

Diagonal Superior

Diagonal Principal

Para trabajar con vectores del mismo tamaño, a los vectores que almacenan los valores de las diagonales secundarias, se los completa con un 0, según corresponda. Esto simplifica enormemente el algoritmo.

Diagonal Inferior

Variables en el método de Thomas


SUBROUTINE Thomas(u_orig, d_orig, l_orig, term_indep)! Metodo de Thomas para matrices Tri-DiagonalesREAL(8), DIMENSION(:,:), INTENT(IN) :: term_indepREAL(8), DIMENSION(:), INTENT(IN) :: u_orig, d_orig, l_origREAL(8), DIMENSION(:,:), ALLOCATABLE :: bREAL(8), DIMENSION(:), ALLOCATABLE :: u, d, l INTEGER orden, cant_vec, i orden = SIZE(term_indep, DIM=1)cant_vec = SIZE(term_indep, DIM=2)

! Realiza copias del originalALLOCATE(u(orden), d(orden), l(orden))u = u_origd = d_origl = l_origALLOCATE(b(orden, cant_vec))b = term_indep Sique en la

próxima página

Implementación del método de Thomas (I)


! Aqui comienza el algoritmoDO i=1, orden-1 u(i) = u(i) / d(i) b(i,:) = b(i,:) / d(i) d(i) = 1.0 d(i+1) = d(i+1) - l(i+1)*u(i) b(i+1,:) = b(i+1,:) - l(i+1)*b(i,:) l(i+1) = 0.0ENDDO

! Obtencion de la solucion por Sustitucion Inversa b(orden,:) = b(orden,:)/d(orden) DO i=orden-1, 1, -1 b(i,:) = b(i,:) - u(i)*b(i+1,:) / d(i) END DO

CALL imprimeMatriz(b) DEALLOCATE(u, d, l, b)END SUBROUTINE

Implementación del método de Thomas (II)


La factorización LU (o factorización triangular) de una matriz A, consiste en reescribir dicha matriz, como producto de dos matrices triangulares L y U. Donde la matriz L es triangular inferior y la matriz U es triangular superior.

Los métodos de factorización triangular más conocidos son:

La Reducción de Crout ( LU1 )

El método de Doolittle ( L1U )

Método de Choleski para matrices A positivas definidas.

Las matrices L1 y U

1 son matrices triangulares que poseen

valores unitarios en la diagonal.

Métodos de Factorización Triangular


Métodos de Factorización

Factorización de CroutSe especifica u

11= u

22= u

33= … = 1

Factorización de DoolittleSe especifica l

11= l

22= l

33= … = 1

Factorización de CholeskiSe especifica l

ii= u

ii

No los veremos en el presente curso


Factorización de Crout (LU1)

[l11 0 0 0l21 l22 0 0l31 l32 l33 0l41 l42 l43 l44

]⋅[1 u12 u13 u14

0 1 u23 u24

0 0 1 u34

0 0 0 1]=[

a11 a12 a13 a14

a21 a22 a23 a24

a31 a32 a33 a34

a41 a42 a43 a44]

Si multiplicamos las filas de L por la 1er. columna de U1 podemos

observar que los elementos de la 1er. columna de L son iguales a los elementos de la 1er. columna de A.

l11=a11 ; l21=a21 ; l31=a31 ; l41=a41

Elementosa calcular


Factorización de Crout (LU1) (I)

[l11 0 0 0l21 l22 0 0l31 l32 l33 0l41 l42 l43 l44

]⋅[1 u12 u13 u14

0 1 u23 u24

0 0 1 u34

0 0 0 1]=[

a11 a12 a13 a14

a21 a22 a23 a24

a31 a32 a33 a34

a41 a42 a43 a44]

Si multiplicamos la primera fila de L por las columnas de U1 obtenemos:

l11=a11 ; l11⋅u12=a12 ; l11⋅u13=a13 ; l11⋅u14=a14

Por lo tanto podemos calcular los valores de la primera fila de U1 como:

u12=a12

l11

; u13=a13

l11

; u14=a14

l11

Elementosa calcular


Factorización de Crout (LU1) (II)

[l11 0 0 0l21 l22 0 0l31 l32 l33 0l41 l42 l43 l44

]⋅[1 u12 u13 u14

0 1 u23 u24

0 0 1 u34

0 0 0 1]=[

a11 a12 a13 a14

a21 a22 a23 a24

a31 a32 a33 a34

a41 a42 a43 a44]

Si multiplicamos las filas desde 2 hasta n, de L por las 2da. columna de U1 obtenemos:

l21⋅u12l22=a22 ; l31⋅u12l32=a32 ; l41⋅u12l42=a42

Por lo tanto podemos calcular los valores de la segunda columna de L1 como:

l22=a22−l21⋅u12 ; l32=a32−l31⋅u12 ; l42=a42−l41⋅u12

Elementosa calcular


[l11 0 0 0l21 l22 0 0l31 l32 l33 0l41 l42 l43 l44

]⋅[1 u12 u13 u14

0 1 u23 u24

0 0 1 u34

0 0 0 1]=[

a11 a12 a13 a14

a21 a22 a23 a24

a31 a32 a33 a34

a41 a42 a43 a44]

Repitiendo el mismo esquema, obtenemos:

l21⋅u13l22⋅u23=a23 ; l21⋅u14l22⋅u24=a24

Por lo tanto podemos calcular los valores de la segunda fila de U1 como:

u23=a23−l21⋅u13

l22

; u24=a24−l21⋅u14

l22

Factorización de Crout (LU1) (III)


Factorización de Crout (LU1) (IV)

[l11 0 0 0l21 l22 0 0l31 l32 l33 0l41 l42 l43 l44

]⋅[1 u12 u13 u14

0 1 u23 u24

0 0 1 u34

0 0 0 1]=[

a11 a12 a13 a14

a21 a22 a23 a24

a31 a32 a33 a34

a41 a42 a43 a44]

Si multiplicamos las filas desde 3 hasta n, de L por las 3ra. columna de U1 obtenemos:

l31⋅u13l32⋅u23l33=a33 ; l41⋅u13l42⋅u23l43=a43

Por lo tanto podemos calcular los valores de la tercera columna de L1 como:

l33=a33−l31⋅u13−l32⋅u23 ; l43=a43−l41⋅u13−l42⋅u23

Elementosa calcular


[l11 0 0 0l21 l22 0 0l31 l32 l33 0l41 l42 l43 l44

]⋅[1 u12 u13 u14

0 1 u23 u24

0 0 1 u34

0 0 0 1]=[

a11 a12 a13 a14

a21 a22 a23 a24

a31 a32 a33 a34

a41 a42 a43 a44]

Repitiendo el mismo esquema, obtenemos:

l31⋅u14l32⋅u24l33⋅u34=a34

Por lo tanto podemos calcular los valores de la segunda fila de U1 como:

u34=a34−l31⋅u14−l32⋅u24

l33

Elementosa calcular

Factorización de Crout (LU1) (V)


Factorización de Crout (LU1) (VI)

[l11 0 0 0l21 l22 0 0l31 l32 l33 0l41 l42 l43 l44

]⋅[1 u12 u13 u14

0 1 u23 u24

0 0 1 u34

0 0 0 1]=[

a11 a12 a13 a14

a21 a22 a23 a24

a31 a32 a33 a34

a41 a42 a43 a44]

Si multiplicamos la filas 4, de L por las 4ta. columna de U1 obtenemos:

l41⋅u14l42⋅u24l43⋅u34l44=a44

Por lo tanto podemos calcular el último valor de L1 como:

l44=a44−l41⋅u14−l42⋅u24−l43⋅u34

Elementosa calcular


Dado que la matriz A se ha transformado en el producto L.U1

A⋅x=b L⋅U 1⋅x=b L⋅U 1⋅x=b

L⋅U 1⋅x =b L⋅c=b

Si llamamos c a (U1.x) , el sistema se transforma en:

Como la matriz L es una matriz triangular inferior, podemos obtener los valores de c, aplicando sustitución progresiva.

c1=b1

l11

; ci=

bi−∑k=1

i−1

lik⋅ck

lii; i=2,3 , , n

Solución a partir de Crout (I)


Una vez obtenidos los valores de c, como:

L⋅U 1⋅x =b L⋅c=b

Y como la matriz U1 es una matriz triangular superior, podemos

obtener los valores de x, aplicando sustitución regresiva.

U 1⋅x=c

Nos queda el siguiente sistema:

xn=cn ; xi=ci− ∑k=i1

n

uik⋅xk ; i=n−1,n−2, ,1

Solución a partir de Crout (II)


PREGUNTAS ...

Análisis Numérico para Ingeniería - fi.mdp.edu.ar · Factorización de Doolittle ... de la...

Documents

Transcript of Análisis Numérico para Ingeniería - fi.mdp.edu.ar · Factorización de Doolittle ... de la...