python样条插值（二）

python实现

样条插值在MATLB中有自带的库，在python中也有，位于scipy库中，具体定义如下：

scipy官方文档：

class scipy.interpolate.interp1d(x, y, kind='linear', axis=-1, copy=True, bounds_error=None, fill_value=nan, assume_sorted=False)

 上为定义，下为参数解释：

• interp1d ，scipy.interpolate包里有很多的模块可以实现对一些已知的点进行插值，即找到一个合适的函数，例如模块 interp1d。当得到插值函数后便可用这个插值函数计算其他xj对应的的yj值了,这也就是插值的意义所在。

from scipy.interpolate import interp1d
import numpy as np

noise = np.random.normal(0, 0.1, 100)
x = np.linspace(0, 10, 100)
y = np.sin(x) + noise
f = interp1d(x, y)
f2 = interp1d(x, y, kind='cubic')
print x[10], np.sin(x[10]), y[10], f(x[10]), f2(x[10])
print x[11], np.sin(x[11]), y[11], f(x[11]), f2(x[11])
xm = (x[10] + x[11]) / 2
print xm, np.sin(xm), (y[10] + y[11]) / 2, f(xm), f2(xm)
print f
xnew = np.linspace(0, 10, 40)
import matplotlib.pyplot as plt
plt.plot(x,y,'o',xnew,f(xnew),'-', xnew, f2(xnew),'--', xnew, np.sin(xnew),linewidth=2)
plt.legend(['data', 'linear', 'cubic', "$cos(x)$"], loc='best')
plt.show()

上为参考文档的代码片，其中介绍了如何使用该库的插值函数，运行结果如下：

该函数可以直接输入一个一维数组，并返回一位数组所对应的一维数组结果。按照此特性，我对原题的点进行了拟合尝试：

x=np.array([0,45,75,105,135,165,225,255])
y=np.array([0,20,60,60,20,-60,-100,20])
xi=np.arange(0,max(x),0.1)
xc=x[:8]
yc=y[:8]
f = interp1d(xc,yc, kind = 2)
ff = interp1d(x,y,kind = 3)
plt.plot(f(xi),xi,'g',label='2 derivative')
plt.plot(ff(xi),xi,'m',label='3 derivative')

plt.xlim(-135,75)
plt.ylim(0,260)
plt.plot(y,x,'kx')
plt.xlabel('y')
plt.ylabel('x')
plt.title('Spline interpolation')
plt.legend()
plt.show()

下图为所给的点集： 

​所给点的坐标显示

下图为拟合结果： 

拟合后的图像

 

绿色的是二次样条插值，紫色的是三次样条插值。

很明显，在这两个插值函数中，开头这个部分有明显的区别，并且从曲率的角度看，二次样条插值明显优于三次样条插值。

但是由于自带库的样条插值函数并没有给我们输入参数的机会，查阅资料，在计算机插值拟合时，数值计算软件如Matlab都把非扭结边界条件作为默认的边界条件。

所以理论上来说，我们在给定特定的点后无法将边界条件作为输入条件进行更改。

那么，能不能通过加入适当的点使得曲率减小呢？

通过加入新的系列点集，使得插值函数的曲率减小，能不能采取梯度下降的方法寻找新点的最佳坐标？

我经过尝试，发现当我加入新的点之后，会导致其他部分的插值函数也发生变化（必然。。），所以加入点的做法并不太适合，需要考虑的问题变得更加复杂了。

import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import interp1d


x=np.array([0,45,75,105,135,165,225,255,202.7])
y=np.array([0,20,60,60,20,-60,-100,20,-106.2])
xi=np.arange(0,max(x),0.1)
xc=x[:8]
yc=y[:8]
f = interp1d(xc,yc, kind = 3)
ff = interp1d(x,y,kind = 3)
for i in range(len(x)-1):
    x_i=np.linspace(x[i],x[i+1],4)
    # x_=np.arange(x[i],x[i+1],0.1)
    y_i=f(x_i)
    fi=np.polyfit(x_i,y_i,3)
    fi_np=np.poly1d(fi)
    y_i_new=ff(x_i)
    fi_new=np.polyfit(x_i,y_i_new,3)
    fi_new_np=np.poly1d(fi_new)
    # fi_yi=fi_np(x_)
    # plt.plot(fi_yi,x_,'b')
    print('第',i+1,'段表达式')
    print(fi_np,'fx')
    print(fi_new_np,'fx_new')
# print(f._kind)
plt.plot(f(xi),xi,'g',label='Original derivative')
plt.plot(ff(xi),xi,'m',label='New derivative')

plt.xlim(-135,75)
plt.ylim(0,260)
plt.plot(y,x,'kx')
plt.xlabel('y')
plt.ylabel('x')
plt.title('Spline interpolation')
plt.legend()
plt.show()

加入了一个新点（202.7，-106.2），效果图如下：

发现其他部分的拟合差别对比并不明显，但是当我们分析各段的解析式：

我们会发现每个解析式都会发生变化，所以相应的曲率也会发生改变，由此来看，插入点的方法并不适用。

所以我决定自己来写样条插值的程序。