Mascota 50 años

Ejercicio de interpolación

  • Desarrollar por grupos
  • integrar soluciones en un mismo gráfico
  • mostrar polinomios y rangos de existencia

# caparazon superior
xi = [113., 117, 134, 153, 169, 184, 194, 203]
fi = [127., 141, 161, 166, 160, 155, 140, 132]
puntos = [[xi,fi]]

# caparazon inferior 1
xi = [113., 123, 130, 149, 182, 197, 208, 211]
fi = [127., 116, 112, 107, 110, 114, 112, 108]
puntos.append([xi,fi])

# caparazon inferior 2
xi = [107., 114, 120, 143, 170, 192, 210]
fi = [124., 114, 108, 99, 99, 106, 107]
puntos.append([xi,fi])

# Patas 01
xi = [110., 116, 120, 121, 134, 143, 153, 168, 173, 177, 181, 188, 194, 201, 205, 210, 214, 218]
fi = [92., 86, 89, 86, 87, 95, 91, 93, 85, 85, 83, 86, 88, 86, 91, 87, 92, 92]
puntos.append([xi,fi])

# Patas 02
xi = [109., 115, 117, 120, 125, 130, 134, 138, 142, 145, 150, 152, 154]
fi = [91., 87, 82, 82, 78, 78, 81, 78, 76, 79, 77, 80, 86]
puntos.append([xi,fi])

# Patas 03
xi = [172., 175, 182, 186, 191, 195, 201, 205, 210, 211, 217]
fi = [86., 79, 79, 77, 82, 78, 81, 79, 82, 85, 90]
puntos.append([xi,fi])

# Patas 04
xi = [113., 118, 122, 127, 132, 136, 140, 144, 149, 152]
fi = [88., 90, 88, 86, 85, 83, 87, 84, 86, 85]
puntos.append([xi,fi])

# Rabito 01
xi = [97., 102, 108, 111, 117, 121]
fi = [120., 113, 108, 105, 102, 102]
puntos.append([xi,fi])

# cabeza 01
xi = [194., 196, 203, 210, 211, 209]
fi = [177., 167, 157, 149, 138, 135]
puntos.append([xi,fi])

# cabeza 02
xi = [195., 199, 207, 214, 219, 224, 229, 234, 239, 242, 244]
fi = [177., 185, 190, 190, 188, 193, 192, 185, 175, 172, 168]
puntos.append([xi,fi])

# cabeza 03
xi = [220., 226, 234, 239, 241]
fi = [164., 162, 163, 163, 163]
puntos.append([xi,fi])

# cabeza 04
xi = [203., 211, 214, 219, 223, 224, 225, 230, 236, 241]
fi = [115., 119, 125, 124, 127, 137, 146, 149, 154, 162]
puntos.append([xi,fi])

# cabeza 05
xi = [208., 212, 215, 219, 221, 225, 228]
fi = [174., 177, 178, 179, 182, 183, 178]
puntos.append([xi,fi])

# cabeza 06
xi = [206., 210, 214, 218, 220, 224, 229, 233, 232]
fi = [179., 182, 182, 181, 174, 171, 170, 175, 180]
puntos.append([xi,fi])

Mascota descansando

Referencia: Burden 9th Ed. Ejercicio 32 p164

La parte superior de una mascota descansando se describe con los puntos presentados en la tabla.  Los puntos se usan para ejercicios de interpolación.

xi = [1,2,5,6,7,8,10,13,17,20,23,24,25,27,27.27,28,29,30]
fi = [3.0,3.7,3.9,4.2,5.7,6.6,7.1,6.7,4.5,7.0,6.1,5.6,5.8,5.2,4.1,4.3,4.1,3.0]

5.4.1 Trazadores Cúbicos

Referencia:  Burden 3.4 p.141 pdf.151, Chapra 18.6.3 p.532 pdf.556

Se busca obtener un polinómio de tercer grado para cada intervalo entre nodos:
S_j(x_j) = a_j + b_j (x-x_j) + c_j(x-xj)^2 + d_j(x-x_j)^3

para cada j = 0,1,…, n-1

Para n datos, existen n-1 intervalos, cuatro incógnitas que evaluar por cada intervalo, como resultado 4(n-1) incógnitas.

