4.3 Interpolación polinómica de Lagrange con Python

Referencia: Chapra 18.2 pdf 540, Burden 9Ed 3.1 p108, Rodriguez 6.2 p195

El polinomio de interpolación de Lagrange es una reformulación del polinomio de interpolación de Newton que el método evita el cálculo de las diferencias divididas. El método tolera las diferencias entre distancias x de los puntos de muestra.

El polinomio de Lagrange se construye a partir de las fórmulas:

f_{n} (x) = \sum_{i=0}^{n} L_{i} (x) f(x_{i}) L_{i} (x) = \prod_{j=0, j \neq i}^{n} \frac{x-x_j}{x_i - x_j}

Donde una vez que se han seleccionado los puntos a usar que generan la misma cantidad de términos que puntos.

Ejemplo: Polinomio Interpolante de Lagrange con grado 3

Dados los 4 puntos en la tabla se requiere generar un polinomio de grado 3 de la forma:

p(x) = a_0 x^3 + a_1 x^2 + a_2 x^1 + a_3
xi 0 0.2 0.3 0.4
fi 1 1.6 1.7 2.0

Revisando la aplicación de la fórmula se tiene que calcular cuatro términos,

término 1

L_{0} (x) = \frac{(x-0.2)(x-0.3)(x-0.4)}{(0-0.2)(0-0.3)(0-0.4)}

término 2

L_{1} (x) = \frac{(x-0)(x-0.3)(x-0.4)}{(0.2-0)(0.2-0.3)(0.2-0.4)}

término 3

L_{2} (x) = \frac{(x-0)(x-0.2)(x-0.4)}{(0.3-0)(0.3-0.2)(0.3-0.4)}

término 4

L_{3} (x) = \frac{(x-0)(x-0.2)(x-0.3)}{(0.4-0)(0.4-0.2)(0.4-0.3)}

se construye el polinomio usando la fórmula para fn(x) para cada valor fi,

p_3(x) = 1 L_{0} (x) + 1.6 L_{1} (x) + 1.7 L_{2} (x) + 2 L_{3} (x) p_3(x) = 1 \frac{(x-0.2)(x-0.3)(x-0.4)}{(0-0.2)(0-0.3)(0-0.4)} + + 1.6 \frac{(x-0)(x-0.3)(x-0.4)}{(0.2-0)(0.2-0.3)(0.2-0.4)} + 1.7 \frac{(x-0)(x-0.2)(x-0.4)}{(0.3-0)(0.3-0.2)(0.3-0.4)} + 2 \frac{(x-0)(x-0.2)(x-0.3)}{(0.4-0)(0.4-0.2)(0.4-0.3)}

La simplificación de la expresión del polinomio se puede dejar como tarea.

Una forma de verificar que el polinomio es correcto es usar un punto original xi[i] y comprobar que la evaluación tiene como resultado fi[i].

xi = np.array([0, 0.2, 0.3, 0.4])
fi = np.array([1, 1.6, 1.7, 2.0])

Algoritmo en Python

En el algoritmo en Python, para construir las expresiones de cada término se usa la forma simbólica con Sympy.

En cada término Li(x) se usan todos los elementos i , excepto el mismo elemento i, en el numerador y denominador de la expresión.

En polinomio se agrupan todos los términos multiplicados por su respectivo valor fi[i].

    valores de fi:  [1.  1.6 1.7 2. ]
divisores en L(i):  [-0.024  0.004 -0.003  0.008]

Polinomio de Lagrange, expresiones
400.0*x*(x - 0.4)*(x - 0.3) 
  - 566.666666666667*x*(x - 0.4)*(x - 0.2) 
  + 250.0*x*(x - 0.3)*(x - 0.2) 
  - 41.6666666666667*(x - 0.4)*(x - 0.3)*(x - 0.2)

Polinomio de Lagrange: 
41.666666666667*x**3 - 27.5*x**2 + 6.8333333333334*x + 1.0

Las expresiones del polinomio contiene los binomios en su forma básica, para resolver y simplificar las ecuaciones se usa polinomio.expand().

Para realizar la gráfica del polinomio es conveniente convertirlo a la forma lambda con Numpy, de esta forma se evalúa en una sola linea todos los puntos para el intervalo [a,b] en x.

# Interpolacion de Lagrange
# divisoresL solo para mostrar valores
import numpy as np
import sympy as sym
import matplotlib.pyplot as plt

# INGRESO , Datos de prueba
xi = np.array([0, 0.2, 0.3, 0.4])
fi = np.array([1, 1.6, 1.7, 2.0])

# PROCEDIMIENTO
# Polinomio de Lagrange
n = len(xi)
x = sym.Symbol('x')
polinomio = 0
divisorL = np.zeros(n, dtype = float)
for i in range(0,n,1):
    
    # Termino de Lagrange
    numerador = 1
    denominador = 1
    for j  in range(0,n,1):
        if (j!=i):
            numerador = numerador*(x-xi[j])
            denominador = denominador*(xi[i]-xi[j])
    terminoLi = numerador/denominador

    polinomio = polinomio + terminoLi*fi[i]
    divisorL[i] = denominador

# simplifica el polinomio
polisimple = polinomio.expand()

# para evaluación numérica
px = sym.lambdify(x,polisimple)

# Puntos para la gráfica
muestras = 101
a = np.min(xi)
b = np.max(xi)
pxi = np.linspace(a,b,muestras)
pfi = px(pxi)

# SALIDA
print('    valores de fi: ',fi)
print('divisores en L(i): ',divisorL)
print()
print('Polinomio de Lagrange, expresiones')
print(polinomio)
print()
print('Polinomio de Lagrange: ')
print(polisimple)

