4.4 Videos – acople de aeronaves

Para tener un contexto sobre el cual revisar el tema se tiene los siguientes videos:

1. Acoplamiento de aviones para recarca de combustible . www.AiirSource.com. 30/diciembre/2015.
KC-135 Stratotanker in Action – Aircraft Air Refueling

2. Acoplamiento con estación espacial internacional ISS. RT en español . 2/Julio/2010. 
El carguero ruso Progress M-06M pasó de largo la Estación Espacial Internacional fracasado en su intento de acoplarse

4.1 Normas vs distancias en 3D

Referencia: Burden Cap7.1 9Ed. p432

Normas de un vector en 3D

La norma de un vector se interpreta como una distancia entre la coordenada definida por el vector [xi, yi, zi] y el origen; también se puede realizar respecto a otra referencia.

Para el caso de soluciones matriciales, lo que se usa en un método iterativo es la diferencia entre una solucion xi vs xi+1 que se usa como criterio del error.

Si se observa como el error entre vectores x1 y xi+1se tendría:

x1     =  [ 1  2  3]
x2     =  [ 2  4 -1]
errado =  [-1 -2  4]

se genera con simples instrucciones de python:

# error como la distancia entre dos puntos
# para el caso en 3D
import numpy as np

x0 = np.array([0,0, 0])
x1 = np.array([1,2, 3])
x2 = np.array([2,4,-1])

errado = x1-x2

print('x1 = ', x1)
print('x2 = ', x2)
print('errado = ', errado)

Sin embargo, para visualizar el concepto se obtiene una gráfica:

La norma de x1 se interpreta como la distancia al origen, es decir:

||x|| = \Big[ \sum_{i=0}^{n} x_i ^2 \Big] ^{1/2} ||x_1|| = \sqrt{x^2+y^2+z^2}
errado =  [-1 -2  4]
||errado|| =  4.58257569495584


Observar en el contexto de los temas en video:

La final, si debe existir acoplamiento, es necesario calcular la raiz cuadrada de los cuadrados?

¿ o solo toman en cuenta las diferencias entre las coordenadas?

¿cuál de las formas tiene menos operaciones?


Norma infinito

||x|| = max\Big[ |x_i| \Big]

Extender el proceso a matrices


# error como la distancia entre dos puntos
# para el caso en 3D
import numpy as np

x0 = np.array([0,0, 0])
x1 = np.array([1,2, 3])
x2 = np.array([2,4,-1])

errado = x1-x2

erradonorma = np.linalg.norm(errado)

print('x1 = ', x1)
print('x2 = ', x2)
print('errado = ', errado)
print('||errado|| = ', erradonorma)

# Grafica
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
figura = plt.figure()
ejes = figura.add_subplot(111,projection = '3d')

punto = x0
[x, y , z] = punto
ejes.scatter(x,y,z, c = 'blue',
             marker='o', label = 'x0')

punto = x1
[x, y , z] = punto
ejes.scatter(x,y,z, c = 'red',
             marker='o', label = 'x1')

punto = x2
[x, y , z] = punto
ejes.scatter(x,y,z, c = 'green',
             marker='o', label = 'x2')

linea = np.concatenate(([x0],[x1]),axis = 0)
x = linea[:,0]
y = linea[:,1]
z = linea[:,2]
ejes.plot(x,y,z, label = '||x1||')

linea = np.concatenate(([x0],[x2]),axis = 0)
x = linea[:,0]
y = linea[:,1]
z = linea[:,2]
ejes.plot(x,y,z, label = '||x2||')

linea = np.concatenate(([x1],[x2]),axis = 0)
x = linea[:,0]
y = linea[:,1]
z = linea[:,2]
ejes.plot(x,y,z, label = '||error||')

ejes.set_xlabel('eje x')
ejes.set_ylabel('eje y')
ejes.set_zlabel('eje z')
ejes.legend()
plt.show()

4.1.2 Normas funciones

Referencia: Chapra 10.3, p299, pdf323

Algunas normas vectoriales y matriciales. Cálculo del número de condición.

# Normas vectoriales y matriciales
# Referencia: Chapra 10.3, p299, pdf323
import numpy as np

def norma_p(X,p):
    Xp = (np.abs(X))**p
    suma = np.sum(Xp)
    norma = suma**(1/p)
    return(norma)

def norma_euclidiana(X):
    norma = norma_p(X,2)
    return(norma)

def norma_filasum(X):
    sfila = np.sum(np.abs(X),axis=1)
    norma = np.max(sfila)
    return(norma)

def norma_frobenius(X):
    tamano = np.shape(X)
    n = tamano[0]
    m = tamano[1]
    norma = 0
    for i in range(0,n,1):
        for j in range(0,m,1):
            norma =  norma + np.abs(X[i,j])**2
    norma = np.sqrt(norma)
    return(norma)

def num_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)