Para los términos (xj+1xj) que se usan varias veces en el desarrollo, se simplifican a:

h_j = x_{j+1} - x_{j}

Si an = f(xn):

a_{j+1} = a_j + b_j h_j + c_j h_j^2 + d_jh_j^3

Las condiciones para los sistemas de ecuaciones son:

1.  Los valores de la función deben ser iguales en los nodos interiores

S_j(x_{j+1}) = f(x_{j+1}) = S_{j+1}(x_{j+1})

2. Las primeras derivadas en los nodos interiores deben ser iguales

S'_j(x_{j+1}) = S'_{j+1}(x_{j+1})

3. Las segundas derivadas en los nodos interiores deben ser iguales

S''_j(x_{j+1}) = S''_{j+1}(x_{j+1})

4. La primera y última función deben pasar a través de los puntos extremos

5. Una de las condiciones de frontera se satisface:

5a. frontera libre o natural: Las segundas derivadas en los nodos extremos son cero

S''(x_0) = S''(x_n) = 0

5b. frontera sujeta: las primeras derivadas en los nodos extremos son conocidas

S'(x_0) = f'(x_0)

S'(x_n) = f'(x_n)

Siguiendo el desarrollo propuesto, se obtiene el siguiente algoritmo:

# Trazador cúbico natural
import numpy as np
import sympy as sym

def traza3natural(xi,yi):
    n = len(xi)
    
    # Valores h
    h = np.zeros(n-1, dtype = float)
    for j in range(0,n-1,1):
        h[j] = xi[j+1] - xi[j]
    
    # Sistema de ecuaciones
    A = np.zeros(shape=(n-2,n-2), dtype = float)
    B = np.zeros(n-2, dtype = float)
    S = np.zeros(n, dtype = float)
    A[0,0] = 2*(h[0]+h[1])
    A[0,1] = h[1]
    B[0] = 6*((yi[2]-yi[1])/h[1] - (yi[1]-yi[0])/h[0])
    for i in range(1,n-3,1):
        A[i,i-1] = h[i]
        A[i,i] = 2*(h[i]+h[i+1])
        A[i,i+1] = h[i+1]
        B[i] = 6*((yi[i+2]-yi[i+1])/h[i+1] - (yi[i+1]-yi[i])/h[i])
    A[n-3,n-4] = h[n-3]
    A[n-3,n-3] = 2*(h[n-3]+h[n-2])
    B[n-3] = 6*((yi[n-1]-yi[n-2])/h[n-2] - (yi[n-2]-yi[n-3])/h[n-3])
    
    # Resolver sistema de ecuaciones
    r = np.linalg.solve(A,B)
    # S
    for j in range(1,n-1,1):
        S[j] = r[j-1]
    S[0] = 0
    S[n-1] = 0
    
    # Coeficientes
    a = np.zeros(n-1, dtype = float)
    b = np.zeros(n-1, dtype = float)
    c = np.zeros(n-1, dtype = float)
    d = np.zeros(n-1, dtype = float)
    for j in range(0,n-1,1):
        a[j] = (S[j+1]-S[j])/(6*h[j])
        b[j] = S[j]/2
        c[j] = (yi[j+1]-yi[j])/h[j] - (2*h[j]*S[j]+h[j]*S[j+1])/6
        d[j] = yi[j]
    
    # Polinomio trazador
    x = sym.Symbol('x')
    polinomio = []
    for j in range(0,n-1,1):
        ptramo = a[j]*(x-xi[j])**3 + b[j]*(x-xi[j])**2 + c[j]*(x-xi[j])+ d[j]
        ptramo = ptramo.expand()
        polinomio.append(ptramo)
    
    return(polinomio)

# PROGRAMA de prueba
# INGRESO , Datos de prueba
xi = np.array([0.1 , 0.2, 0.3, 0.4])
fi = np.array([1.45, 1.8, 1.7, 2.0])
resolucion = 10 # entre cada par de puntos

# PROCEDIMIENTO
n = len(xi)
# Obtiene los polinomios por tramos
polinomio = traza3natural(xi,fi)

# SALIDA
print('Polinomios por tramos: ')
for tramo in range(1,n,1):
    print(' x = ['+str(xi[tramo-1])+','+str(xi[tramo])+']')
    print(str(polinomio[tramo-1]))