# Gráfica
plt.plot(xi,fi,'o', label = 'Puntos')
plt.plot(pxi,pfi, label = 'Polinomio')
plt.legend()
plt.xlabel('xi')
plt.ylabel('fi')
plt.title('Interpolación Lagrange')
plt.show()

4.2.2 Diferencias divididas de Newton

Referencia: Chapra 18.1.3 p508 pdf532, Rodriguez 6.7 p223, Burden 9Ed 3.3 p124.

El método se usa en el caso que los puntos en el «eje x» se encuentran espaciados de forma arbitraria y provienen de una función desconocida pero supuestamente diferenciable.

La n-ésima diferencia dividida finita es:

f[x_{n}, x_{n-1}, \text{...}, x_{1}, x_{0}] = \frac{f[x_{n}, x_{n-1}, \text{...}, x_{1}]- f[x_{n-1}, x_{n-2}, \text{...}, x_{0}]}{x_{n}-x_{0}}

Para lo cual se debe interpretar la tabla de diferencias divididas, también como una aproximación a una derivada:

i xi f[xi] Primero Segundo Tercero
0 x0 f[x0] f[x1,x0] f[x2,x1,x0] f[x3,x2,x1,x0]
1 x1 f[x1] f[x2,x1] f[x3,x2,x1]
2 x2 f[x2] f[x3,x2]
3 x2 f[x3]

En la fórmula del polinomio, las diferencias divididas sirven para evaluar los coeficientes de cada término adicional para aumentar el grado. Semejante al proceso realizado para  «diferencias finitas divididas»:

f_n(x) = f(x_0)+(x-x_0)f[x_1,x_0] + + (x-x_0)(x-x_1)f[x_2,x_1,x_0] + \text{...}+ + (x-x_0)(x-x_1)\text{...}(x-x_{n-1})f[x_n,x_{n-1},\text{...},x_0]

Este polinomio de interpolación se conoce como de Newton en diferencias divididas.


Ejercicio

El ejercicio perminte mostrar con mayor detalle las operaciones del procedimiento. Para el ejercicio, se toma como referencia el descrito en el enlace, observando que las distancias entre puntos para el «eje x» NO son equidistantes,

Referencia: 1Eva_IIT2008_T3_MN Ganancia en inversión

Los datos del problema muestran varios valores de inversión xi que generan cada uno ganancias fi.

Considere que los valores invertidos en materia prima para producción, dependen de la demanda de del producto en el mercado, motivo por el que los valores de inversión no guardan distancias equidistante entre si.

Siguiento del formato usado en los algoritmos, los datos son:

xi = np.array([3.2 , 3.8 , 4.2 , 4.5 ])
fi = np.array([5.12, 6.42, 7.25, 6.85])

Con los datos se llena la tabla de  diferencias divididas, donde por simplicidad, se escriben las operaciones en cada casilla. La última de deja como tarea.:

i xi f[xi] Primero Segundo Tercero
0 3.2 5.12 \frac{6.42-5.12}{3.8-3.2} =2.1667 \frac{2.075-2.1667}{4.2-3.2} =-0.0917 f[x3,x2,x1,x0]
1 3.8 6.42 \frac{7.25-6.42}{4.2-3.8} =2.075 \frac{-1.3333-2.075}{4.5-3.8} =-4.869
2 4.2 7.25 \frac{6.85-7.25}{4.5-4.2} =-1.3333
3 4.5 6.85

Las diferencias divididas de la primera fila son los valores usados para la expresión del polinomio de interpolacion:

ddividida = tabla[0,3:] = [ 2.1667, -0.0917, -3.6749, 0. ]

La expresión del polinomio inicia con el primer valor de la función f(x0).

p_3(x) = f_0 = 5.12

se calcula el primer término usando el factor de la primera diferencia dividida y se añade a la expresión del polinomio.

término = (x-x_0)f[x1,x0] = (x-3.2)2.1667 p_3(x) = 5.12 + 2.1667(x-3.2)

con el siguiente valor de ddividida[] se procede de manera semejante:

término = (x-x_0)(x-x_1)f[x1,x0] = = (x-3.2)(x-3.8)(-0.0917) p_3(x) = 5.12 + 2.1667*(x-3.2)+ + (-0.0917)(x-3.2)(x-3.8)

Realizando el cálculo del tercer y último termino con diferencias divididas,  el polinomio obtenido es:

p_3(x) = 5.12 + 2.1667(x-3.2) + + (-0.0917)(x-3.2)(x-3.8) + + (-3.6749)(x - 3.2)(x - 3.8)(x - 4.2)

que simplificando al multiplicar entre los términos (x-xi):

p_3(x) = 184.7569 - 149.9208 x + 41.0673 x^2 - 3.6749 x^3

El polinomio graficado en el intervalo de xi muestra que pasa por todos los puntos.


Algoritmo en Python

El algoritmo para  interpolación de diferencias divididas o Newton, considera reutilizar el procedimiento de cálculo de diferencias finitas incorporando la parte del denominador.

La creación de la expresión del polinomio también es semejante a la usada para diferencias finitas avanzadas.

Se incorpora la parte gráfica para observar los resultados en el intervalo xi, con el número de muestras = 101 para tener una buena resolución de la línea del polinomo.

# Polinomio interpolación
# Diferencias Divididas de Newton
# Tarea: Verificar tamaño de vectores,
#        verificar puntos equidistantes en x
import numpy as np
import sympy as sym
import matplotlib.pyplot as plt

# INGRESO , Datos de prueba
xi = np.array([3.2, 3.8, 4.2, 4.5])
fi = np.array([5.12, 6.42, 7.25, 6.85])

