Approximation polynomial fit to given data

Question

I'm facing a problem with aproximation polynomial fit. A more detailed review of my problem is shown HERE. Basically I want to smooth the mid section of the polynomial that I got through long mathematical manipulation and is shown below. The idea that I got is to take a range from the beginning and a range from the end of the polynomial (its points) and make an approximation polynomial that hasn't got the waving in the mid section; selected points are shown below. However by using the polyfit function I didn't quite get the polynomial that I've hoped for.

This the crucial part of my code:

from sympy import*
import numpy as np
init_printing()
%matplotlib inline
import matplotlib.pyplot as plt

S = np.linspace(25, 400, 1000)
SS = np.log10(S)
def f(logS): # the polynomial that I got
    return 10**(-57.2476193068601*logS**5 + 585.900632193053*logS**4 - 2384.35277925916*logS**3 + 4821.25582425353*logS**2 - 4845.47249368281*logS + 1943.75303313331)  

xdata = f(SS)
ydata = S

plt.figure(figsize=(12, 10), dpi= 600, facecolor='w', edgecolor='k')

plt.loglog(xdata, ydata,'k-')
# Data for new polynomial (approximation)
XX = xdata[:40].tolist()
X1 = xdata[len(xdata)-350:].tolist()
YY = ydata[:40].tolist()
Y1 = ydata[len(ydata)-350:].tolist()
for i in range(len(X1)):
    XX.append(X1[i])
    YY.append(Y1[i])

# Approximation polynomial fit
plt.loglog(XX,YY,'mo')
x_t = np.linspace(np.min(xdata), np.max(xdata),1000)
p = np.poly1d(np.polyfit(XX,YY, deg=5))

plt.loglog(x_t, p(x_t),'g-')
print(p)

plt.ylim((np.min(ydata), np.max(ydata)))
plt.xlim((np.min(xdata), np.max(xdata)*0.7))  

plt.xlabel('xdata')
plt.ylabel('ydata')
plt.grid(True, which="both")

I would like some advice on how to fix this approximation polynomial problem or if there is another way to solve the waving in the mid section - a better way of smoothing. Any help is highly appreciated!

Answer 1

You are trying to fit a polynomial on regular data which is in exp-exp scale, and only plot it in log-log, where it looks like a polynomial. You will not be able to represent such relation with a polynomial. Preprocess everything to be in log scale in the first place, fit polynomial there, and if you want to go back to your exp-exp world, post-process again.

from sympy import*
import numpy as np
import matplotlib.pyplot as plt

S = np.linspace(25, 400, 1000)
SS = np.log10(S)
def f(logS): # the polynomial that I got
    return 10**(-57.2476193068601*logS**5 + 585.900632193053*logS**4 - 2384.35277925916*logS**3 + 4821.25582425353*logS**2 - 4845.47249368281*logS + 1943.75303313331)  



xdata = np.log(f(SS))
ydata = np.log(S)

plt.figure(figsize=(12, 10),facecolor='w', edgecolor='k')
plt.plot(xdata, ydata, 'k-')
#plt.loglog(xdata, ydata,'k-')
# Data for new polynomial (approximation)
XX = xdata[:40].tolist()
X1 = xdata[len(xdata)-350:].tolist()
YY = ydata[:40].tolist()
Y1 = ydata[len(ydata)-350:].tolist()
for i in range(len(X1)):
    XX.append(X1[i])
    YY.append(Y1[i])


# Approximation polynomial fit
#plt.loglog(XX,YY,'mo')
x_t = np.linspace(np.min(xdata), np.max(xdata),1000)
p = np.poly1d(np.polyfit(XX,YY, deg=5))

#plt.loglog(x_t, p(x_t),'g-')
plt.plot(x_t, p(x_t),'g-')
print(p)

plt.ylim((np.min(ydata), np.max(ydata)))
plt.xlim((np.min(xdata), np.max(xdata)))  

plt.xlabel('xdata')
plt.ylabel('ydata')
plt.grid(True, which="both")
plt.show()