con lo que se obtiene el resultado por cada tramo

Polinomios por tramos: 
 x = [0.1,0.2]
-146.666666666667*x**3 + 44.0*x**2 + 0.566666666666666*x + 1.1
 x = [0.2,0.3]
283.333333333333*x**3 - 214.0*x**2 + 52.1666666666667*x - 2.34
 x = [0.3,0.4]
-136.666666666667*x**3 + 164.0*x**2 - 61.2333333333333*x + 9.0
>>>

Para el trazado de la gráfica se añade al final del algoritmo:

# GRAFICA
# Puntos para grafica en cada tramo
xtrazado = np.array([])
ytrazado = np.array([])
tramo = 1
while not(tramo>=n):
    a = xi[tramo-1]
    b = xi[tramo]
    xtramo = np.linspace(a,b,resolucion)
    
    ptramo = polinomio[tramo-1]
    pxtramo = sym.lambdify('x',ptramo)
    ytramo = pxtramo(xtramo)
    
    xtrazado = np.concatenate((xtrazado,xtramo))
    ytrazado = np.concatenate((ytrazado,ytramo))
    tramo = tramo + 1

# Gráfica
import matplotlib.pyplot as plt
plt.title('Trazador cúbico natural (splines)')
plt.plot(xtrazado,ytrazado)
plt.plot(xi,fi,'o')
plt.show()

5.4 Trazadores lineales (Splines) grado1

Referencia: Chapra 18.6.1 p.525 pdf.549

El concepto de trazador se originó en la técnica de dibujo que usa una cinta delgada y flexible (spline) para dibujar curvas suaves a través de un conjunto de puntos.

La unión más simple entre dos puntos es una línea recta. Los trazadores de primer grado para un grupo de datos ordenados pueden definirse como un conjunto de funciones lineales.

f(x) = f(x_0) + m_0(x-x_0), x0\leq x\leq x_1 f(x) = f(x_1) + m_1(x-x_1), x1\leq x\leq x_2

f(x) = f(x_{n-1}) + m_{n-1}(x-x_{n-1}), x{n-1}\leq x\leq x_n

donde

m_i = \frac{f(x_{i+1}) - f(x_i)}{(x_{i+1}-x_i)}

Las ecuaciones se pueden usar para evaluar la función en cualquier punto entre x0 y xn localizando primero el intervalo dentro del cual está el punto.

Aplicando el método como algoritmo Python:

# Trazador (spline) lineal, grado 1
import numpy as np
import sympy as sp

def trazalineal(xi,fi):
    n = len(xi)
    x = sp.Symbol('x')
    polinomio = []
    tramo=1
    while not(tramo>=n):
        m =(fi[tramo]-fi[tramo-1])/(xi[tramo]-xi[tramo-1])
        inicio = fi[tramo-1]-m*xi[tramo-1]
        ptramo = inicio + m*x
        polinomio.append(ptramo)
        tramo = tramo + 1
    return(polinomio)

# PROGRAMA
# INGRESO , Datos de prueba
xi = [0.1 , 0.2, 0.3, 0.4]
fi = [1.45, 1.8, 1.7, 2.0]
resolucion = 10 # entre cada par de puntos

# PROCEDIMIENTO
n = len(xi)
# Obtiene los polinomios por tramos
polinomio = trazalineal(xi,fi)

# SALIDA
print('Polinomios por tramos: ')
for tramo in range(1,n,1):
    print(' x = ['+str(xi[tramo-1])
          +','+str(xi[tramo])+']')
    print(str(polinomio[tramo-1]))

Se obtiene como resultado:

Polinomios por tramos: 
 x = [0.1,0.2]
3.5*x + 1.1
 x = [0.2,0.3]
-1.0*x + 2.0
 x = [0.3,0.4]
3.0*x + 0.8
>>> 

Para la gráfica se añade:

# GRAFICA
# Puntos para graficar cada tramo
xtrazado = np.array([])
ytrazado = np.array([])
tramo = 1
while not(tramo>=n):
    a = xi[tramo-1]
    b = xi[tramo]
    xtramo = np.linspace(a,b,resolucion)
    
    ptramo = polinomio[tramo-1]
    pxtramo = sp.lambdify('x',ptramo)
    ytramo = pxtramo(xtramo)
    
    xtrazado = np.concatenate((xtrazado,xtramo))
    ytrazado = np.concatenate((ytrazado,ytramo))
    tramo = tramo + 1

