El método de Jacobi también usa el vector inicial X0, la diferencia consiste en que la actualización del vector X en cada iteración se realiza cuando se ha calculado el vector nuevo completo.
Tarea
Realice el algoritmo de Jacobi en Python, las modificaciones al algoritmo de Gauss-Seidel, incluyendo el pivoteo por filas de la matriz.
El ejemplo de referencia, ya presenta una matriz pivoteada por filas, por lo que no fué implementado el procedimiento. Para generalizar el algoritmo, incluya como tarea aumentar el procedimiento de pivoteo por filas.
# Método de Gauss-Seidel# solución de sistemas de ecuaciones# por métodos iterativosimport numpy as np
# INGRESO
A = np.array([[3. , -0.1, -0.2],
[0.1, 7 , -0.3],
[0.3, -0.2, 10 ]])
B = np.array([7.85,-19.3,71.4])
X0 = np.array([0.,0.,0.])
tolera = 0.00001
iteramax = 100
# PROCEDIMIENTO# Gauss-Seidel
tamano = np.shape(A)
n = tamano[0]
m = tamano[1]
# valores iniciales
X = np.copy(X0)
diferencia = np.ones(n, dtype=float)
errado = 2*tolera
itera = 0
whilenot(errado<=tolera or itera>iteramax):
# por filafor i inrange(0,n,1):
# por columna
suma = 0
for j inrange(0,m,1):
# excepto diagonal de Aif (i!=j):
suma = suma-A[i,j]*X[j]
nuevo = (B[i]+suma)/A[i,i]
diferencia[i] = np.abs(nuevo-X[i])
X[i] = nuevo
errado = np.max(diferencia)
itera = itera + 1
# Respuesta X en columna
X = np.transpose([X])
# revisa si NO convergeif (itera>iteramax):
X=0
# revisa respuesta
verifica = np.dot(A,X)
# SALIDAprint('respuesta X: ')
print(X)
print('verificar A.X=B: ')
print(verifica)
La analogía presentadas entre la «norma como distancia 3D» y el «error de acoplamiento de aeronaves», pertimite considerar desde un punto de partida o inicial las aproximaciones o iteraciones sucesivas hacia una solución del sistema de ecuaciones. Las iteraciones pueden ser convergentes o no.
Los métodos iterativos para sistemas de ecuaciones, son semejantes al método de punto fijo para búsqueda de raíces, requieren un punto inicial para la búsqueda de la raiz o solución que satisface el sistema.
Para describir el método iterativo de Gauss-Seidel, se usa un sistema de 3 incógnitas y 3 ecuaciones, diagonalmente dominante.
Para facilitar la escritura del agoritmo, note el uso de índices ajustado a la descripción de arreglos en Python para la primera fila i=0 y primera columna j=0.
Semejante a despejar una variable de la ecuación para representar un plano, se plantea despejar una variable de cada ecuación. Se obtiene así los valores de cada xi, por cada por cada fila i:
La parte de la sumatoria se realiza para cada término de A[i,j] en la fila i, excepto para el término de la diagonal A[i,i].
Si se tiene conocimiento del problema planteado y se puede «intuir o suponer» una solución para el vector X. Por ejemplo, iniciando con el vector cero, es posible calcular un nuevo vector X usando las ecuaciones para cada X[i] encontradas.
X = \begin{pmatrix} 0 \\ 0 \\ 0 \end{pmatrix}
Con cada nuevo valor se calcula el vector diferencias entre el vector X y cada nuevo valor calculado X[i] .
El error que llama la atención es al mayor valor de las diferencias; se toma como condición para repetir la evaluación de cada vector.
Se observa los resultados de errado para cada iteración, relacionados con la convergencia. Si luego de «muchas» iteraciones se encuentra que (errado>tolera), se detiene el proceso.
El número de condición de una matriz se usa para cuantificar su nivel de mal condicionamiento.
Sea AX=B un sistema de ecuaciones lineales, entonces:
cond(A) = ||A|| ||A-1||
es el número de condición de la matriz A
Tarea
Usando como base los procedimientos desarrollados en Python, elabore un algoritmo para encontrar el número de condición de una matriz.
«el error relativo de la norma de la solución calculada puede ser tan grande como el error relativo de la norma de los coeficientes de [A], multiplicada por el número de condición.»
Por ejemplo,
si los coeficientes de [A] se encuentran a t dígitos de precisión
(esto es, los errores de redondeo son del orden de 10–t) y
Cond [A] = 10c,
la solución [X] puede ser válida sólo para t – c dígitos
(errores de redondeo ~ 10c–t).
verifique el resultado obtenido con el algoritmo, comparando con usar la instrucción
Ejercicio 1. Usando los conceptos de normas mostradas, para el siguiente vector:
x= [5, -3, 2]
a) calcule las normas mostradas (en papel),
b) Realice los respectivos algoritmos en python,
c) Determine los tiempos de ejecución de cada algoritmo. ¿Cúál es el más rápido?
Ejercicio 2. Usando los conceptos de normas mostradas, para la siguiente matriz:
A = [[5, -3, 2],
[4, 8,-4],
[2, 6, 1]]
Repita los literales del ejercicio anterior.
Nota: para convertir una lista X en arreglo use: np.array(X)
Normas con Python y Numpy
Algunas normas vectoriales y matriciales con Python. Cálculo del número de condición.
Se presenta un ejemplo usando la matriz A y el vector B en un programa de prueba.
A = np.array([[3,-0.1,-0.2],
[0.1,7,-0.3],
[0.3,-0.2,10]])
B = np.array([7.85,-19.3,71.4])
Instrucciones en Python:
# Normas vectoriales y matriciales# Referencia: Chapra 10.3, p299, pdf323import numpy as np
defnorma_p(X,p):
Xp = (np.abs(X))**p
suma = np.sum(Xp)
norma = suma**(1/p)
return(norma)
defnorma_euclidiana(X):
norma = norma_p(X,2)
return(norma)
defnorma_filasum(X):
sfila = np.sum(np.abs(X),axis=1)
norma = np.max(sfila)
return(norma)
defnorma_frobenius(X):
tamano = np.shape(X)
n = tamano[0]
m = tamano[1]
norma = 0
for i inrange(0,n,1):
for j inrange(0,m,1):
norma = norma + np.abs(X[i,j])**2
norma = np.sqrt(norma)
return(norma)
defnum_condicion(X):
M = np.copy(X)
Mi = np.linalg.inv(M)
nM = norma_filasum(M)
nMi= norma_filasum(Mi)
ncondicion = nM*nMi
return(ncondicion)
# Programa de prueba ######## INGRESO
A = np.array([[3,-0.1,-0.2],
[0.1,7,-0.3],
[0.3,-0.2,10]])
B = np.array([7.85,-19.3,71.4])
p = 2
# PROCEDIMIENTO
normap = norma_p(B, p)
normaeucl = norma_euclidiana(B)
normafilasuma = norma_filasum(A)
numerocondicion = num_condicion(A)
# SALIDAprint('vector:',B)
print('norma p: ',2)
print(normap)
print('norma euclididana: ')
print(normaeucl)
print('******')
print('matriz: ')
print(A)
print('norma suma fila: ',normafilasuma)
print('número de condición:')
print(numerocondicion)
La norma de un vector se interpreta como una distancia entre la coordenada definida por el vector [xi, yi, zi] y el origen [0,0,0]. También se puede realizar respecto a otro punto de referencia, se conoce como Norma Ecuclidiana o Norma p=2.
Previamente, en busqueda de raíces, el error de aproximación se considera como la diferencia entre dos iteraciones sucesivas:.
error = |x_{i+1} - x_{i}|
Si el concepto se extiende a vectores en tres dimensiones, se observa como el error entre vectores |Xi+1 – Xi| que se interpreta mejor como una distancia.
Por ejemplo, si se usan los puntos y se visualizan en un gráfico:
X1 = [ 1 2 3]
X2 = [ 2 4 -1]
errado = |[ 1 2 -4]|
La diferencia entre los vectores |[ 1, 2, -4]| es más simple de observar como un número escalar. Una forma de convertir el punto a un escalar es usar la distancia al origen.
errado = \sqrt{1^2 + 2^2 + (-4)^2} = 4.58
El resultado también se interpreta como la distancia entre los puntos X1 y X2
De ampliar el concepto a n dimensiones se conoce como norma de un vector o matriz.
¿Considera importante que exista acoplamiento para iniciar la carga el combustible? o ¿para abrir la escotilla del transbordador?
Si el error de acoplamiento entre artefactos se calcula como la Norma entre los puntos de «contacto»,
¿es necesario calcular la raiz cuadrada de los cuadrados de las diferencias? o,
¿Solo toma en cuenta el mínimo o el máximo de las diferencias entre las coordenadas?,
¿cuál de las formas tiene menos operaciones, es más rápida de realizar?
considere sus respuestas para el siguiente concepto.
Norma infinito
Se determina como el valor mas grande entre los elementos del vector o matriz.
||x|| = max\Big[ |x_i| \Big]
Es más sencilla de calcular que la Norma Ecuclidiana, Norma P=2, mostrada al principio.
Se interpreta como si alguna de las diferencia entre las coordenadas de los puntos de acople entre aeronaves es mayor que lo tolerado, no se debe permitir abrir la válvula de combustible o la escotilla del transbordador. No es prioritario calcular la suma de los cuadrados de las diferencias para saber que no existe acople aún.
Existen otras formas de calcular la Norma, que se presentan en la siguiente página web.
Algoritmo en Python
Principalmente se usa para generar las gráficas 3D, que ayudan a la interpretación del concepto con los puntos de coordenadas:
El método de Gauss-Jordan es semejante al método de Gauss en los procedimientos para obtener:
la matriz aumentada,
pivoteada por filas
eliminación hacia adelante
El cambio se presenta a partir de la matriz triangular superior, donde se aplica el procedimiento para obtener la solución:
eliminación hacia atrás
Método de Gauss – Jordan
El método de Gauss-Jordan presenta un procedimiento alterno al de «sustitución hacia atrás» realizado para el método de Gauss.
A partir de haber terminado el procedimiento de «eliminación hacia adelante» y haber obtenido la «matriz triangular superior» aumentada, se aplica el procedimiento de. «eliminación hacia atrás».
Se continúa con el ejercicio desde la «matriz triangular superior» aumentada:
La solución del sistema de ecuaciones se presenta como una matriz identidad concatenada a un vector columa de constantes.
solución X es:
[[2.8]
[4.5]
[8.1]]
X= \begin{pmatrix} 2.8\\ 4.5 \\ 8.1 \end{pmatrix}
Observación: en la matriz hay unos valores del orden de 10-16, que corresponden a errores de operaciones en computadora (truncamiento y redondeo) que pueden ser descartados por ser casi cero. Hay que establecer entonces un parámetro para controlar los casos en que la diferencia entre los ordenes de magnitud son por ejemplo menores a 15 ordenes de magnitud 10-15. e implementarlo en los algoritmos.
Algoritmo en Python
Esta sección reutiliza el algoritmo desarrollado para el Método de Gauss, por lo que los bloques de procedimiento son semejantes hasta #eliminación hacia adelante». Se añade el procedimiento de eliminación hacia atras para completar la solución al sistema de ecuaciones.
El algoritmo desarrollado por partes:
# Método de Gauss-Jordan# Solución a Sistemas de Ecuaciones# de la forma A.X=Bimport numpy as np
# INGRESO
A = np.array([[4,2,5],
[2,5,8],
[5,4,3]])
B = np.array([[60.70],
[92.90],
[56.30]])
# PROCEDIMIENTO
casicero = 1e-15 # Considerar como 0# Evitar truncamiento en operaciones
A = np.array(A,dtype=float)
# Matriz aumentada
AB = np.concatenate((A,B),axis=1)
AB0 = np.copy(AB)
# Pivoteo parcial por filas
tamano = np.shape(AB)
n = tamano[0]
m = tamano[1]
# Para cada fila en ABfor i inrange(0,n-1,1):
# columna desde diagonal i en adelante
columna = abs(AB[i:,i])
dondemax = np.argmax(columna)
# dondemax no está en diagonalif (dondemax !=0):
# intercambia filas
temporal = np.copy(AB[i,:])
AB[i,:] = AB[dondemax+i,:]
AB[dondemax+i,:] = temporal
AB1 = np.copy(AB)
# eliminacion hacia adelantefor i inrange(0,n-1,1):
pivote = AB[i,i]
adelante = i + 1
for k inrange(adelante,n,1):
factor = AB[k,i]/pivote
AB[k,:] = AB[k,:] - AB[i,:]*factor
AB2 = np.copy(AB)
# elimina hacia atras
ultfila = n-1
ultcolumna = m-1
for i inrange(ultfila,0-1,-1):
pivote = AB[i,i]
atras = i-1
for k inrange(atras,0-1,-1):
factor = AB[k,i]/pivote
AB[k,:] = AB[k,:] - AB[i,:]*factor
# diagonal a unos
AB[i,:] = AB[i,:]/AB[i,i]
X = np.copy(AB[:,ultcolumna])
X = np.transpose([X])
# SALIDAprint('Matriz aumentada:')
print(AB0)
print('Pivoteo parcial por filas')
print(AB1)
print('eliminacion hacia adelante')
print(AB2)
print('eliminación hacia atrás')
print(AB)
print('solución de X: ')
print(X)
Tarea: implementar caso cuando aparecen ceros en la diagonal para dar respuesta, convertir a funciones cada parte
Una matriz A puede separarse en dos matrices triangulares:
L de tipo triangular inferior
U de tipo triangular superior
que entre ellas tienen la propiedad que: A = L.U
La matriz U se obtiene después de aplicar el proceso de «eliminación hacia adelante» del método de Gauss.
La matriz L contiene los factores usados en el proceso de «eliminación hacia adelante» del método de Gauss, escritos sobre una matriz identidad en las posiciones donde se calcularon.
Ejercicio
Ejemplo Chapra 10.1 p285, pdf309
Presente las matrices LU de la matriz A siguiente:
Realizado a partir del algoritmo de la sección «método de Gauss» y modificando las partes necesarias para el algoritmo.
Para éste algoritmo, se procede desde el bloque de «pivoteo por filas, continuando con el algoritmo de «eliminación hacia adelante» del método de Gauss. Procedimientos que dan como resultado la matriz U.
La matriz L requiere iniciar con una matriz identidad, y el procedimiento requiere que al ejecutar «eliminación hacia adelante» se almacene cada factor con el que se multiplica la fila para hacer cero. El factor se lo almacena en la matriz L, en la posición de dónde se determinó el factor.
# Matrices L y U# Modificando el método de Gaussimport numpy as np
# INGRESO
A = np.array([[ 3. , -0.1, -0.2],
[ 0.1, 7. , -0.3],
[ 0.3, -0.2, 10. ]], dtype=float)
B = np.array([7.85,-19.3,71.4], dtype=float)
# PROCEDIMIENTO
B = np.transpose([B])
AB = np.concatenate((A,B),axis=1)
AB = np.copy(AB)
# Pivoteo parcial por filas
tamano = np.shape(AB)
n = tamano[0]
m = tamano[1]
# Para cada fila en ABfor i inrange(0,n-1,1):
# columna desde diagonal i en adelante
columna = abs(AB[i:,i])
dondemax = np.argmax(columna)
# dondemax no está en diagonalif (dondemax !=0):
# intercambia filas
temporal = np.copy(AB[i,:])
AB[i,:] = AB[dondemax+i,:]
AB[dondemax+i,:] = temporal
AB1 = np.copy(AB)
A1 = np.copy(AB[:,:m-1])
B1 = np.copy(AB[:,m-1])
# eliminacion hacia adelante# se inicializa L
L = np.identity(n,dtype=float)
for i inrange(0,n-1,1):
pivote = AB[i,i]
adelante = i+1
for k inrange(adelante,n,1):
factor = AB[k,i]/pivote
AB[k,:] = AB[k,:] - AB[i,:]*factor
L[k,i] = factor
U = np.copy(AB[:,:m-1])
# SALIDAprint('Pivoteo parcial por filas')
print(AB1)
print('eliminación hacia adelante')
print('Matriz U: ')
print(U)
print('matriz L: ')
print(L)
Si se requiere una respuesta unificada en una variable, se puede convertir en una arreglo de matrices [L,U].
Las matrices L y U se pueden leer como L=LU[0] y U=LU[1]
LU = np.array([L,U])
# SALIDAprint('triangular inferior L')
print(LU[0])
print('triangular superior U')
print(LU[1])
Tarea
Verificar los resultados, y considerar las divisiones para cero y «casicero».