s3Eva_IT2018_T2 Drenaje de estanque

Ejercicio: 3Eva_IT2018_T2 Drenaje de estanque

literal a

Se usa interpolación para encontrar los polinomios que pasan por los puntos seleccionados.

El error de A(5) se obtiene como la diferencia entre el valor de la tabla y el polinomio del tramo [4,6] evaluado en el punto.

ordenado:  [6 5 4 3 2 1 0]
hi:  [0 1 2 3 4 5 6]
Ai:  [ 0.02  0.18  0.32  0.45  0.67  0.97  1.17]

puntos seleccionados:
h1:  [0, 2, 4, 6]
A1:  [ 0.02  0.32  0.67  1.17]

Polinomios por tramos: 
 x = [0,2]
0.000416666666666669*x**3 + 0.148333333333333*x + 0.02
 x = [2,4]
0.00416666666666666*x**3 - 0.0224999999999999*x**2 + 0.193333333333333*x - 0.00999999999999984
 x = [4,6]
-0.00458333333333333*x**3 + 0.0824999999999999*x**2 - 0.226666666666666*x + 0.549999999999999

error en px(5):  0.0637499999999998

se observa que la evaluación se realiza para el polinomio entre [4,6]

Desarrollo en Python

# 3ra Evaluación I Término 2018
# Tema 2. Drenaje de Estanque

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

def traza3natural(xi,yi):
    # Trazador cúbico natural, splines
    # resultado: polinomio en forma simbólica
    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 -------------------------

hi = np.array([6, 5, 4, 3, 2, 1, 0])
Ai = np.array([1.17, 0.97, 0.67, 0.45, 0.32, 0.18, 0.02])
xk = 5

# PROCEDIMIENTO LITERAL a
# reordena en forma ascendente
ordenado = np.argsort(hi)
hi = hi[ordenado]
Ai = Ai[ordenado]

# Selecciona puntos
xi = [0,2,4,6]
fi = Ai[xi]
n = len(xi)

polinomio = traza3natural(xi,fi)

# literal a, estima error
px = polinomio[2]
pxk = px.subs('x',xk)
errado = np.abs(Ai[xk] - pxk)

# SALIDA
print('ordenado: ', ordenado)
print('hi: ', hi)
print('Ai: ', Ai)
print('puntos seleccionados:')
print('h1: ', xi)
print('A1: ', fi)

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

print('error en px(5): ', errado)

# GRAFICA
# Puntos para grafica en cada tramo
resolucion = 10 # entre cada par de puntos
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

# GRAFICA
# puntos originales
plt.plot(hi,Ai,'o',label = 'Ai')
# Trazador cúbico
plt.plot(xtrazado,ytrazado, label = 'p(h)')
plt.plot(xi,fi,'o', label = 'Apx')
plt.title('Trazador cúbico natural (splines)')
plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.grid()
plt.show()

Literal b

TAREA