# Gráfica
import matplotlib.pyplot as plt
plt.title('Trazadores lineales (splines)')
plt.plot(xtrazado,ytrazado)
plt.plot(xi,fi,'o')
plt.show()

5.5 Interpolación paramétrica

Referencia: Rodriguez 6.9.2 pdf236

Si los datos (x,y) no tienen una relación de tipo funcional y(x), entonces no se pueden aplicar directamente los métodos de interpolación revisados.

Sin embargo si las coordenadas (x,y) se expresan como funciones de otra variable t denomiada parámetro, entonces los puntos x(t), y(t) tienen relación funcional, y se pueden construir polinomios de interpolación.

Ejemplo: las coordinadas x(t) y y(t) del recorrido de un cohete registradas en los instantes t fueron:

# interpolación paramétrica
import numpy as np
import matplotlib.pyplot as plt

ti  = [0,1,2,3]
xti = [2,1,3,4]
yti = [0,4,5,0]

# PROCEDIMIENTO
# interpolando con lagrange
px = lambda t: (-2/3)*(t**3) + (7/2)*(t**2) + (-23/6)*t + 2
py = lambda t: (-1/2)*(t**3) + (9/2)*t

t =  np.arange(0,3,0.01)
puntosx = px(t)
puntosy = py(t)

# Salida
plt.plot(puntosx,puntosy)
plt.show()

5.1.2 Diferencias divididas – Newton

Referencia: Chapra 10.1.3 Pdf.532, Rodriguez 6.7 Pdf.223

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-x_{0}}

Para lo cual debe interpretar la tabla de diferencias divididas:

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]

Las diferencias sirven para evaluar los coeficientes y obtener el polinomio de interpolacion:

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]

Se conoce como el polinomio de interpolación de Newton en diferencias divididas.

Tarea: Realizar el algoritmo para implementar el polinomio usando la interpolación de Newton en diferencias divididas

5.1.1 Diferencias finitas avanzadas polinomio

Referencia: Rodriguez 6.6.4 Pdf.221

En la tabla de diferencias finitas,  si los valores en el eje x se encuentran igualmente espaciados, la diferencia es una constante denominada 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 z en el intervalo [x0,xn]

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

Polinomio de interpolación de diferencias finitas avanzadas:

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

El polinomio se puede realizar usando el ejercicio de diferencias finitas. Se crea la variable simbólica con sympy, se genera el polinomio añadiendo los términos de la operación.

Para simplificar se expanden las operaciones de multiplicación, y se obtiene el polinomio.

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

# Diferencias finitas avanzadas para polinomio interpolación
# Referencia Rodriguez 6.6.4 Pdf.221
# Tarea: Verificar tamaño de vectores
#        considerar puntos no equidistantes en eje x
import numpy as np
import matplotlib.pyplot as plt
import sympy as sym

# 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)
# cambia a forma de columnas
i = np.arange(1,n+1,1)
i = np.transpose([i])
xi = np.transpose([xi])
fi = np.transpose([fi])
# Añade matriz de diferencias
dfinita = np.zeros(shape=(n,n),dtype=float)
tabla = np.concatenate((i,xi,fi,dfinita), axis=1)
# Sobre matriz de diferencias, por columnas
[n,m] = np.shape(tabla)
c = 3
diagonal = n-1
while (c<m):
    # Aumenta el título para cada columna
    titulo.append('df'+str(c-2))
    # calcula cada diferencia por fila
    f = 0
    while (f < diagonal):
        tabla[f,c] = tabla[f+1,c-1]-tabla[f,c-1]
        f = f+1
    
    diagonal = diagonal - 1
    c = c+1

# POLINOMIO con diferencias finitas
# caso: puntos en eje x equidistantes
dfinita = tabla[:,3:]
n = len(dfinita)
x = sym.Symbol('x')
h = xi[1,0]-xi[0,0]
polinomio = fi[0,0]
for c in range(1,n,1):
    denominador = np.math.factorial(c)*(h**c)
    factor = dfinita[0,c-1]/denominador
    termino=1
    for f  in range(0,c,1):
        termino = termino*(x-xi[f])
    polinomio = polinomio + termino*factor
