Получите уравнение из полифита, чтобы применить дальнейшие производные - PullRequest
Новая
       24

Получите уравнение из полифита, чтобы применить дальнейшие производные

0 голосов
/ 06 ноября 2019

Мне удалось выполнить полиномиальное сопоставление с набором сглаженных данных. Теперь мне нужно получить полиномиальное уравнение, чтобы я мог выполнить первую и вторую производные для полинома и найти максимумы и минимумы. Хотя я могу легко получить коэффициенты с помощью polyfit, я не знаю, как получить весь многочлен, чтобы я мог применить к нему производные. Я думал об использовании SimPy для выполнения производных, но я не думаю, что смогу без построения полного полиномиального уравнения.

Вот мой код: 

import matplotlib.pyplot as plt
import numpy as np
from pandas import *
from math import factorial
from scipy import interpolate
import numpy.polynomial.polynomial as poly
import itertools

# savitzky_golay definition
def savitzky_golay(y, window_size, order, deriv=0, rate=1):
    try:
        window_size = np.abs(np.int(window_size))
        order = np.abs(np.int(order))
    except ValueError:
        raise ValueError("window_size and order have to be of type int")
    if window_size % 2 != 1 or window_size < 1:
        raise TypeError("window_size size must be a positive odd number")
    if window_size < order + 2:
        raise TypeError("window_size is too small for the polynomials order")
    order_range = range(order+1)
    half_window = (window_size -1) // 2
    # precompute coefficients
    b = np.mat([[k**i for i in order_range] for k in range(-half_window, half_window+1)])
    m = np.linalg.pinv(b).A[deriv] * rate**deriv * factorial(deriv)
    # pad the signal at the extremes with
    # values taken from the signal itself
    firstvals = y[0] - np.abs( y[1:half_window+1][::-1] - y[0] )
    lastvals = y[-1] + np.abs(y[-half_window-1:-1][::-1] - y[-1])
    y = np.concatenate((firstvals, y, lastvals))
    return np.convolve( m[::-1], y, mode='valid')


import os, sys
#get the xlsx files
files= os.listdir(r'C:\Users\Dani\Desktop\REACTIVITY')
file=[f for f in files if f[-4:]=='xlsx']

for file in file:
    #read the files
    data=pandas.read_excel(file, sheet_name=0)

    y1=np.array(data["Viscosity"])
    x=np.array(data["Time (s)"])
    y2=np.array(data["Temperature (oC)"])

    fig, ax1 = plt.subplots()
    #create initial axis
    ax1.set_xlabel("Time (s)", labelpad=10)
    ax1.set_ylabel("Viscosity", color='blue', labelpad=10)
    ax1.tick_params(axis='y', labelcolor='blue')

    #smooth data by using savitzky_golay
    yhat=savitzky_golay(y1, 51, 3)
    ax1.plot(x,yhat, color='blue',linewidth=2)
    #create temperature graph
    ax2=ax1.twinx()
    ax2.set_ylabel(u'Temperature (C)' , color='red', labelpad=10)
    ax2.plot(x, y2, color='red',linewidth=2)
    ax2.tick_params(axis='y', labelcolor='red')

    #fit the Viscosity smoothed data to a polynomial
    for i in itertools.count():
        #create the polynomial fit for different degrees
        coefs=poly.polyfit(x, yhat, i)
        fit=poly.polyval(x, coefs)
        #calculate rsquare
        ybar = np.sum(yhat)/len(yhat)
        ssreg = np.sum((fit-ybar)**2)
        sstot = np.sum((yhat - ybar)**2)
        square = ssreg / sstot

        if square > 0.995:
            #get the desired degree
            ax1.plot(x, fit, color='black', linestyle='--',linewidth=2)
            print(square)
            print(coefs)
            print(i)
            break

    plt.show()
...