多项式插值拟合（二）

多项式二阶插值公式：y=ax^2+bx+c

interp1d(x,y,kind='quadratic')中kind参数值为'quadratic'

import numpy as np
from scipy.interpolate import interp1d,interp2d
import pylab as pl

x=np.array([1,2,3,4,5])
y=np.array([1,4,2,4,1])

fc=interp1d(x,y,kind='quadratic')
xint=np.linspace(x.min(),x.max(),20)
print("x:",xint)
print("y:",fc(xint))
pl.scatter(xint,fc(xint),marker='o',c='r')
pl.plot(xint,fc(xint),color='blue')
pl.show()

插值点：

x: [1.         1.21052632 1.42105263 1.63157895 1.84210526 2.05263158
 2.26315789 2.47368421 2.68421053 2.89473684 3.10526316 3.31578947
 3.52631579 3.73684211 3.94736842 4.15789474 4.36842105 4.57894737
 4.78947368 5.        ]
y: [1.         2.1966759  3.09196676 3.68587258 3.97839335 3.96952909
 3.65927978 3.04764543 2.37894737 2.04210526 2.04210526 2.37894737
 3.04764543 3.65927978 3.96952909 3.97839335 3.68587258 3.09196676
 2.1966759  1.        ]

自定义函数实现：

3点确定方程参数

a=-((y2-y3)*x1-(x2-x3)*y1+x2*y3-x3*y2)/((x2-x3)*(x1-x2)*(x1-x3))
b = ((y2 - y3) * x1*x1 +x2*x2*y3-x3*x3*y2-(x2**2-x3**2)*y1) / ((x2 - x3) * (x1 - x2) * (x1 - x3))
c = ((x2*y3-x3*y2)*x1*x1  - (x2*x2*y3-x3*x3*y2) * x1 +(x2*x2*x3-x2*x3*x3)*y1) / ((x2 - x3) * (x1 - x2) * (x1 - x3))

import numpy as np
from scipy.interpolate import interp1d,interp2d
import pylab as pl


x=np.array([1,2,3,4,5])
y=np.array([1,4,2,4,1])


xint=np.linspace(x.min(),x.max(),20)

# 2、二次
def quadratic_inter(x1,y1,x2,y2,x3,y3,nx):
    a=-((y2-y3)*x1-(x2-x3)*y1+x2*y3-x3*y2)/((x2-x3)*(x1-x2)*(x1-x3))
    b = ((y2 - y3) * x1*x1 +x2*x2*y3-x3*x3*y2-(x2**2-x3**2)*y1) / ((x2 - x3) * (x1 - x2) * (x1 - x3))
    c = ((x2*y3-x3*y2)*x1*x1  - (x2*x2*y3-x3*x3*y2) * x1 +(x2*x2*x3-x2*x3*x3)*y1) / ((x2 - x3) * (x1 - x2) * (x1 - x3))

    ny=a*nx*nx+b*nx+c

    return ny

#
for i in range(0,len(x),2):
    if i + 2 < len(x):
        print("---------")
        # print(ny)
        nx=[]
        for val in xint:
            if  val>=x[i] and val<=x[i+2]:
                nx.append(val)
    
        print(nx)
        nx=np.array(nx)
        ny = quadratic_inter(x[i], y[i], x[i+1], y[i+1],x[i+2],y[i+2], nx)
        print(ny)

        pl.scatter(nx,ny,c='b',marker='o')
        pl.plot(nx,ny,color='red')

pl.show()

分段连接存在不连贯，可优化。

插值点：

nx: [1.0, 1.2105263157894737, 1.4210526315789473, 1.631578947368421, 1.8421052631578947, 2.052631578947368, 2.263157894736842, 2.473684210526316, 2.6842105263157894, 2.894736842105263]
ny: [1.         2.04709141 2.87257618 3.47645429 3.85872576 4.01939058
 3.95844875 3.67590028 3.17174515 2.44598338]
--------
nx: [3.1052631578947367, 3.3157894736842106, 3.526315789473684, 3.7368421052631575, 3.9473684210526314, 4.157894736842105, 4.368421052631579, 4.578947368421052, 4.789473684210526, 5.0]
ny: [2.44598338 3.17174515 3.67590028 3.95844875 4.01939058 3.85872576
 3.47645429 2.87257618 2.04709141 1.    ]