Python数值计算（22）——akima插值-优快云博客

本文链接：https://blog.youkuaiyun.com/cdinten/article/details/140934667

前面介绍到的很多插值算法，在多项式阶次较高时，都会存在龙格现象，最直观的表现就是导数的变化很急剧，所以一般采用的方式是使用分段的三次插值，例如前面介绍到的三次样条曲线，这里再介绍Akima插值算法，通过该方法生成的插值曲线，具有更加平滑的特性。

1. 一点数学知识

该算法最初由Hiroshi Akima在1970年提出，论文名称是《A new method of interpolation and smooth curve fitting based on local procedures》，传送门在这里。

Aikma插值法规定在两个实测点之间进行内插,除需要用到这两个实测值外，还要用这两个点相近邻的四个实测点上的观测值，也就是说在两个实测点之间进行内插共需六个实测点。Akima方法和其他拟合方法类似，假设需要拟合的一组数据点是 $(x_i,y_i)$ ，其中 $i=0,1,2,...,n$ ，需要确定的就是 $y=f(x)$ 的多项式，具有连续的一阶导数，在任意相邻数据点之间使用阶次最高为3的多项式。

与样条曲线类似，可用方程如下：

$\left\{\begin{matrix} y_i=f(x_i)\\ y_{i+1}=f(x_{i+1})\\ y'_i=t_i\\ y'_{i+1}=t_{i+1} \end{matrix}\right.$

上述公式可以唯一的确定一个三次多项式，问题的关键是如何得到后面两个方程中的 $t_i,t_{i+1}$ 。

如果只取六个点 $(x_i,y_i)$ ，其中 $i=1,2,3,4,5,6$ ，插值点位于第3,4个点之间，即 $x_3<x<x_4$ ，则插值表达式如下（公式0）：

$y=p_0+p_1(x-x_3)+p_2(x-x_3)^2+p_3(x-x_3)^3$

可得：

$\left\{\begin{matrix} p_0=y_3\\ p_1=t_3\\ p_2=\frac{3\frac{y_4-y_3}{x_4-x_3}-2t_3-t_4}{x_4-x_3}\\ p_3=\frac{t_3+t_4-2\frac{y_4-y_3}{x_4-x_3}}{(x_4-x_3)^2} \end{matrix}\right.$

其中按如下公式1和2给出：

$t_i=\frac{|m_{i+1}-m_i|m_{i-1}+|m_{i-1}-m_{i-2}|m_i}{|m_{i+1}-m_i|+|m_{i-1}-m_{i-2}|},i=3,4$

$m_i=\frac{y_{i+1}-y_i}{x_{i+1}-x_i}$

这就可以解决该点处的斜率，但是随之而来的问题是，起始的两个点 $(i=1,2)$ ，和最后的两个点 $(i=5,6)$ ，它们的斜率该如何计算呢？它们的值通过如下公式3计算：

$t_2=2*t_3-t_4\\ t_1=2*t_2-t_3\\ t_5=2*t_4-t_3\\t_6=2*t_5-t_4$

推广一下，所以对于 $(x_i,y_i)$ ， $i=0,1,2,...,n$ 共计 $(n+1)$ 个点，可以先通过公式1和2计算出 $t_2,t_3,...,t_{n-2}$ ，然后再通过公式3计算 $t_0,t_1,t_{n-1},t_n$ 。

另外一个问题是，如果公式1中分母为零时如何处理？即有 $m_{i+1}==m_i \cap m_{i-1}==m_{i-2}$ 时的情况，此时规定 $t_i=(m_{i-1}+m_i)/2$

2. 算法实现

在搞清楚了其算法原理以后，我们可以很容易写出算法，这里采用和前面类似的方式，构造一个类，实现对akima插值算法的管理以及对值的计算，代码如下：

from scipy.interpolate import PPoly
from scipy.interpolate import Akima1DInterpolator
import numpy as np
from numpy.polynomial import Polynomial
import matplotlib.pyplot as plt

class akimaIntp:
    __coef:np.ndarray
    __bpx:np.ndarray
    __bpy:np.ndarray
    def __init__(self,x:np.ndarray,y:np.ndarray):
        n,=x.shape        
        m=np.zeros(n+3)
        m[2:-2]=np.diff(y)/np.diff(x)
        m[1]=2*m[2]-m[3]
        m[0]=2*m[1]-m[2]
        m[-2]=2*m[-3]-m[-4]
        m[-1]=2*m[-2]-m[-3]
        d=np.zeros(n)
        for i in range(n):
            w1=np.abs(m[i+3]-m[i+2])
            w2=np.abs(m[i+1]-m[i])
            if w1==0 and w2==0:
                d[i]=(m[i+1]+m[i+2])/2
            else:
                d[i]=w1/(w1+w2)*m[i+1]+w2/(w1+w2)*m[i+2]

        self.__coef=np.zeros((n-1,4))
        self.__coef[:,0]=y[0:-1]
        self.__coef[:,1]=d[0:-1]
        for i in range(n-1):
            self.__coef[i,2]=(3*(y[i+1]-y[i])/(x[i+1]-x[i])-2*d[i]-d[i+1])/(x[i+1]-x[i])
            self.__coef[i,3]=(d[i+1]+d[i]-2*(y[i+1]-y[i])/(x[i+1]-x[i]))/(x[i+1]-x[i])/(x[i+1]-x[i])

        self.__bpx=x.copy()
        self.__bpy=y.copy()


    def __call__(self,w:np.ndarray):
        n,=w.shape
        ret= np.zeros(n)
        for i in range(n):
            if self.__bpx[0]>=w[i]:
                ret[i]=self.__bpy[0]
                continue
            if self.__bpx[-1]<=w[i]:
                ret[i]=self.__bpy[-1]
                continue
            j=0
            while self.__bpx[j]<w[i]:
                j+=1
            cp=self.__coef[j-1,:]
            p=Polynomial(0)
            for t in range(len(cp)):
                p+=cp[t]*Polynomial([-self.__bpx[j-1],1])**t
            ret[i]=p(w[i])
        return ret

    @property
    def c(self):
        return self.__coef
    @property
    def x(self):
        return self.__bpx

3. 算法测试

我们以HIROSHI AKIMA当初论文中的一组数据为例，测试如下：

np.set_printoptions(precision=4,linewidth=100)
# test data
x=np.array([0,1,2,3,4,5,6,7,8,9,10])
y=np.array([10,10,10,10,10,10,10.5,15,50,60,85])

aip=akimaIntp(x,y)
print(aip.c)
'''
[[ 10.       0.       0.       0.    ]
 [ 10.       0.       0.       0.    ]
 [ 10.       0.       0.       0.    ]
 [ 10.       0.       0.       0.    ]
 [ 10.       0.       0.       0.    ]
 [ 10.       0.       0.9355  -0.4355]
 [ 10.5      0.5645   3.6641   0.2714]
 [ 15.       8.7069  69.3444 -43.0513]
 [ 50.      18.2418 -25.8585  17.6168]
 [ 60.      19.375    3.75     1.875 ]]
'''

其中每一行就是公式0中每个分段中的系数 $p_0,p_1,p_2,p_3$ 。

4. 现有算法

Scipy中也有现成的工具包供使用，scipy.interpolate中的Akima1DInterpolator类，可完成该项工作，使用起来非常方便：

akima=Akima1DInterpolator(x,y)
print(akima.c)
'''
[[  0.       0.       0.       0.       0.      -0.4355   0.2714 -43.0513  17.6168   1.875 ]
 [  0.       0.       0.       0.       0.       0.9355   3.6641  69.3444 -25.8585   3.75  ]
 [  0.       0.       0.       0.       0.       0.       0.5645   8.7069  18.2418  19.375 ]
 [ 10.      10.      10.      10.      10.      10.      10.5     15.      50.      60.    ]]
'''

注意，系数和前面的很多工具包类似，它是按列升幂排列的，即每一列表示每个分段区间的系数，而每列的第一个系数是最高次幂。

对比两者拟合的图形，将其绘制在同一个图表中：

X=np.linspace(x[0],x[-1],100)
Y1=aip(X)
plt.plot(X,Y1,'b-.')
Y2=akima(X)
plt.plot(X,Y2,'r--')
plt.grid()
plt.show()

效果如下：

其中蓝色点画线是我们自己实现的算法，红色虚线是软件包提供的算法，可见两者吻合的非常好。