Taller 3 Análisis Numérico

Resultado de aprendizaje

Explica la necesidad de métodos numéricos, la condición de problemas tal como la estabilidad de

algoritmos, y todas las aproximaciones involucradas en el proceso de resolución computacional de
problemas de ingenierı́a, obteniendo soluciones numéricas adecuadas y controlando los errores cometi-
dos.

Contenido temático

Métodos numéricos para sistemas lineales. Normas vectoriales y matriciales, condicionamiento de una
matriz. Matrices triangulares. Factorización LU. Cholesky, Jacobi y Gauss-Seidel.

Definiciones y contenidos previas

i) Definiremos el número de condición de A:

cond(A) = ∥A∥ · ∥A−1 ∥

ii) Los métodos directos entregan una solución exacta en un número finito de pasos; sin embargo,
son ineficientes en matrices grandes o mal condicionadas.

iii) Los métodos iterativos entregan soluciones aproximadas, pero son más eficientes para matrices
grandes o mal condicionadas.

Métodos

1. Métodos directos:

1.1 Factorización LU.

1.2 Factorización de Cholesky.

2. Métodos iterativos:

2.1 Jacobi.
2.2 Gauss-Seidel.
 Actividad 1

Consideremos el sistema Ax = b dado por:

     
1 1 x 2
· =
1 0.0001 y 2.0001

a) Obtenga la inversa de A y luego mulitpı́quela por b para obtener una solución al sistema, llámela
x1

b) Considere ahora b2 = (2, 2.0002)t y resuelva nuevamente el sistema.

c) Compare la soluciones.

De manera similar, esta diferencia la podemos ver como:

d) Considere la aproximación x2 = (1, 1.1)t de la solución del sistema original. Esta solución parece
razonable; para argumentar esto, calcule ||x1 − x2||.

e) Finalmente, calcule A · x2. ¿Le parece esta una buena aproximación de la solución?

f) Calcule el número de condición de A y justifique las diferencias antes encontradas.

Funciones y operadores a utilizar

• Para obtener A−1 · b :

A\b

• Para calcular el número de condición de A:

cond(A)

• Para calcular la norma de un vector ||v||:

norm(v)
 
a b
• Para ingresar una matriz :
c d

[a, b; c, d];

• El producto de matrices (o matriz x vector) se hace simplemente con el sı́mbolo ∗

 Método 1: Factorización LU

Este método consiste consta de dos partes:

• Factorización: Consiste en descomponer la matriz A (o, de ser necesario, habiendo permutado

filas previamente a A) como el producto de dos matrices triangulares; superior e inferior.

• Resolución del sistema: Utilizando la factorización encontrada, resolver dos sistemas más sen-
cillos, y menos costosos computacionalmente, que el original.

Factorización

Si A es una matriz cuadrada, invertible, podemos escribirla como producto de una matriz triangu-
lar inferior (L) y una matriz triangular superior (U ) y, si es necesario, incluyendo una matriz de
permutaciones al comienzo:
P ·A=L·U

Actividad 2

Considere la matriz A  
1 2 3
A = 4 5 6 
7 8 9

a) Obtener la descomposición LU de A obteniendo, además, la matriz de permutación de A que

permite hacer esta descomposición.

b) Obtenga explı́citamente P, L, U .

c) Usando las matrices entregadas: L, U, P , verificar que P · A = L · U

Funciones y operadores a utilizar

• Para hacer la descomposición LU de A:

[L,U,P] = lu(A)

Esto automáticamente crea las matrices L, U y P de A.

 Resolución de sistemas

Una vez que tenemos la descomposición P · A = L · U , un sistema de la forma A · x = b se puede

reescribir como:
P −1 · L · U · x = b
o, de manera equivalente,
L·U ·x=P ·b
haciendo el cambio de variables, y = U · x, tenemos:

L·y =P ·b

Observemos que, al ser L una matriz triangular inferior, es posible resolver fácilmente obteniendo los
valores de y desde arriba hacia abajo.

Una vez obtenida la solución y de este nuevo sistema, se resuelve U · x = y el cual, nuevamente es
un sistema sencillo de resolver puesto que U es ahora triangular superior, por lo que la resolución del
sistema será de abajo hacia arriba.

Actividad 3

Resuelva el sistema:      
0 2 1 x 4
1 1 0 · y  = 2
2 1 1 z 6

Por el método P A = LU indicando las matrices P, L, U y los vectores y y la solución x.

No es necesario hacer a mano la descomposición LU, comience con lu(A) y luego, utilizando esas
matrices, escriba el código para resolver el sistema.

Funciones y operadores a utilizar

Esta actividad se puede realizar sin funciones u operadores diferentes a los ya vistos.
 Algoritmo en pseudo lenguaje para la resolución de sistemas

En este algoritmo se presenta cómo resolver un sistema “a mano ” a partir de la descomposición

P A = LU entregada por matlab, sin embargo es posible hacerlo directamente invirtiendo las matrices
L y U.

A ← matriz (n x n)
b ← vector (n x 1)
Hacer la descomposición PA=LU obteniendo la matriz L con unos en la diagonal:
L, U, P] = lu(A)
Establecer Pb ← P*b
Establecer n ← largo(b)
Inicialiar y ← vector de n x 1 con solo ceros
Para i desde 1 hasta n
y(i) ← Pb(i) - Suma de j desde 1 hasta i-1 L(i,j)*y(j)
FinPara
Para i desde n hasta 1
x(i) ← (y(i) - Suma de j desde i+1 hasta n U(i,j)*x(j))/U(i,i)
FinPara
Mostrar "Solución final "x""
Verificar Ax=b
 Método 2: Factorización de Cholesky