# PROCEDIMIENTO

# Tabla de Diferencias Divididas Avanzadas
titulo = ['i   ','xi  ','fi  ']
n = len(xi)
ki = np.arange(0,n,1)
tabla = np.concatenate(([ki],[xi],[fi]),axis=0)
tabla = np.transpose(tabla)

# diferencias divididas vacia
dfinita = np.zeros(shape=(n,n),dtype=float)
tabla = np.concatenate((tabla,dfinita), axis=1)

# Calcula tabla, inicia en columna 3
[n,m] = np.shape(tabla)
diagonal = n-1
j = 3
while (j < m):
    # Añade título para cada columna
    titulo.append('F['+str(j-2)+']')

    # cada fila de columna
    i = 0
    paso = j-2 # inicia en 1
    while (i < diagonal):
        denominador = (xi[i+paso]-xi[i])
        numerador = tabla[i+1,j-1]-tabla[i,j-1]
        tabla[i,j] = numerador/denominador
        i = i+1
    diagonal = diagonal - 1
    j = j+1

# POLINOMIO con diferencias Divididas
# caso: puntos equidistantes en eje x
dDividida = tabla[0,3:]
n = len(dfinita)

# expresión del polinomio con Sympy
x = sym.Symbol('x')
polinomio = fi[0]
for j in range(1,n,1):
    factor = dDividida[j-1]
    termino = 1
    for k in range(0,j,1):
        termino = termino*(x-xi[k])
    polinomio = polinomio + termino*factor

# simplifica multiplicando entre (x-xi)
polisimple = polinomio.expand()

# polinomio para evaluacion numérica
px = sym.lambdify(x,polisimple)

# Puntos para la gráfica
muestras = 101
a = np.min(xi)
b = np.max(xi)
pxi = np.linspace(a,b,muestras)
pfi = px(pxi)

# SALIDA
np.set_printoptions(precision = 4)
print('Tabla Diferencia Dividida')
print([titulo])
print(tabla)
print('dDividida: ')
print(dDividida)
print('polinomio: ')
print(polinomio)
print('polinomio simplificado: ' )
print(polisimple)

# Gráfica
plt.plot(xi,fi,'o', label = 'Puntos')
##for i in range(0,n,1):
##    plt.axvline(xi[i],ls='--', color='yellow')
plt.plot(pxi,pfi, label = 'Polinomio')
plt.legend()
plt.xlabel('xi')
plt.ylabel('fi')
plt.title('Diferencias Divididas - Newton')
plt.show()

El resultado del algoritmo se muestra a continuación:

Tabla Diferencia Dividida
[['i   ', 'xi  ', 'fi  ', 'F[1]', 'F[2]', 'F[3]', 'F[4]']]
[[ 0.      3.2     5.12    2.1667 -0.0917 -3.6749  0.    ]
 [ 1.      3.8     6.42    2.075  -4.869   0.      0.    ]
 [ 2.      4.2     7.25   -1.3333  0.      0.      0.    ]
 [ 3.      4.5     6.85    0.      0.      0.      0.    ]]
dDividida: 
[ 2.1667 -0.0917 -3.6749  0.    ]
polinomio: 
2.16666666666667*x - 3.67490842490842*(x-4.2)*(x-3.8)*(x-3.2)
 - 0.0916666666666694*(x - 3.8)*(x - 3.2) - 1.81333333333334
polinomio simplificado: 
-3.67490842490842*x**3 + 41.0673076923077*x**2 
- 149.920860805861*x + 184.756923076923
>>> 

4.2.1 Diferencias finitas avanzadas – polinomio de interpolación

Referencia: Rodriguez 6.6.4 Pdf.221, Burden 9Ed p129

Se usa en interpolación cuando los puntos en el «eje x» se encuentran igualmente espaciados, la diferencia entre puntos consecutivos xi es una constante denominada h.

h = xi+1 – xi


Una relación entre derivadas y diferencias finitas se establece mediante:

f^{(n)}(z) = \frac{\Delta ^{n} f_{0}}{h^{n}} para algún k en el intervalo [x0,xn]

\frac{\Delta ^{n} f_{0}}{h^{n}} es una aproximación para la n-ésima derivada f(n) en el intervalo [x0,xn]

El polinomio de interpolación se puede contruir por medio de diferencias finitas avanzadas con las siguiente fórmula:

p_n (x) = f_0 + \frac{\Delta f_0}{h} (x - x_0) + + \frac{\Delta^2 f_0}{2!h^2} (x - x_0)(x - x_1) + + \frac{\Delta^3 f_0}{3!h^3} (x - x_0)(x - x_1)(x - x_2) + \text{...} + \frac{\Delta^n f_0}{n!h^n} (x - x_0)(x - x_1) \text{...} (x - x_{n-1})

Observe que en la medida que se toman más puntos de muestra, el grado del polinomio puede ser mayor.


Ejercicio

Se toman los datos del ejercicio de «diferencias finitas» y se grafican los resultados. Se resalta en la gráfica la distancia equidistante h entre los puntos consecutivos xi como los mostrados en la gráfica anterior.

xi = np.array([0.10, 0.2, 0.3, 0.4])
fi = np.array([1.45, 1.6, 1.7, 2.0])

y se usan los resulados previos de diferencias finitas:

Tabla Diferencia Finita: 
[['i', 'xi', 'fi', 'df1', 'df2', 'df3', 'df4']]
[[ 0.    0.1   1.45  0.15 -0.05  0.25  0.  ]
 [ 1.    0.2   1.6   0.1   0.2   0.    0.  ]
 [ 2.    0.3   1.7   0.3   0.    0.    0.  ]
 [ 3.    0.4   2.    0.    0.    0.    0.  ]]

