Estaba intentando graficar un Método de Splines (trazadores cúbicos) y obtuve las distintas funciones para cada segmento. El tema es que cuando intento graficar creo que la función S sólo grafica una de las funciones Spline (creo que es S0) Y necesito que los puntos y y las curvas Spline coincidan. ¿Que error estoy cometiendo?
Aquí va el planteo de los Trazadores cúbicos:
Interpolación segmentada: Interpolación por "splines" En el segundo caso, aplicamos la Interpolación por Splines, es decir, interpolar cada segmento de puntos con trazadores cúbicos (polinomios de grado 3), donde hay n polinomios por cada n+1 por definición.
Para este polinomio existen cuatro incógnitas: ai,bi,ci,di.
Si(x) = ai + bi(x − xi) + ci(x − xi)2 + di(x − xi)3
A las curvas se les impone ciertas condiciones:
1) Condición de interpolación: que el trazador valga lo mismo que la función en el punto.
Si(xi) = f(xi) para cada i = 0; 1; . . . ; n
2) Continuidad en la unión de trazadores.
Si+1(xi+1) = Si(xi+1) para cada i = 0; 1; . . . ; n − 2
3) Continuidad de derivadas.
S'i+1(xi+1) = S'i(xi+1) para cada i = 0; 1; . . . ; n − 2
4) Continuidad de la segunda derivada: que tenga la misma concavidad y no haya puntos e inflexión.
S"i+1(xi+1) = S"i(xi+1) para cada i = 0; 1; . . . ; n − 2;
5) Condiciones de borde: porque faltan datosdn los puntos extremos.
a) S"0(x0) = S"n−1(xn) = S"n(xn) = 0 (frontera libre);
b) S'0(x0) = f'(x0) y S'n−1(xn) = S'n(xn) = f'(xn) (frontera sujeta).
En este caso tomamosla condición de Frontera Sujeta pues se pueden calcular las derivadas en los puntos tomados a partir del trabajo práctico anterior.
Se obtendrá una matriz A con estos datos y condiciones. A será rala, simétrica con la condición de Frontera Sujeta y bien condicionada si los datos del problema son de órdenes similares."""
Justificación de la elección de la Interpolación por "splines"
La interpolación con spline de frontera sujeta es más precisa que la natural, pero requiere los valores de la derivada en los extremos o al menos buenas aproximaciones. Tenemos la ventaja que es posible calcular las derivadas en los puntos tomados en la tabla a partir de la expresión del trabajo práctico anterior.
import matplotlib.pyplot as plt
import numpy as np
xd=np.array([0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1,1.1,1.2,1.3,1.4,1.5])
yd=np.array([-1022.845661,
-30.14890729,
0.1208137489
,3.79802413,
4.594543763,
4.833148038,
4.921240953,
4.958896418,
4.976838489,
4.986134744,
4.991283432,
4.994294217,
4.996136119,
4.997306702,
4.998075298])
def S(x):
x = np.atleast_1d(x) # Convertimos a array
sgn = np.zeros_like(x)
for i in range(len(sgn)):
if x [i]>= 0 : #S0 0,1 y 0,2
return -1022.845661 + 40162.20606605946*(x-0.1) +40162.20606605946 *(x- 0.1 ) + (-650711.8459371064)* (x- 0.1 )**2 + -2169039.4213116993 *(x- 0.1 )**3
elif x [i]<= 0.3: #S1 usa 0,3
return -30.14890729 + (-5561.768437415327) *(x- 0.2 ) + (-5561.768437415327)* (x- 0.2 ) + 97550.20664613906 * (x-0.2)**2 + 2494206.841944152 *(x- 0.2 )**3
elif x[i] <=0.4:
return 0.1208137489 + 1636.5818718766814 *(x- 0.3 ) + 1636.5818718766814* (x- 0.3 ) + (-28217.09044877938) * (x- 0.3 )**2 + -419224.323649728 *(x- 0.3 )**3
elif x[i] <=0.5:
return 3.79802413 + (-422.3760933812574)* (x- 0.4 ) + (-422.3760933812574)* (x- 0.4 ) + 7340.401951639437 * (x- 0.4 )**2 + 118524.97466806276 *(x- 0.4 )**3
elif x [i]<=0.6:
return 4.594543763 + 116.35219996324318 *(x- 0.5 ) + 116.35219996324318* (x- 0.5 ) + (-2008.72458220815)* (x- 0.5 )**2 + (-31163.7551128253)* (x- 0.5 )**3
elif x [i]<= 0.7:
return 4.833148038 + (-30.307984723510824) *(x- 0.6 ) + (-30.307984723510824) *(x- 0.6 ) + 527.1217697940036 * (x- 0.6 )**2 + 8452.821173340515 *(x- 0.6 )**3
elif x [i]<= 0.8:
return 4.921240953 + 8.415295010765005* (x- 0.7 ) + 8.415295010765005 *(x- 0.7 ) + (-144.91590496768234) * (x- 0.7 )**2 + (-2240.125582538951)* (x- 0.7 )**3
elif x [i]<= 0.9:
return 4.958896418 + (-2.1387571513999424) *(x- 0.8 ) + (-2.1387571513999424)* (x- 0.8 ) + 37.410615077414704 *(x- 0.8 )**2 + 607.7550668169903 *(x- 0.8 )**3
elif x [i]<= 1:
return 4.976838489 + 0.6244169682741899 *(x- 0.9 ) + 0.6244169682741899 *(x- 0.9 ) + (-10.640573541831118) * (x- 0.9 )**2 + (-160.17062873081943) *(x- 0.9 )**3
elif x [i]<= 1.1:
return 4.986134744 + (-0.142509539102246)* (x- 1 ) + (-0.142509539102246) *(x- 1 ) + 2.5579342904365348 * (x- 1 )**2 + 43.99502610755881 *(x-1)**3
elif x [i]<=1.2:
return 4.991283432 + 0.051048284514616245* (x- 1.1 ) + 0.051048284514616245* (x- 1.1 ) + (-0.8354337198056924)* (x- 1.1 )**2 + (-11.311226700807438) *(x- 1.1 )**3
elif x [i]<=1.3:
return 4.994294217 + (-0.0066022513180368685)*(x- 1.2 ) + (-0.0066022513180368685) *(x- 1.2 ) + 0.14242968917289545* (x- 1.2 )**2 + 3.25954469659529 *(x- 1.2 )**3
else:
return 4.996136119 + 0.0058365366090152426* (x- 1.3 ) + 0.0058365366090152426* (x- 1.3 ) + (-0.08494993680468522) * (x- 1.3 )**2 + (-0.7579320865919366) *(x- 1.3 )**3
xs=np.linspace(xd[0],xd[-1],5)
print('xs sería igual a:',xs)
xs2=np.linspace(xd[0],xd[-1],100)
print('xs2 sería igual a:',xs2)
plt.figure(figsize=(10,10))
plt.plot(xd,yd,'o')
plt.plot(xs2,S(xs2),'blue')
plt.show()
'''