# simplifica polinomio, multiplica los (x-xi)
polinomio = polinomio.expand()
# para evaluacion numérica
px = sym.lambdify(x,polinomio)

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

# SALIDA
print([titulo])
print(tabla)
print('polinomio:')
print(polinomio)

# Gráfica
plt.title('Interpolación polinómica')
plt.plot(xi,fi,'o', label = 'Puntos')
plt.plot(xi_p,fi_p, label = 'Polinomio')
plt.legend()
plt.show()

teniendo como resultado:

[['i', 'xi', 'fi', 'df1', 'df2', 'df3', 'df4']]
[[ 1.    0.1   1.45  0.15 -0.05  0.25  0.  ]
 [ 2.    0.2   1.6   0.1   0.2   0.    0.  ]
 [ 3.    0.3   1.7   0.3   0.    0.    0.  ]
 [ 4.    0.4   2.    0.    0.    0.    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 expandido
41.6666666666667*x**3 - 27.5*x**2 + 6.83333333333335*x + 0.999999999999999
>>> 

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

>>> p(0.1)
1.4500000000000004
>>> p(0.2)
1.6000000000000025
>>> p(0.3)
1.7000000000000042
>>> 

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

5.1 Diferencias finitas

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)

Por ejemplo, si los puntos son:

xi = [0.10, 0.2, 0.3, 0.4]
fi = [1.45, 1.6, 1.7, 2.0]

la tabla será:

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

Cada casila de diferencia finita se obtiene restando los dos valores consecutivos de la columna anterior.

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.

Para crear la tabla mediante un algoritmo en Python, se crean primero las columnas que contienen, el índice i, los valores xi y fi.

Con las tres columnas iniciales, se crea una matriz aumentada con la matriz de diferencias, que tiene tamaño nxn.

Se calculan las diferencias para cada columna, realizando la operación entre los valores de cada fila.
Finalmente se muestra el resultado.

# Diferencias finitas
# Chapra 18.1.3 pdf.509; Rodriguez 6.5 Pdf.211
# Genera la 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)
# cambia a forma de columnas
i = np.arange(1,n+1,1)
i = np.transpose([i])
xi = np.transpose([xi])
fi = np.transpose([fi])
# Añade matriz de diferencias
dfinita = np.zeros(shape=(n,n),dtype=float)
tabla = np.concatenate((i,xi,fi,dfinita), axis=1)
# Sobre matriz de diferencias, por columnas
[n,m] = np.shape(tabla)
diagonal = n-1

c = 3
while (c<m):
    # Aumenta el título para cada columna
    titulo.append('df'+str(c-2))
    # calcula cada diferencia por fila
    f = 0
    while (f < diagonal):
        tabla[f,c] = tabla[f+1,c-1]-tabla[f,c-1]
        f = f+1
    
    diagonal = diagonal - 1
    c = c+1

# SALIDA
print([titulo])
print(tabla)

el resultado de aplicar el algoritmo es:

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

>>> 

Pato en pleno vuelo

Referencia: Burden Cap 3. Ejemplo 1.

La figura muestra un joven pato en pleno vuelo. Para aproximar el perfil de la parte superior del pato, se presentan 21 puntos a lo largo de la curva de aproximación relativos a un sistema de coordenadas sobrepuestas.

xi = [0.9, 1.3, 1.9, 2.1, 2.6, 3.0, 3.9, 4.4, 4.7, 5, 6.0, 7.0, 8.0, 9.2, 10.5, 11.3, 11.6, 12.0, 12.6, 13.0, 13.3]
fi = [1.3, 1.5, 1.85, 2.1, 2.6, 2.7, 2.4, 2.15, 2.05, 2.1, 2.25, 2.3, 2.25, 1.95, 1.4, 0.9, 0.7, 0.6, 0.5, 0.4, 0.25]


a) Realice la interpolación entre puntos usando la interpolación de Lagrange para la parte superior del pato en vuelo.

b) ¿Es posible crear un solo polinomio para los 21 puntos presentados? Explique su respuesta en pocas palabras.

