2

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

'''
0

1 respuesta 1

1

Una lista de problemas:

  • El test if x [i]>= 0 siempre sale verdadero. ¿Quizás cambiarlo por if x [i]<= 0.2?
  • El for i in range... que ya hace un return para la primera i que encuentra. Como x es un array de numpy, hace el cálculo de la primera expresión con todo el array x.
  • En vez de a0 + a1(x-0.1) + a2(x-0.1)**2 + a3(x-0.1)**3 se repite el termino a1(x-0.1) en cada una de las expresiones.

Aún cambiando todo eso, no consigo el gráfico deseado. El siguiente paso es dibujar todos los splines por separado, y solapando un poco.

Haciendo un zoom a la parte interesante: zoom de los diferentes splines

Se nota que la curva para 0.2 funciona entre 0.1 y 0.2, y así las otras. Pues, la versión corregida de los splines tiene que tener en cuenta estos límites.

El código modificado para los 2 plots:

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_k(x, k):
    x = np.atleast_1d(x)
    if k == 1:  # S0 0,1 y 0,2
        return -1022.845661 + 40162.20606605946 * (x - 0.1) + (-650711.8459371064) * (x - 0.1) ** 2 + -2169039.4213116993 * (x - 0.1) ** 3
    elif k == 2:
        return -30.14890729 + (-5561.768437415327) * (x - 0.2) + 97550.20664613906 * (x - 0.2) ** 2 + 2494206.841944152 * (x - 0.2) ** 3
    elif k == 3:
        return 0.1208137489 + 1636.5818718766814 * (x - 0.3) + (-28217.09044877938) * (x - 0.3) ** 2 + -419224.323649728 * (x - 0.3) ** 3
    elif k == 4:
        return 3.79802413 + (-422.3760933812574) * (x - 0.4) + 7340.401951639437 * (x - 0.4) ** 2 + 118524.97466806276 * (x - 0.4) ** 3
    elif k == 5:
        return 4.594543763 + 116.35219996324318 * (x - 0.5) + (-2008.72458220815) * (x - 0.5) ** 2 + (-31163.7551128253) * (x - 0.5) ** 3
    elif k == 6:
        return 4.833148038 + (-30.307984723510824) * (x - 0.6) + 527.1217697940036 * (x - 0.6) ** 2 + 8452.821173340515 * (x - 0.6) ** 3
    elif k == 7:
        return 4.921240953 + 8.415295010765005 * (x - 0.7) + (-144.91590496768234) * (x - 0.7) ** 2 + (-2240.125582538951) * (x - 0.7) ** 3
    elif k == 8:
        return 4.958896418 + (-2.1387571513999424) * (x - 0.8) + 37.410615077414704 * (x - 0.8) ** 2 + 607.7550668169903 * (x - 0.8) ** 3
    elif k == 9:
        return 4.976838489 + 0.6244169682741899 * (x - 0.9) + (-10.640573541831118) * (x - 0.9) ** 2 + (-160.17062873081943) * (x - 0.9) ** 3
    elif k == 10:
        return 4.986134744 + (-0.142509539102246) * (x - 1) + 2.5579342904365348 * (x - 1) ** 2 + 43.99502610755881 * (x - 1) ** 3
    elif k == 11:
        return 4.991283432 + 0.051048284514616245 * (x - 1.1) + (-0.8354337198056924) * (x - 1.1) ** 2 + (-11.311226700807438) * (x - 1.1) ** 3
    elif k == 12:
        return 4.994294217 + (-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.08494993680468522) * (x - 1.3) ** 2 + (-0.7579320865919366) * (x - 1.3) ** 3


def S(x):
    x = np.atleast_1d(x)
    y = np.zeros_like(x)
    for i in range(len(x)):
        if x[i] <= 0.1:
            y[i] = -1022.845661 + 40162.20606605946 * (x[i] - 0.1) + (-650711.8459371064) * (
                        x[i] - 0.1) ** 2 + -2169039.4213116993 * (x[i] - 0.1) ** 3
        elif x[i] <= 0.2:
            y[i] = -30.14890729 + (-5561.768437415327) * (x[i] - 0.2) + 97550.20664613906 * (
                        x[i] - 0.2) ** 2 + 2494206.841944152 * (x[i] - 0.2) ** 3
        elif x[i] <= 0.3:
            y[i] = 0.1208137489 + 1636.5818718766814 * (x[i] - 0.3) + (-28217.09044877938) * (
                        x[i] - 0.3) ** 2 + -419224.323649728 * (x[i] - 0.3) ** 3
        elif x[i] <= 0.4:
            y[i] = 3.79802413 + (-422.3760933812574) * (x[i] - 0.4) + 7340.401951639437 * (
                        x[i] - 0.4) ** 2 + 118524.97466806276 * (x[i] - 0.4) ** 3
        elif x[i] <= 0.5:
            y[i] = 4.594543763 + 116.35219996324318 * (x[i] - 0.5) + (-2008.72458220815) * (x[i] - 0.5) ** 2 + (
                -31163.7551128253) * (x[i] - 0.5) ** 3
        elif x[i] <= 0.6:
            y[i] = 4.833148038 + (-30.307984723510824) * (x[i] - 0.6) + 527.1217697940036 * (
                        x[i] - 0.6) ** 2 + 8452.821173340515 * (x[i] - 0.6) ** 3
        elif x[i] <= 0.7:
            y[i] = 4.921240953 + 8.415295010765005 * (x[i] - 0.7) + (-144.91590496768234) * (x[i] - 0.7) ** 2 + (
                -2240.125582538951) * (x[i] - 0.7) ** 3
        elif x[i] <= 0.8:
            y[i] = 4.958896418 + (-2.1387571513999424) * (x[i] - 0.8) + 37.410615077414704 * (
                        x[i] - 0.8) ** 2 + 607.7550668169903 * (x[i] - 0.8) ** 3
        elif x[i] <= 0.9:
            y[i] = 4.976838489 + 0.6244169682741899 * (x[i] - 0.9) + (-10.640573541831118) * (x[i] - 0.9) ** 2 + (
                -160.17062873081943) * (x[i] - 0.9) ** 3
        elif x[i] <= 1.0:
            y[i] = 4.986134744 + (-0.142509539102246) * (x[i] - 1) + 2.5579342904365348 * (
                        x[i] - 1) ** 2 + 43.99502610755881 * (x[i] - 1) ** 3
        elif x[i] <= 1.1:
            y[i] = 4.991283432 + 0.051048284514616245 * (x[i] - 1.1) + (-0.8354337198056924) * (x[i] - 1.1) ** 2 + (
                -11.311226700807438) * (x[i] - 1.1) ** 3
        elif x[i] <= 1.2:
            y[i] = 4.994294217 + (-0.0066022513180368685) * (x[i] - 1.2) + 0.14242968917289545 * (
                        x[i] - 1.2) ** 2 + 3.25954469659529 * (x[i] - 1.2) ** 3
        else:
            y[i] = 4.996136119 + 0.0058365366090152426 * (x[i] - 1.3) + (-0.08494993680468522) * (x[i] - 1.3) ** 2 + (
                -0.7579320865919366) * (x[i] - 1.3) ** 3
    return y


plt.figure(figsize=(10, 10))
plt.plot(xd, yd, 'o')
for k in range(1, 14):
    x = np.linspace(k * 0.1 - 0.12, k * 0.1 + 0.12, 50)
    plt.plot(x, S_k(x, k), label=f'{k * 0.1:.1f}')
plt.legend()
plt.show()

plt.figure(figsize=(10, 10))
xs = np.linspace(xd[0], xd[-1], 200)
plt.plot(xd, yd, 'o', color='tomato')
plt.plot(xs, S(xs), 'blue')
plt.show()

Así se queda el resultado. el resultado final

1
  • Muchas gracias!!
    – SilverSong
    Commented el 18 dic. 2019 a las 15:35

Tu Respuesta

By clicking “Publica tu respuesta”, you agree to our terms of service and acknowledge you have read our privacy policy.

¿No es la respuesta que buscas? Examina otras preguntas con la etiqueta o formula tu propia pregunta.