4.4 Trazadores Cúbicos Natural o libre con Python

Referencia:  Burden 3.5 p105, Chapra 18.6.3 p532, Rodríguez 6.11.1 p244

Tiene como objetivo incorporar condiciones adicionales para la función en los extremos del intervalo donde se encuentran los puntos o «nodos».

Por ejemplo, si los datos son los de una trayectoria en un experimento de física, se podría disponer de la aceleración el punto inicial de las muestras (izquierda) y de salida (derecha) del intervalo de las muestras.

El método busca obtener un polinomio de tercer grado para cada sub-intervalo o tramo entre «nodos» consecutivos (j,  j+1), de la forma:

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 tramos, cuatro incógnitas (coeficientes) que evaluar por cada tramo, como resultado 4(n-1) incógnitas para todo el intervalo .

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

h_j = x_{j+1} - x_{j} S_j(x_j) = a_j + b_j h_j + c_j h_j^2 + d_jh_j^3

Para generar el sistema de ecuaciones, se siguen los siguientes criteros:

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. El primer y último polinomio deben pasar a través de los puntos extremos

f(x_0) = S_0(x_0) = a_0 f(x_n) = S_n(x_n) = a_n

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)

Tarea: Revisar en los textos el planteamiento de las ecuaciones para resolver el sistema que se genera al plantear el polinomio.

Ubicar en el texto las ecuaciones resultantes, que son las que se aplicarán en el algoritmo.


Algoritmo en Python

El algoritmo parte de lo desarrollado para «trazadores lineales», donde se presentaron los bloques para:

  • construir el trazador en una tabla de polinomios por tramos
  • graficar los trazadores en cada tramo,
  • evaluar cada uno de ellos en cada tramo con muestras suficientes para presentar una buena precisión en la gráfica

Del algoritmo básico se modifica entonces el bloque del cálculo de los polinomios de acuerdo al planteamiento de formulas y procedimientos para trazadores cúbicos naturales (enviado a revisar como tarea).

# Trazador cúbico natural
# Condición: S''(x_0) = S''(x_n) = 0
import numpy as np
import sympy as sym
import matplotlib.pyplot as plt

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]
        factor21 = (yi[i+2]-yi[i+1])/h[i+1]
        factor10 = (yi[i+1]-yi[i])/h[i]
        B[i] = 6*(factor21 - factor10)
        
    A[n-3,n-4] = h[n-3]
    A[n-3,n-3] = 2*(h[n-3]+h[n-2])
    factor12 = (yi[n-1]-yi[n-2])/h[n-2]
    factor23 = (yi[n-2]-yi[n-3])/h[n-3]
    B[n-3] = 6*(factor12 - factor23)
    
    # Resolver sistema de ecuaciones S
    r = np.linalg.solve(A,B)
    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
        factor10 = (yi[j+1]-yi[j])/h[j]
        c[j] = factor10 - (2*h[j]*S[j]+h[j]*S[j+1])/6
        d[j] = yi[j]
    
    # Polinomio trazador
    x = sym.Symbol('x')
    px_tabla = []
    for j in range(0,n-1,1):

        pxtramo = a[j]*(x-xi[j])**3 + b[j]*(x-xi[j])**2
        pxtramo = pxtramo + c[j]*(x-xi[j])+ d[j]
        
        pxtramo = pxtramo.expand()
        px_tabla.append(pxtramo)
    
    return(px_tabla)

# PROGRAMA -----------------------
# 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])
muestras = 10 # entre cada par de puntos

# PROCEDIMIENTO
# Tabla de polinomios por tramos
n = len(xi)
px_tabla = 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(px_tabla[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 graficar cada tramo
xtraza = np.array([])
ytraza = np.array([])
tramo = 1
while not(tramo>=n):
    a = xi[tramo-1]
    b = xi[tramo]
    xtramo = np.linspace(a,b,muestras)
    
    # evalua polinomio del tramo
    pxtramo = px_tabla[tramo-1]
    pxt = sym.lambdify('x',pxtramo)
    ytramo = pxt(xtramo)

    # vectores de trazador en x,y
    xtraza = np.concatenate((xtraza,xtramo))
    ytraza = np.concatenate((ytraza,ytramo))
    tramo = tramo + 1

# Gráfica
plt.plot(xi,fi,'ro', label='puntos')
plt.plot(xtraza,ytraza, label='trazador'
         , color='blue')
plt.title('Trazadores Cúbicos Naturales')
plt.xlabel('xi')
plt.ylabel('px(xi)')
plt.legend()
plt.show()

Si los polinomios no se igualan entre los nodos, tampoco sus velocidades y aceleraciones; podría considerar la experiencia semejante a variaciones de velocidad y aceleración en una trayectoria como la mostrada en los siguientes videos.

Car sales woman scares customers. maxman.tv. 4 enero 2016