valores con los que se inicia el cálculo de la distancia constante entre puntos xi:

h = xi[1] – xi[0] = 0.2-0.1 = 0.1

Los valores de diferencias finitas a usar en la fórmula son de la primera fila desde la columna 3 en adelante.

dfinita = tabla[0,3:] = [0.15, -0.05, 0.25, 0.]

La construción de la expresión del polinomio empieza con el valor del primer punto de la función f(0.1)

p_3(x) = f_0 = 1.45

para añadirle el primer término j=1, hay que calcular el factor:

factor = \frac{\Delta^1 f_0}{1!h^1} = \frac{0.15}{0.1} = 1.5

con lo que el término se completa como:

término = factor(x-x_0) = 1.5(x-0.1) p_3(x) = 1.45 + 1.5(x-0.1)

Para el segundo término j=2, se repite el proceso:

factor = \frac{\Delta^2 f_0}{2!h^2} = \frac{-0.05}{2!(0.1)^2} = -2.5 término = factor(x-x_0) = -2.5(x-0.1)(x-02) p_3(x)= 1.45 + 1.5(x-0.1) +(-2.5)(x-0.1)(x-02)

Finalmente, añadiendo el tercer término j=3 cuyo cálculo es semejante a los anteriores, se deja como tarea.

El resultado del método es:

p_3(x)= 1.45 + 1.5(x-0.1) -2.5(x-0.1)(x-02) + + 41.667(x - 0.3)(x - 0.2)(x - 0.1)

Se puede seguir simpificando la respuesta, por ejemplo usando solo el término de grado con 1.5(x-0-1) se tiene que:

p_3(x) = 1.3 + 1.5x -2.5(x-0.1)(x-02)+ + 41.667(x - 0.3)(x - 0.2)(x - 0.1)

Seguir simplificando la expresión en papel, se deja como tarea.

La gráfica del polinomio obtenido comprueba que pasa por cada uno de los puntos dados para el ejercicio.


Algoritmo en Python

El polinomio se contruye usando el ejercicio de diferencias finitas.

Para contruir la expresión del polinomio añadiendo los términos de la fórmula, se define la variable simbólica x con sympy.

Para simplificar el polinomio resultante con las expresiones de multiplicación, se utiliza la instrucción sym.expand().

En caso de requerir evaluar la fórmula con un vector de datos se la convierte a la forma lambda.

Se añaden las instrucciones para realizar la gráfica en el intervalo [a,b] dado por los valores xi, definiendo el número de muestras = 101 que son suficientes para que la gráfica se observe sin distorsión.

# Polinomio interpolación
# Diferencias finitas avanzadas
# Tarea: Verificar tamaño de vectores,
#        verificar puntos equidistantes en x

import numpy as np
import sympy as sym
import matplotlib.pyplot as plt

# INGRESO , Datos de prueba
xi = np.array([0.10, 0.2, 0.3, 0.4])
fi = np.array([1.45, 1.6, 1.7, 2.0])

# PROCEDIMIENTO

# Tabla de Diferencias Finitas
titulo = ['i','xi','fi']
n = len(xi)
ki = np.arange(0,n,1)
tabla = np.concatenate(([ki],[xi],[fi]),axis=0)
tabla = np.transpose(tabla)

# diferencias finitas vacia
dfinita = np.zeros(shape=(n,n),dtype=float)
tabla = np.concatenate((tabla,dfinita), axis=1)

# Calcula tabla, inicia en columna 3
[n,m] = np.shape(tabla)
diagonal = n-1
j = 3
while (j < m):
    # Añade título para cada columna
    titulo.append('df'+str(j-2))
    # cada fila de columna
    i = 0
    while (i < diagonal):
        tabla[i,j] = tabla[i+1,j-1]-tabla[i,j-1]
        i = i+1
    diagonal = diagonal - 1
    j = j+1

# POLINOMIO con diferencias Finitas avanzadas
# caso: puntos equidistantes en eje x
h = xi[1] - xi[0]
dfinita = tabla[0,3:]
n = len(dfinita)

# expresión del polinomio con Sympy
x = sym.Symbol('x')
polinomio = fi[0]
for j in range(1,n,1):
    denominador = np.math.factorial(j)*(h**j)
    factor = dfinita[j-1]/denominador
    termino = 1
    for k in range(0,j,1):
        termino = termino*(x-xi[k])
    polinomio = polinomio + termino*factor

# simplifica multiplicando entre (x-xi)
polisimple = polinomio.expand()

# polinomio para evaluacion numérica
px = sym.lambdify(x,polisimple)

# Puntos para la gráfica
muestras = 101
a = np.min(xi)
b = np.max(xi)
pxi = np.linspace(a,b,muestras)
pfi = px(pxi)

# SALIDA
print('Tabla Diferencia Finita')
print([titulo])
print(tabla)
print('dfinita: ')
print(dfinita)
print('polinomio: ')
print(polinomio)
print('polinomio simplificado: ' )
print(polisimple)

# Gráfica
plt.plot(xi,fi,'o', label = 'Puntos')
##for i in range(0,n,1):
##    plt.axvline(xi[i],ls='--', color='yellow')
plt.plot(pxi,pfi, label = 'Polinomio')

plt.legend()
plt.xlabel('xi')
plt.ylabel('fi')
plt.title('Interpolación polinómica')
plt.show()

teniendo como resultado:

Tabla Diferencia Finita
[['i', 'xi', 'fi', 'df1', 'df2', 'df3', 'df4']]
[[ 0.    0.1   1.45  0.15 -0.05  0.25  0.  ]
 [ 1.    0.2   1.6   0.1   0.2   0.    0.  ]
 [ 2.    0.3   1.7   0.3   0.    0.    0.  ]
 [ 3.    0.4   2.    0.    0.    0.    0.  ]]
dfinita: 
[ 0.15 -0.05  0.25  0.  ]
polinomio: 
1.5*x + 41.6666666666667*(x - 0.3)*(x - 0.2)*(x - 0.1) 
- 2.50000000000001*(x - 0.2)*(x - 0.1) + 1.3
polinomio simplificado: 
41.6666666666667*x**3 - 27.5*x**2 
+ 6.83333333333335*x + 0.999999999999999

el polinomio de puede evaluar como px(valor) una vez que se convirtió a lambda:

>>> px(0.1)
1.4500000000000004
>>> px(0.2)
1.6000000000000025
>>> px0.3)
1.7000000000000042
>>> 

Como tarea se recomienda realizar las gráficas comparativas entre métodos:

4.2 Diferencias finitas con Python

Referencia: Rodriguez 6.5 Pdf.211, Chapra 18.1.3 pdf.509

Establece relaciones simples entre los puntos dados que describen la funcion f(x). Las tabla de diferencias finitas es un elemento base para varios métodos de interpolación, por lo que se trata como un tema inicial.

Para un ejemplo se toman los siguientes puntos:

xi 0.10 0.2 0.3 0.4
fi 1.45 1.6 1.7 2.0

La tabla de diferencias finitas se construye tomando los datos y sus índices como parte de las primeras columnas.

i xi fi Δ1fi Δ2fi Δ3fi Δ4fi
0 0.1 1.45 0.15 -0.05 0.25 0.
1 0.2 1.6 0.1 0.2 0. 0.
2 0.3 1.7 0.3 0. 0. 0.
3 0.4 2.0 0. 0. 0. 0.

Cada casila de diferencia finita Δjfi se obtiene restando los dos valores consecutivos de la columna anterior. Por ejemplo, para la primera columna:

\Delta ^1 f_0 = 1.6-1.45 = 0.15 \Delta ^1 f_1 = 1.7-1.6 = 0.1 \Delta ^1 f_2 = 2.0-1.7 = 0.3

y la siguiente Δ1f3 no es posible calcular.

Siguiendo el mismo procedimiento se calculan los valores de la siguiente columna Δ2fi como las diferencias de la columna anterior, así sucesivamente.

Si la función f(x) de donde provienen los datos es un polinonio de grado n, entonces la n-ésima diferencia finita será una constante, y las siguientes diferencias se anularán.


Algoritmo en Python

Para crear la tabla de diferencias finitas, las primeras columnas requieren concatenar los valores de los índice ixi y fi.

xi = np.array([0.10, 0.2, 0.3, 0.4])
fi = np.array([1.45, 1.6, 1.7, 2.0])

Los índices i se crean en un vector ki, pues la variable i es usada como fila en matrices, por lo que evitamos confundirlas al usar la variable.

A la matriz con las tres columnas descritas, se le añade a la derecha una matriz de nxn para calcular las diferencias.

Se calculan las diferencias para cada columna, realizando la operación entre los valores de cada fila. Considere que no se realizan cálculos desde la diagonal hacia abajo en la tabla, los valores quedan como cero.

Al final se muestra el título y el resultado de la tabla.

# Tabla de Diferencias finitas
# resultado en: [título,tabla]
# Tarea: verificar tamaño de vectores
import numpy as np

# INGRESO, Datos de prueba
xi = np.array([0.10, 0.2, 0.3, 0.4])
fi = np.array([1.45, 1.6, 1.7, 2.0])

# PROCEDIMIENTO

# Tabla de Diferencias Finitas
titulo = ['i','xi','fi']
n  = len(xi)
ki = np.arange(0,n,1)
tabla = np.concatenate(([ki],[xi],[fi]),axis=0)
tabla = np.transpose(tabla)

# diferencias finitas vacia
dfinita = np.zeros(shape=(n,n),dtype=float)
tabla   = np.concatenate((tabla,dfinita), axis=1)

# Calcula tabla, inicia en columna 3
[n,m] = np.shape(tabla)
diagonal = n-1
j = 3
while (j < m):
    # Añade título para cada columna
    titulo.append('df'+str(j-2))

    # cada fila de columna
    i = 0
    while (i < diagonal):
        tabla[i,j] = tabla[i+1,j-1]-tabla[i,j-1]
        i = i + 1

    diagonal = diagonal - 1
    j = j + 1

# SALIDA
print('Tabla Diferencia Finita: ')
print([titulo])
print(tabla)

el resultado de aplicar el algoritmo es:

Tabla Diferencia Finita: 
[['i', 'xi', 'fi', 'df1', 'df2', 'df3', 'df4']]
[[ 0.    0.1   1.45  0.15 -0.05  0.25  0.  ]
 [ 1.    0.2   1.6   0.1   0.2   0.    0.  ]
 [ 2.    0.3   1.7   0.3   0.    0.    0.  ]
 [ 3.    0.4   2.    0.    0.    0.    0.  ]]

>>> 

Gráfica de puntos con líneas horizontales

Para tener una referencia visual sobre las primeras diferencias finitas, en una gráfica se trazan las líneas horizontales que pasan por cada punto. Para las segundas diferencias se debe graficar las primeras diferencias finitas vs xi repitiendo el proceso. Las líneas de distancia se añadieron con un editor de imágenes.

las intrucciones adicionales al algoritmo anterior para añadir la gráfica son:

# Gráfica
import matplotlib.pyplot as plt