La descomposición de Cholesky consiste en encontrar una matriz L, triangular inferior, tal que:

A = L · Lt

Dado que L·Lt es tanto simétrica como definida positiva, la matriz A debe cumplir ambas condiciones.

La forma de resolver un sistema Ax = b utilizando la descomposición de Cholesky es similar al utilizado

en la descomposición P A = LU , es decir:

A·x = b
L · Lt · x = b
L·y = b

Obteniendo el valor de y, resolvemos para x:

Lt · x = y

Actividad 4

Resolver el siguiente sistema utilizando la descomposición de Cholesky indicando, además, la matriz

utilizada L.
     
6 3 4 x 1
3 6 5  · y  = 2
4 5 10 z 3

Funciones y operadores a utilizar

• Para obtener la descomposición de Cholesky de una matriz A (simétrica y definida positiva):

chol(A,’lower’)

Observación: Esta función calcula L en la descomposición de Cholesky de la matriz A. Si se

omite el argumento ’lower’, en vez de entregar L entregará Lt .

Algoritmo en pseudo lenguaje

Puede reutilizar y adecuar el algoritmo y el código de la resolución por la descomposición LU realizado

en la actividad anterior o, simplemente utilizar L y Lt entregadas y calcular sus inversas.
 Métodos 3 y 4: Jacoby y Gauss-Seidel

Tanto el método de Jacoby como el de Gauss-Seidel son métodos iterativos, por lo que se obtiene
una aproximación de la solución según un nivel de tolerancia previamente definido. En este caso,
utilizaremos la norma infinito. En cada paso de xn → xn+1 calcularemos ||xn − xn+1 ||∞ y, si esta
distancia es menor a la tolerancia, detendremos la iteración.
Estos métodos requieren que la matriz A del sistema lineal sea diagonalmente dominante o definida
positiva.

Descripción de los métodos

Ambos métodos se basan en descomponer la matriz A como:

A=M −N

e iterar:
M · xk+1 = N · xk + b
Si ρ(M −1 N ) < 1 entonces el método converge.

• En el método de Jacoby, se utiliza M = D (D=diag(A)) y N = −(L + U ).

• En el método de Gauss-Seidel, se utiliza M = D + L y N = −U

Actividad 5

• Resolver el siguiente sistema utilizando primero el método de Jacoby y luego con el método de
Gauss-Seidel (hacerlo, en ambos casos, para distintos valores de tolerancia y cantidad máxima
de iteraciones).      
10 −1 2 0 x1 6
−1 11 −1 3  x2   25 
 2 −1 10 −1 · x3  = −11
     

0 3 −1 8 x4 15
 Algoritmo en pseudo lenguaje - Jacoby

Entrada:
A → matriz (n x n)
b → vector (n x 1)
x → vector inicial (por ejemplo, todos ceros)
tol → tolerancia para detener el método
maxIter → número máximo de iteraciones
Inicio del algoritmo:

xNuevo ← vector de ceros del mismo tama~no que b

k ← 0 (contador de iteraciones)
Mientras k < maxIter hacer:
Para i desde 1 hasta n hacer:
suma ← 0
Para j desde 1 hasta n hacer:
Si j != i entonces:
suma ← suma + A[i][j] * x[j]
Fin Para
xNuevo[i] ← (b[i] - suma) / A[i][i]
Fin Para
Si normaInfinita(xNuevo - x) < tol entonces:
Terminar programa
x ← xNuevo
k ← k + 1
Fin Mientras
Salida:
xNuevo → solución aproximada
k → número de iteraciones realizadas
 Algoritmo en pseudo lenguaje Gauss-Seidel

A → matriz de coeficientes (n x n)
b → vector independiente (n x 1)
x → vector inicial (por ejemplo, todos ceros)
tol → tolerancia para detener el método
maxIter → número máximo de iteraciones
Inicializar:
xAnterior ← copia de x
k ← 0 (contador de iteraciones)
Mientras k < maxIter hacer:
Para i desde 1 hasta n hacer:
suma ← 0
Para j desde 1 hasta n hacer:
Si j != i entonces:
suma ← suma + A[i][j] * x[j]
Fin Para
x[i] ← (b[i] - suma) / A[i][i]
Fin Para
Si normaInfinita(x - xAnterior) < tol entonces:
Terminar programa
xAnterior ← copia de x
k ← k + 1
Fin Mientras
Salida:
x → solución aproximada
k → número de iteraciones realizadas

Indicaciones y sugerencias

• Utilice a libre disposición Google o algún asistente de inteligencia artificial de forma responsable
y útil; pregunte cómo definir una función, cómo derivarla, etc. No lo utilice para crear el código
automáticamente; un objetivo de este laboratorio es reforzar los contenidos de clases viendo los
resultado sin necesidad de hacer cálculos a mano.

• Si tienen inquietudes, no duden en llamar a su ayudante.

• Si guardan variables y corre nuevamente su código, recuerden usar clear para que, en caso
de haber algún error en la redefinición de una variable, este se informe y no utilice los valores
anteriores.