Python数值计算（12）——线性插值_业界新闻

发布时间:2024-08-03 03:28

阅读量:0

1. 概述

插值是根据已知的数据序列（可以理解为你坐标中一系列离散的点），找到其中的规律，然后根据找到的这个规律，来对其中尚未有数据记录的点进行数值估计的方法。最简单直观的一种插值方式是线性插值，它是数学、计算机图形学等领域广泛使用的一种简单插值方法，根据已知数据点来估某一点的函数值。线性插值法的特点是简单、方便，适用于函数在某一段区间内是线性的情况，即函数在该区间内可以用一条直线来近似表示。线性插值法的优点是计算量小，缺点是对于。多个点而言在连接点出出在明显的间断点。

2. 一点数学知识

线性插值是简单直观的，假设已知坐标与，要得到区间内某一位置在直线上的值，根据点斜式直线方程，我们只需要求出斜率即可，显然有：

$k=\frac{y_1-y_0}{x_1-x_0}$

则：

$y-y_0=\frac{y_1-y_0}{x_1-x_0}*(x-x_0)$

如果写成斜截式 $y=kx+b$ ：

$y-y_0=\frac{y_1-y_0}{x_1-x_0}*x-\frac{y_1x_0-y_0x_1}{x_1-x_0}$

当然如果时，就违背了函数单一性原则（不讨论广义上的曲线）。

对于多个点，可以在每个区间段内采用分段线性插值。

3. 算法实现

3.1 单段线性插值

单段线性插值就是确定一个线性函数，然后通过函数估算其他点处的函数值。

以下定义个函数，返回一个多维数组，表示多项式的系数，并通过该返回值创建多项式，随后计算处的值，并在图形中表示：

import numpy as np from numpy.polynomial import Polynomial import matplotlib.pyplot as plt  def linear_2p(x:np.ndarray,y:np.ndarray):     k=(y[1]-y[0])/(x[1]-x[0])     b=y[0]-k*x[0]     return np.array([b,k])  x=np.array([2,5]) y=np.array([3,7])  coef= linear_2p(x,y) # 线性插值函数的系数 p=Polynomial(coef) # 构造多项式 print(p) # 0.33333333 + 1.33333333 x print(p(2.5)) # 3.6666666666666665  # 绘制图形 X=np.linspace(x[0],x[1],100) Y=p(X) plt.plot(X,Y,'r')  #plt.plot(x,y,’r’)  #注:这两种在线性时绘图效果相同， # 但是实际含义不同，前者是用多个点绘制拟合曲线， # 后者仅用线段连接起止点 plt.plot(2.5,p(2.5),'b*') plt.grid() plt.show()

绘制图形如下：

如果需要计算单个点或者多个点的插值，其实不需要去计算该段的直线，也可以使用定比分点的概念：

$P=(1-t)P_1+tP_2,0\leqslant t \leqslant 1$

因此，插值可以使用如下函数：

def linear_2p_intp(x:np.ndarray,y:np.ndarray,w):     t=(w-x[0])/(x[1]-x[0]) return (1-t)*y[0]+t*y[1]

例如，实现5个点的插值并显示：

x=np.array([2,5]) y=np.array([3,7]) w=np.linspace(3,4,10) plt.plot(x,y,'r') # 简化绘图 yw=linear_2p_intp(x,y,w) plt.plot(w,yw,'b*') plt.grid() plt.show()

运行效果如下：

动态演示效果如下：

3.2 多段线性插值

多段线性插值在每一段上都是两点线性插值，假设点序列(X,Y)中X为单调递增，即具有如下特点：

$x_i<x_j (0 \leqslant i < j \leqslant n )$

估算x=w点处的值时，首先需要定位w属于哪个区间段， $x_i \leqslant w<x_{i+1}$

为此，在构造该算法时，除了需要各段的函数外，还要有各段的区间信息，定义类如下：