for i in range(0,n,1):
    plt.axhline(fi[i],ls='--', color='yellow')

plt.plot(xi,fi,'o', label = 'Puntos')

plt.legend()
plt.xlabel('xi')
plt.ylabel('fi')
plt.title('Diferencia Finita')
plt.show()

4.1 El polinomio de interpolación con Sympy-Python

Referencia: Rodriguez 6.1 p191, Chapra 18.3 p520 pdf544

Unir n puntos en el plano usando un polinomios de interpolación se puede realizar con un polinomio de grado (n-1) y la misma cantidad de puntos en el plano.

Por ejemplo, dados los 4 puntos en la tabla se requiere generar un polinomio de grado 3 de la forma:

p(x) = a_0 x^3 + a_1 x^2 + a_2 x^1 + a_3
xi 0 0.2 0.3 0.4
fi 1 1.6 1.7 2.0

Es posible generar un sistema de ecuaciones en base p(x) que pase por cada uno de los puntos. Por ejemplo si se toma el primer punto con xi=0 y fi=1 se establece la ecuación:

p(0) = 1 a_0 (0)^3 + a_1 (0)^2 + a_2 (0)^1 + a_3 = 1

Note que ahora las incógnitas son los coeficientes ai. Luego se continúa con los otros puntos seleccionados hasta completar las ecuaciones necesarias para el grado del polinomio seleccionado.

a_0 (0.0)^3 + a_1 (0.0)^2 + a_2 (0.0)^1 + a_3 = 1.0 a_0 (0.2)^3 + a_1 (0.2)^2 + a_2 (0.2)^1 + a_3 = 1.6 a_0 (0.3)^3 + a_1 (0.3)^2 + a_2 (0.3)^1 + a_3 = 1.7 a_0 (0.4)^3 + a_1 (0.4)^2 + a_2 (0.4)^1 + a_3 = 2.0

El sistema obtenido se resuelve con alguno de los métodos conocidos.

Los métodos requieren usar la matriz A y vector B:

\begin{pmatrix} 0.0^3 & 0.0^2 & 0.0^1 & 1 \\ 0.2^3 & 0.2^2 & 0.2^1 & 1 \\ 0.3^3 & 0.3^2 & 0.3^1 & 1 \\ 0.4^3 & 0.4^2 & 0.4^1 & 1 \end{pmatrix} . \begin{pmatrix} a_0 \\ a_1 \\ a_2 \\ a_3 \end{pmatrix} = \begin{pmatrix} 1.0 \\ 1.6 \\ 1.7 \\ 2.0 \end{pmatrix}

La matriz A también se conoce como Matriz Vandermonde D, se construye observando que los coeficientes se elevan al exponente referenciado al índice columna pero de derecha a izquierda. La última columa da 1 por tener como coeficiente el valor de xi0.

Por los métodos para resolver sistemas de ecuaciones, se observa que hay que realizar la matriz aumentada, pivotear por filas, etc.

Para la solución se propone un algoritmo realizado en Python, que para enfocarnos en la interpolación, se obtiene el siguente resultado:

los coeficientes del polinomio: 
[[ 41.66666667]
 [-27.5       ]
 [  6.83333333]
 [  1.        ]]

con lo que se contruye el polinomio con los valores ai.

Polinomio de interpolación p(x): 
41.6666666666667*x**3 - 27.5*x**2 + 6.83333333333333*x + 1.0

que re-escrito es

p(x) = 41.6666 x^3 - 27.5 x^2 + 6.8333 x + 1

Polinomio que se grafica dentro del intervalo de los datos del problema y se obtiene la línea que pasa por los puntos.


Algoritmo en Python

Para desarrollar el ejercicio, se realiza un bloque para construir la matriz A y vector B, usando los vectores que representan los puntos de muestra de la función o experimento.

xi = [0,0.2,0.3,0.4]
fi = [1,1.6,1.7,2.0]

Observe que en éste caso los datos se dan en forma de lista y se deben convertir hacia arreglos.

La matriz A o Matriz Vandermonde D, se construye observando que los coeficientes del vector xi se elevan al exponente referenciado al índice columna pero de derecha a izquierda.

Construir el polinomio consiste en resolver el sistema de ecuaciones indicado en la sección anterior. Para simplificar la solución del sistema, se usa numpy para algebra lineal, que entrega el vector solución que representan los coeficientes del polinomio de interpolación.

Para contruir la expresión del polinomio, se usa la forma simbólica con Sympy, de forma semejante a la usada para construir el polinomio de Taylor.

Para facilitar la evación del polinomio de forma numérica, se convierte el polinomio a la forma Lambda. Con un número de muestras suficientes para suavizar la curva dentro del intervalo [a,b] del problema se puede realizar la gráfica del polinomio.

Finalmente, con la gráfica se comprueba que la curva pase por todos los puntos del intervalo.

# El polinomio de interpolación
import numpy as np
import sympy as sym
import matplotlib.pyplot as plt

# INGRESO
xi = [0,0.2,0.3,0.4]
fi = [1,1.6,1.7,2.0]
# muestras = tramos+1
muestras = 101

# PROCEDIMIENTO
# Convierte a arreglos numpy 
xi = np.array(xi)
B = np.array(fi)
n = len(xi)

# Matriz Vandermonde D
D = np.zeros(shape=(n,n),dtype =float)
for i in range(0,n,1):
    for j in range(0,n,1):
        potencia = (n-1)-j # Derecha a izquierda
        D[i,j] = xi[i]**potencia

# Aplicar métodos Unidad03. Tarea
# Resuelve sistema de ecuaciones A.X=B
coeficiente = np.linalg.solve(D,B)

