Как построить производную кубического сплайна в Python 3?

Question

Я использую Python 3 для задачи, связанной с численным анализом.

Я должен построить несколько точек, полученных из функции синуса. Более того, мне нужно сделать кубическую интерполяцию этих точек (кубический сплайн).

Итак, эти задачи выполнены. Выходная картинка отличная и код работает. Однако мне нужно проверить, не является ли производная кубического сплайна функцией косинуса.

Посмотрите на это изображение:

В оранжевом вы видите функцию косинуса. Синим цветом вы видите функцию синуса. В красном, 5 точек, которые я выбрал. В фиолетовом вы видите линейную интерполяцию. А в тире вы видите кубическую интерполяцию.

Мне нужно построить производную пунктирной кривой и сравнить ее с оранжевой.

Интуитивно, я знаю , они будут очень похожи. Однако я не смог доказать , что с помощью графа.

Это код:

import math
import random
from numpy import array 
import numpy as np 
import matplotlib.pyplot as plot
from scipy.interpolate import interp1d
from scipy import interpolate
from scipy.misc import derivative as deriv

def random_sine():

    lista_f_x = []
    lista_x = []

    for i in range(1,6):

        aleatorio = random.uniform(0,360)
        aleatorio = math.radians(aleatorio)
        lista_x.append(aleatorio)

        sine_random = math.sin(aleatorio)
        lista_f_x.append(sine_random)

    lista_x = array(lista_x)
    lista_f_x = array(lista_f_x)

    return ("x",lista_x,"f(x)", lista_f_x) 

# para ter um grupo controle melhor, deixe esses valores aleatórios, gerados uma vez, como fixos
fixed_x = array([5.80990031, 1.7836885,  4.62073799, 0.89337425, 5.62219906])
fixed_y = array([-0.45581264,  0.97742392, -0.99580299,  0.77919112, -0.61389568])

x = fixed_x
y = fixed_y

"""
caso deseje usar os valores fixos
basta inserir o comentário "#" nas linhas
39, 40 e 41 abaixo
"""
#teste_dinamico = random_sine()
#x = teste_dinamico[1]
#y = teste_dinamico[3]

time = np.arange(0,10,0.1)

amplitude = np.sin(time)

amplitude_cosine = np.cos(time)

plot.plot(time, amplitude, time, amplitude_cosine)

plot.title('Função Seno')

plot.xlabel('Coordenadas de X')

plot.ylabel('Seno(x)')

plot.grid(True, which='both')

plot.axhline(y=0, color='k')

pares_x_y = list(zip(x,y))

sort_pares_x_y = sorted(pares_x_y)

x_ordenado = []
y_ordenado_simetric = []

for i in sort_pares_x_y:

    x_ordenado.append(i[0])
    y_ordenado_simetric.append(i[1])

x_ordenado = array(x_ordenado)
y_ordenado_simetric = array(y_ordenado_simetric)

f = interp1d(x_ordenado, y_ordenado_simetric)

f2 = interp1d(x_ordenado, y_ordenado_simetric, kind="cubic")

plot.plot(x_ordenado, f2(x_ordenado))

minimo = min(x_ordenado)
maximo = max(x_ordenado)

xnew = np.linspace(minimo, maximo, num=400, endpoint=True)

plot.plot(x_ordenado, y_ordenado_simetric, 'o', xnew, f(xnew), '-', xnew, f2(xnew), '--')

plot.scatter(x,y)

plot.show()

Я пытался подписаться на эту запись . Тем не менее, пост предлагает другой пакет. Я использую interp1d до сих пор, и пост предлагает использовать interpolate. Я пытался преобразовать то, что я сделал до сих пор, но это не сработало.

Что я могу сделать?

Если я строю кривую с помощью matplotlib, есть ли более простой способ построить производную кривой?

Какова лучшая стратегия для решения этой проблемы?

1 - Должен ли я повторить попытку и изменить все на другой пакет, предложенный для SO?

2 - Должен ли я использовать какой-либо численный метод для дифференциации?

Заранее спасибо.

Answer 1

Другие пакеты или функции могут, конечно, быть более удобными, но для понимания вы можете просто вычислить производную (dy / dx) и построить ее.

plot.plot(xnew[:-1], np.diff(f2(xnew))/np.diff(xnew), color="red")

Добавление этой строки приведет к

где красная линия является производной от f2.

Answer 2

Я бы предпочел использовать interpolate.splev вместо interp1d, поскольку, помимо обеспечения интерполяции, первое также позволяет легко вычислять производные с помощью простого аргумента der=n, где n - это порядок производной. Я должен был сделать только незначительных изменений в вашем коде, чтобы все заработало.

Ниже я показываю только соответствующие строки кода, которые я добавил / изменил (выделил комментарием), и полученный график прилагается. Красная пунктирная линия - искомая производная, сравнимая с функцией косинуса.

f2 = interpolate.splrep(x_ordenado, y_ordenado_simetric) # Added

minimo = min(x_ordenado)
maximo = max(x_ordenado)

xnew = np.linspace(minimo, maximo, num=400, endpoint=True)
ynew = interpolate.splev(xnew, f2) # Added 
ynew_der = interpolate.splev(xnew, f2, der=1) # Added to compute first derivative

plot.plot(x_ordenado, y_ordenado_simetric, 'o', xnew, f(xnew), '-', xnew, ynew, '--') # Modified
plot.plot(xnew, ynew_der, '--r') # Derivative Added

выход