c) ¿En concepto, cuando considera necesario crear un nuevo polinomio?

d) Adjunte en los archivos para: la gráfica que interpola los puntos, y los polinomios que la generan.

e) (Puntos extra: 10, para promediar con lección o taller)

Se adjuntan los demás puntos del perfil del pato en pleno vuelo, presente el o los polinomios que interpolan la figura (gráfica y polinomios), las despuestas a literales b,c,d deben también ser válidas para el ejercicio completo.

# Perfil Superior del pato
xiA = [0.9, 1.3, 1.9, 2.1, 2.6, 3.0, 3.9, 4.4, 4.7, 5, 6.0, 7.0, 8.0, 9.2, 10.5, 11.3, 11.6, 12.0, 12.6, 13.0, 13.3]
fiA = [1.3, 1.5, 1.85, 2.1, 2.6, 2.7, 2.4, 2.15, 2.05, 2.1, 2.25, 2.3, 2.25, 1.95, 1.4, 0.9, 0.7, 0.6, 0.5, 0.4, 0.25]

# Perfil inferior cabeza
xiB = [0.817, 0.897, 1.022, 1.191, 1.510, 1.834, 2.264, 2.962, 3.624, 4.202, 4.499, 4.779, 5.109, 5.527]
fiB = [1.180, 1.065, 1.023, 1.010, 1.032, 1.085, 1.192, 1.115, 1.087, 1.100, 0.830, 0.608, 0.350, 0.106]

# Perfil Ala superior
xiC = [4.659, 4.865, 5.085, 5.261, 5.387, 5.478, 5.527]
fiC = [-5.161, -4.741, -3.933, -2.951, -1.970, -0.981, 0.106]

# Perfil Ala inferior
xiD = [4.659, 4.750, 4.990, 5.289, 5.560, 5.839, 6.113, 6.606, 6.916, 7.305, 7.563, 7.802, 7.983]
fiD = [-5.161, -5.259, -5.284, -5.268, -5.161, -4.982, -4.769, -4.286, -3.911, -3.213, -2.670, -2.176, -1.655]

# Perfil inferior posterior
xiE = [8.141, 8.473, 8.832, 9.337, 9.887, 10.572, 10.995, 11.501, 11.923, 12.364, 12.763, 13.300]
fiE = [-1.138, -0.434, -0.514, -0.494, -0.382, -0.005, -0.090, -0.085, -0.030, 0.093, 0.120, 0.250]

# Perfil ojo superior
xiF = [2.663, 2.700, 2.805, 2.886]
fiF = [2.202, 2.279, 2.293, 2.222]

# Perfil ojo inferior
xiG = [2.663, 2.720, 2.826, 2.886]
fiG = [2.202, 2.130, 2.143, 2.222]

5.3 Interpolación de Lagrange

Referencia: Chapra 18.2 pdf 540

El polinomio de interpolación de lagrange es una reformulación del polinomio de Newton que evita el cálculo de las diferencias divididas. Se crea como:

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}

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

p(x) = a0x3 + a1x2 + a2x1 + a3

xi 0 0.2 0.3 0.4
fi 1 1.6 1.7 2.0
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) + 8.33333333333333*(-5.0*x + 1.0)*(x - 0.4)*(x - 0.3)

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

Propuesta de algoritmo usando Python, como tarea realice la gráfica del polinomio resultante.

# Interpolacion de Lagrange
# Polinomio en forma simbólica
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
n = len(xi)
x = sym.Symbol('x')
# Polinomio
polinomio = 0
for i in range(0,n,1):
    # Termino de Lagrange
    termino = 1
    for j  in range(0,n,1):
        if (j!=i):
            termino = termino*(x-xi[j])/(xi[i]-xi[j])
    polinomio = polinomio + termino*fi[i]
# Expande el polinomio
px = polinomio.expand()
# para evaluacion numérica
pxn = sym.lambdify(x,polinomio)

# Puntos para la gráfica
a = np.min(xi)
b = np.max(xi)
muestras = 101
xi_p = np.linspace(a,b,muestras)
fi_p = pxn(xi_p)

# Salida
print('Polinomio de Lagrange, expresiones')
print(polinomio)
print()
print('Polinomio de Lagrange: ')
print(px)

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