# Polinomio en forma simbólica
x = sym.Symbol('x')
polinomio = 0
for i in range(0,n,1):
    potencia = (n-1)-i   # Derecha a izquierda
    termino = coeficiente[i]*(x**potencia)
    polinomio = polinomio + termino

# Polinomio a forma Lambda
# para evaluación con vectores de datos xin
px = sym.lambdify(x,polinomio)

# Para graficar el polinomio en [a,b]
a = np.min(xi)
b = np.max(xi)
xin = np.linspace(a,b,muestras)
yin = px(xin)

# Usando evaluación simbólica
##yin = np.zeros(muestras,dtype=float)
##for j in range(0,muestras,1):
##    yin[j] = polinomio.subs(x,xin[j])
    
# SALIDA
print('Matriz Vandermonde: ')
print(D)
print('los coeficientes del polinomio: ')
print(coeficiente)
print('Polinomio de interpolación: ')
print(polinomio)
print('\n formato pprint')
sym.pprint(polinomio)

# Grafica
plt.plot(xi,fi,'o', label='[xi,fi]')
plt.plot(xin,yin, label='p(x)')
plt.xlabel('xi')
plt.ylabel('fi')
plt.legend()
plt.title(polinomio)
plt.show()

con lo que se obtiene el siguiente resultado:

Matriz Vandermonde: 
[[0.    0.    0.    1.   ]
 [0.008 0.04  0.2   1.   ]
 [0.027 0.09  0.3   1.   ]
 [0.064 0.16  0.4   1.   ]]
los coeficientes del polinomio: 
[[ 41.66666667]
 [-27.5       ]
 [  6.83333333]
 [  1.        ]]
Polinomio de interpolación: 
41.6666666666667*x**3 - 27.5*x**2 + 6.83333333333333*x + 1.0

formato pprint
                  3         2                           
41.6666666666667*x  - 27.5*x  + 6.83333333333333*x + 1.0

Para mostrar el polinomio de una manera más fácil de interpretar se usa la instrucción sym.pprint(), usada al final del algoritmo.

3.6.1 Método de Gauss-Seidel con Python

Referencia: Ejemplo 01 Chapra Ejemplo 11.3  p311 pdf335

\begin{cases} 3 x_0 -0.1 x_1 -0.2 x_2 = 7.85 \\ 0.1 x_0 +7x_1 -0.3 x_2 = -19.3 \\ 0.3 x_0 -0.2 x_1 +10 x_2 = 71.4 \end{cases}

El ejemplo de referencia, ya presenta una matriz pivoteada por filas, por lo que no fué implementado el procedimiento.


Planteamiento

Para el desarrollo con lápiz y papel, se despeja una variable de cada ecuación, se empieza con la primera expresión para obtener x0

3 x_0 -0.1 x_1 -0.2 x_2 = 7.85 3 x_0 = 7.85 + 0.1 x_1 + 0.2 x_2 x_0 = \frac{7.85 + 0.1 x_1 + 0.2 x_2}{3}

con la segunda ecuación se obtiene x1

0.1 x_0 +7x_1 -0.3 x_2 = -19.3 7x_1 = -19.3 - 0.1 x_0 +0.3 x_2 x_1 = \frac{ -19.3 - 0.1 x_0 +0.3 x_2}{7}

usando la tercera ecuación se obtiene x2

0.3 x_0 - 0.2 x_1 + 10 x_2 = 71.4 10 x_2 = 71.4 - 0.3 x_0 + 0.2 x_1 x_2 = \frac{71.4 - 0.3 x_0 + 0.2 x_1}{10}

Vector inicial

Como el vector inicial no se indica en el enunciado, se considera usar el vector de ceros para iniciar  las iteraciones.

X = [0,0,0]

con tolerancia de 0.00001


Iteraciones

Con las ecuaciones obtenidas en el planteamiento, se desarrolla usando el vector inicial presentado, hasta que el |error|<tolera

itera = 0

x_0 = \frac{7.85 + 0.1 (0) + 0.2 (0)}{3} = 2.61 x_1 = \frac{ -19.3 - 0.1 (2.61) +0.3 (0)}{7} = -2.79 x_2 = \frac{71.4 - 0.3 (2.61) + 0.2 (-2.79)}{10} = 7.00

X = [2.61, -2.79, 7.00]

diferencias = [2.61, -2.79, 7.00] – [0,0,0]
diferencias = [2.61, -2.79, 7.00]
error = ||diferencias|| , se usa la norma de max(diferencias)
error = 7.00

itera = 1

X = [2.61, -2.79, 7.00]

x_0 = \frac{7.85 + 0.1 (-2.79) + 0.2 (7.00)}{3} = 2.99 x_1 = \frac{ -19.3 - 0.1 (2.99) +0.3 (7.00)}{7} = -2.49 x_2 = \frac{71.4 - 0.3 (2.99) + 0.2 (-2.49)}{10} = 7.00

X = [2.99, -2.49, 7.00]

diferencias = [2.99, -2.49, 7.00] – [2.61, -2.79, 7.00]

diferencias = [0.38, 0.3 , 0. ]

error = ||diferencias||, se usa la norma de max(diferencias)

error = 0.38

itera = 2

X = [2.99, -2.49, 7.00]

x_0 = \frac{7.85 + 0.1 (-2.49) + 0.2 (7.00)}{3} = 3.00 x_1 = \frac{ -19.3 - 0.1 (3.00) +0.3 (7.00)}{7} = -2.5 x_2 = \frac{71.4 - 0.3 (3.00) + 0.2 (-2.5)}{10} = 7.00

X = [3.00, -2.50, 7.00]