# SALIDA
print('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)

cuyos resultados del ejercicio serán:

vector: [  7.85 -19.3   71.4 ]
norma p:  2
74.3779033047
norma euclididana: 
74.3779033047
******
matriz: 
[[  3.   -0.1  -0.2]
 [  0.1   7.   -0.3]
 [  0.3  -0.2  10. ]]
norma suma fila:  10.5
número de condición:
3.61442432483
>>> 

compare sus resultados con las funciones numpy:

np.linalg.norm(A)
np.linalg.cond(A)

http://www.numpy.org/devdocs/reference/generated/numpy.linalg.norm.html

http://www.numpy.org/devdocs/reference/generated/numpy.linalg.cond.html

4.2 Gauss-Seidel – Algoritmo

El método se describe con un sistema de 3 incógnitas y 3 ecuaciones.
Note el uso de índices ajustado a la descripción de python para la primera fila=0 y primera columna=0:

a00x0 + a01x1+ a02x2 = b0
a10x0 + a11x1+ a12x2 = b1
a20x0 + a21x1+ a22x2 = b2

La forma matricial del sistema original de ecuaciones es: [A]{X} = {B}

Se usa los indices i para fila y j para columna, asi un elemento de la matriz es aij .
Para obtener los valores de cada xi, se despeja uno por cada por cada fila i:

x0 = (b0        - a01x1 - a02x2) /a00
x1 = (b1 - a10x0        - a12x2) /a11
x2 = (b2 - a20x0 - a21x1         ) /a22

Observe el patrón de cada fórmula usada para determinar x[i]:

x_i = \bigg(b_i - \sum_{j=0, j\neq i}^n A_{i,j}X_j\bigg) \frac{1}{A_{ii}}
la sumatoria 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 = [0, 0, 0]

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, y se toma como condición para repetir la evaluación de cada vector.

     nuevo = [ 0, 0,  0]
         X = [x0, x1, x2]
diferencia = [|0 - x0|, |0 - x1|, |0 - x2|]
errado = max(diferencia)

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.

Con lo descrito, se muestra una forma de implementar el algoritmo como una función y un ejemplo de prueba con datos del problema conocido.

# Algoritmo Gauss-Seidel,
# matrices, métodos iterativos
# Referencia: Chapra 11.2, p.310, pdf.334
#      Rodriguez 5.2 p.162
# ingresar iteramax si requiere más iteraciones

import numpy as np

def gauss_seidel(A,B,tolera,X,iteramax=100):
    tamano = np.shape(A)
    n = tamano[0]
    m = tamano[1]
    diferencia = np.ones(n, dtype=float)
    errado = np.max(diferencia)
    
    itera = 0
    while not(errado<=tolera or itera>iteramax):
        for i in range(0,n,1):
            nuevo = B[i]
            for j in range(0,m,1):
                if (i!=j): # excepto diagonal de A
                    nuevo = nuevo-A[i,j]*X[j]
            nuevo = nuevo/A[i,i]
            diferencia[i] = np.abs(nuevo-X[i])
            X[i] = nuevo
        errado = np.max(diferencia)
        itera = itera + 1
    # Vector en columna
    X = np.transpose([X])
    # No converge
    if (itera>iteramax):
        X=0
    return(X)

# 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])

tolera = 0.00001

# PROCEDIMIENTO
n = len(B)
X = np.zeros(n, dtype=float)
respuesta = gauss_seidel(A,B,tolera,X)
verifica = np.dot(A,respuesta)

# SALIDA
print('respuesta de A.X=B : ')
print(respuesta)
print('verificar A.X: ')
print(verifica)

que da como resultado:

respuesta de A.X=B : 
[[ 3. ]
 [-2.5]
 [ 7. ]]
verificar A.X: 
[[  7.84999999]
 [-19.3       ]
 [ 71.4       ]]
>>> 

que es la respuesta del problema obtenida durante el desarrollo del ejemplo con otros métodos.

Número de Condición

Referencia: Chapra 10.3.2 p300/pdf324

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

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).

http://blog.espol.edu.ec/matg1013/normas-funciones/

4.1.1 Normas de Vector o Matriz

Referencia: Burden Cap7.1 9Ed. p432

Es una manera de expresar la magnitud de sus componentes:

Sea X un vector de n componentes:

||X|| = \sum_{i=1}^{n}|X_i| ||X|| = max|X_i| , i=1,2, ...,n ||X|| = \left( \sum_{i=1}^{n}X_i^2 \right) ^{1/2}

Sea una matriz A de nxn componentes:

||A|| = max(\sum_{j=1}^{n}|a_{i,j}|, i = 1,2,...,n) ||A|| = max(\sum_{i=1}^{n}|a_{i,j}|, j = 1,2,...,n) ||A|| = \left( \sum_{i=1}^{n}\sum_{j=1}^{n} a_{i,j}^2 \right) ^{1/2}

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)