from scipy.interpolate import PPoly import numpy as np from numpy.polynomial import Polynomial as P import matplotlib.pyplot as plt  class multiSegLinearIntp:     __coef:np.ndarray # coefficient     __bps:np.ndarray # breakpoints     __seg=0 # segments count     def __init__(self,x:np.ndarray,y:np.ndarray):         n,=x.shape         k=np.diff(y)/np.diff(x)         self.__coef=np.zeros((2,n-1))         self.__coef[0,:]=y[0:-1]-k*x[0:-1]         self.__coef[1,:]=k         self.__bps=x.copy()         self.__seg=n-1          def __call__(self,x:np.ndarray):         n,=x.shape         y=np.zeros(n)         for i in range(n):             w=x[i]             if w<self.__bps[0]:                 y[i]=self.__coef[0,0]+self.__coef[1,0]*w                 continue             if w>=self.__bps[-1]:                 y[i]=self.__coef[0,-1]+self.__coef[1,-1]*w                 continue             j=0             while w>=self.__bps[j]:                 j+=1             y[i]=self.__coef[0,j-1]+self.__coef[1,j-1]*w         return y          @property     def c(self):         return self.__coef     @property     def seg(self):         return self.__seg

该类multiSegLinearIntp的构造函数将numpy.ndarray的实例x和y作为参数，内部保存各分段区间的系数，实例对象在收到传入值x后，查找x其所在的分段，并在该段返回函数的值，测试如下：

x=np.array([0,1,2,3,4,5,6,7,8,9]) y=np.array([2,4,3,3,1,5,6,3,1,0])  z=multiSegLinearIntp(x,y) x1=np.linspace(-2,11,200) y1=z(x1) plt.plot(x,y,'r') w=np.array([-2,11]) print(z(w)) #[-2. -2.] #plt.plot(x1,y1,'b-') #y2=np.interp(x1,x,y) #plt.plot(x1,y2,'g*') plt.grid() plt.show()  x=np.array([0,1,2,3,4,5,6,7,8,9]) y=np.array([2,4,3,3,1,5,6,3,1,0])  z=multiSegLinearIntp(x,y) x1=np.linspace(-2,11,200) y1=z(x1) plt.plot(x,y,'r') w=np.array([-2,4.5,11]) yw=z(w) print(z(w)) #[-2.  3. -2.] plt.plot(w,yw,'b*') plt.grid() plt.show()

运行效果如下：

可以看到，在点处计算的值为，这是相符的，对于不在区间内的点和，可以看到依旧保持了这种线性关系，这就已经是外插了（extrapolation），是否要应用这种关系，需要根据实际情况判断。

动态演示效果如下：

4. 现有工具包

通过前面一节的例子，发现自己实现的多段线性插值还是挺麻烦的，有现场的工具包吗？当然。numpy中，使用numpy. Interp函数实现插值运算，其函数原型为：

numpy.interp(x, xp, fp, left=None, right=None, period=None)

其中x是待估算值的横坐标，xp，fp是已知点序列，left和right是x没有落在插值区间时的值，缺省值是left=x[0]，right=x[-1]，period是周期型，通常用于角坐标的插值，返回值是与x具有同样长度的多维数组。

测试如下：

x=np.array([2,5]) y=np.array([3,7]) w=np.linspace(3,4,5) print(linear_2p_intp(x,y,w))  print(np.interp(w,x,y))  ''' 输出均为[4.33333333 4.66666667 5 5.33333333 5.66666667] '''

另外，在scipy软件包中，scipy.interpolate.interp1d类也可以实现线性插值，但是该类已经作为遗留代码，不在被更新，在以后得升级中可能会被移除，官方给的建议是使用前面提到的numpy.interp函数。以下是简单的一个示例，供参考：

x=np.array([2,5]) y=np.array([3,7]) w=np.linspace(3,4,5) f=interp1d(x,y) plt.plot(w,f(w),'r*-') plt.grid() plt.show()

运行效果为：

5. 双线性插值

双线性插值（Bilinear interpolation）是有两个变量的插值函数的线性插值扩展，其核心思想是在两个方向分别进行一次线性插值。假如我们想得到未知函数在点的值，假设我们已知函数 f 在 $Q_{11} = (x_1, y_1),Q_{12} = (x_1, y_2),Q_{21} = (x_2, y_1) ,Q_{22} = (x_2, y_2)$ 四个点的值。首先在方向进行线性插值，然后在方向进行线性插值。这种插值方法并不是线性的，而是两个线性函数的乘积。线性插值的结果与插值的顺序无关。首先进行方向的插值，然后进行方向的插值，所得到的结果是一样的。

基本数学原理如下：

先在方向插值：

$f(R_1)=\frac{x_2-x}{x_2-x_1}f(Q_{11})+\frac{x-x_1}{x_2-x_1}f(Q_{21}),R_1=(x,y_1)\\ f(R_2)=\frac{x_2-x}{x_2-x_1}f(Q_{12})+\frac{x-x_1}{x_2-x_1}f(Q_{22}),R_2=(x,y_2)$