diferencias = [3.00, -2.50, 7.00] – [2.99, -2.49, 7.00]

diferencias = [ 0.01, -0.01, 0. ]

error = ||diferencias||, se usa la norma de max(diferencias)

error = 0.01

El error aún es mayor que tolera, por lo que es necesario continuar con las iteraciones.

Observaciones

El error disminuye en cada iteración, por lo que el método converge hacia la raiz.


Algoritmo con Python

Con la descripción dada para el método de Gauss-Seidel, se muestra una forma básica de implementar el algoritmo.

Ejemplo 01 Chapra Ejemplo 11.3 p311 pdf335

\begin{cases} 3 x_0 -0.1 x_1 -0.2 x_2 = 7.85 \\ 0.1 x_0 +7x_1 -0.3 x_2 = -19.3 \\ 0.3 x_0 -0.2 x_1 +10 x_2 = 71.4 \end{cases}

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 iterativos

import 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
while not(errado<=tolera or itera>iteramax):
    # por fila
    for i in range(0,n,1):
        # por columna
        suma = 0 
        for j in range(0,m,1):
            # excepto diagonal de A
            if (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 converge
if (itera>iteramax):
    X=0
# revisa respuesta
verifica = np.dot(A,X)

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

que da como resultado:

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

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


Tarea

Completar el algoritmo para usar una matriz diagonal dominante, usando al menos el pivoteo parcial por filas.

3.6 Método de Gauss-Seidel

Referencia: Chapra 11.2 p310 pdf334, Burden 9Ed 7.3 p454, Rodriguez 5.2 p162

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.

\begin{cases} a_{0,0} x_0+a_{0,1}x_1+a_{0,2} x_2 = b_{0} \\ a_{1,0} x_0+a_{1,1}x_1+a_{1,2} x_2 = b_{1} \\ a_{2,0} x_0+a_{2,1}x_1+a_{2,2} x_2 = b_{2} \end{cases}

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:

x_0 = \frac{b_{0} -a_{0,1}x_1 -a_{0,2} x_2 }{a_{0,0}} x_1 = \frac{b_{1} -a_{1,0} x_0 -a_{1,2} x_2}{a_{1,1}} x_2 = \frac{b_{2} -a_{2,0} x_0 - a_{2,1} x_1}{a_{2,2}}

Observe que el patrón de cada fórmula para determinar x[i], tiene la forma:

x_i = \bigg(b_i - \sum_{j=0, j\neq i}^n A_{i,j}X_j\bigg) \frac{1}{A_{ii}}

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.

     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.

3.5.2 Número de Condición

Referencia: Chapra 10.3.2 p300/pdf324, Burden 9Ed 7.5 p470, Rodriguez 4.4.3 p133

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

 np.linalg.cond(A)

3.5.1 Normas de Vector o Matriz

Referencia: Chapra 10.3.1 p298 pdf322, Burden 9Ed Cap7.1 p432, Rodríguez 4.4.1 p132, MATG1049 – Norma, distancias y ángulos

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}

Ejercicios

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, 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

3.5 Normas de vector o matriz como distancias con Python

Referencia: Chapra 10.3.1 p298 pdf322, Burden 9Ed Cap7.1 p432, Rodriguez 4.4.1 p132

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

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

Distancia entre dos puntos en el espacio, error y Norma

Observe los siguientes videos, sobre acople de aeronaves.

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

  • ¿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:

X0 = np.array([0.0, 0, 0])
X1 = np.array([1.0, 2, 3])
X2 = np.array([2.0, 4,-1])

Las instrucciones gráficas principales son:

El resultado de la parte numérica se obtiene con pocas instrucciones en el bloque de procedimiento

X1 =  [1 2 3]
X2 =  [ 2  4 -1]
errado =  [-1 -2  4]
||errado|| =  4.58257569495584
Norma euclidiana :  4.58257569495584

las intrucciones en Python son:

# Norma como error
# o distancia entre dos puntos
# caso 3D
import numpy as np

# INGRESO
X0 = np.array([0.0, 0, 0])
X1 = np.array([1.0, 2, 3])
X2 = np.array([2.0, 4,-1])

# PROCEDIMIENTO
errado = X1 - X2
distancia = np.sqrt(np.sum(errado**2))
# funciones numpy
Nerrado = np.linalg.norm(errado)

# SALIDA
print('X1 = ', X1)
print('X2 = ', X2)
print('errado = ', errado)
print('||errado|| = ', distancia)
print('Norma euclidiana : ',Nerrado)


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

# puntos en el espacio
[x, y , z] = X0
grafica.scatter(x,y,z, c = 'blue',
                marker='o', label = 'X0')

[x, y , z] = X1
grafica.scatter(x,y,z, c = 'red',
                marker='o', label = 'X1')

[x, y , z] = X2
grafica.scatter(x,y,z, c = 'green',
                marker='o', label = 'X2')

# líneas entre puntos
linea = np.concatenate(([X0],[X1]),axis = 0)
x = linea[:,0]
y = linea[:,1]
z = linea[:,2]
grafica.plot(x,y,z, label = '||X1||')

linea = np.concatenate(([X0],[X2]),axis = 0)
x = linea[:,0]
y = linea[:,1]
z = linea[:,2]
grafica.plot(x,y,z, label = '||X2||')

linea = np.concatenate(([X1],[X2]),axis = 0)
x = linea[:,0]
y = linea[:,1]
z = linea[:,2]
grafica.plot(x,y,z, label = '||error||')

grafica.set_xlabel('eje x')
grafica.set_ylabel('eje y')
grafica.set_zlabel('eje z')
grafica.legend()

grafica.view_init(35, 25)
plt.show()