Ejercicio: 1Eva_2021PAOII_T1 Interpolación para perfil de terreno
a. Plantee y desarrolle un polinomio P3(d1) de grado 3, que describa el perfil del terreno en el intervalo [0,1300] de distancias a la primera antena d1.
Para plantear el ejercicio, se observan los punto en el intervalo [0,1300], asi como los tamaños de paso.
Como el polinomio es de grado 3,se requiere 3 tramos o 4 muestras, el intervalo presenta 5. Para incluir el punto de inicio, fin y máximo, las muestras seleccionadas son las que se muestran en la tabla siguiente:
| Δd1 |
distancia d1 |
Perfil de Terreno |
DifDiv1 |
DifDiv2 |
DifDiv3 |
| 350 |
0 |
85 |
DifDiv11 = 0.02857 |
DifDiv12 = -6.122E-05 |
DifDiv13 = -1.123E-05 |
| 350 |
350 |
95 |
DifDiv21 = -0.01428 |
DifDiv22 = -1.127E-05 |
|
| 600 |
700 |
90 |
DifDiv31 =
-0.025 |
|
|
|
1300 |
75 |
|
|
|
Los tamaños de paso, Δd1, son diferentes, por lo que se pueden usar el método de diferencias divididas de Newton o el Método de Lagrange.
Usando diferencias divididas de Newton, se completa la tabla con las expresiones:
DifDiv_{11} = \frac{95-85}{350-0} = 0.02857143
DifDiv_{21} = \frac{90-95}{700-350} = -0.01428
DifDiv_{31} = \frac{75-90}{1300-700} =-0.025
DifDiv_{12} = \frac{-0.01428-0.02857}{700-0} = -6.122E-05
DifDiv_{22} = \frac{-0.025-(-0.01428)}{1300-350} = -1.127E-05
DifDiv_{13} = \frac{(-1.127E-05) -(6.122E-05)}{1300-0} = 1.123E-05
con lo que se puede construir el polinomio
p_{3}(d) = 85 + 0.02857(d-0)+(-6.122E5)(d-0)(d-350)+
-(1.123E-5)(d-0)(d-350)(d-700)
simplificando el polinomio se obtiene:
p_{3}(d) = 85.0+ 0.05941 d - 0.0001015 d^2 +(3.842E-8)d^3
b. Calcule el error sobre el o los datos que no se usaron en el intervalo
El valor no usado que estaba en el intervalo es d=1000, que se evalúa en p3(d) y se compara con el valor muestreado para ver el error del modelo
p_{3}(1000) = 85.0+ 0.05941 (1000) - 0.0001015 (1000)^2
+(3.842E-8)(1000)^3
p_{3}(1000) = 81.26
Error = |81.26 - 80| = 1.26
encontrando que el error es de 1.26m de altitud del terreno a los 1000m de distancia desde la referencia a la izquierda.
c. Desarrolle y justifique una propuesta para disminuir los errores encontrados en el literal anterior, sobre el mismo intervalo, es decir obtiene un nuevo polinomio (use algoritmo).
En caso de requerir una mayor precisión, se puede aumentar el grado del polinomio al incluir el punto d=1000 no usado en el paso anterior.
Usando el algoritmo con Python se obtiene:
Tabla Diferencia Dividida
[['i ', 'xi ', 'fi ', 'F[1]', 'F[2]', 'F[3]', 'F[4]', 'F[5]']]
[[ 0.0000e+00 0.0000e+00 8.5000e+01 2.8571e-02 -6.1224e-05 3.1920e-08 2.1666e-11 0.0000e+00]
[ 1.0000e+00 3.5000e+02 9.5000e+01 -1.4286e-02 -2.9304e-05 6.0086e-08 0.0000e+00 0.0000e+00]
[ 2.0000e+00 7.0000e+02 9.0000e+01 -3.3333e-02 2.7778e-05 0.0000e+00 0.0000e+00 0.0000e+00]
[ 3.0000e+00 1.0000e+03 8.0000e+01 -1.6667e-02 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00]
[ 4.0000e+00 1.3000e+03 7.5000e+01 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00]]
dDividida:
[ 2.8571e-02 -6.1224e-05 3.1920e-08 2.1666e-11 0.0000e+00]
polinomio:
2.16658863275405e-11*x*(x - 1000.0)*(x - 700.0)*(x - 350.0) + 3.19204604918891e-8*x*(x - 700.0)*(x - 350.0) - 6.12244897959184e-5*x*(x - 350.0) + 0.0285714285714286*x + 85.0
polinomio simplificado:
2.16658863275405e-11*x**4 - 1.24946064795689e-8*x**3 - 6.6683650518237e-5*x**2 + 0.0525123706702654*x + 85.0
junto a la gráfica:
d. Escriba sus conclusiones y recomendaciones sobre los resultados obtenidos entre los dos polinomios.
Dado que el error se considera mínimo en el intervalo, y el coeficiente del polinomio para x4 es del orden de 10-11 se recomendaría usar el polinomio de grado 3.
Instrucciones en Python
# 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([0.0,350,700,1000,1300])
fi = np.array([85.0,95,90,80,75])
# 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()
