Python数值计算（23）——modified akima插值

1. 数学原理

在前面的Akima插值中，计算斜率使用如下公式：

如果记：

在出现分母分子同时为零的情况时，会出现NaN的计算结果，Akima他自己也意识到这种问题，因此，在原来的算法上做了修订，也就是Modified Akima算法，有时也被简写为makima算法。修订如下：

这样的话，就解决了在出现超过2个点具有连续值时的无法被零除问题，此时。

2. 算法实现

修改原来的算法如下：

import numpy as np
from numpy.polynomial import Polynomial
import matplotlib.pyplot as plt

class makimaIntp:
    __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])+np.abs(m[i+3]+m[i+2])/2
        w2=np.abs(m[i+1]-m[i])+np.abs(m[i+1]+m[i])/2
        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=makimaIntp(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.9412  -0.4412]
 [ 10.5      0.5588   4.2111  -0.2699]
 [ 15.       8.1713  68.8387 -42.01  ]
 [ 50.      19.8187 -27.1375  17.3187]
 [ 60.      17.5      9.8684  -2.3684]]
'''

同样的，每一行就是每个分段中的系数。

4. 现成算法

Scipy中也有现成的工具包供使用，scipy.interpolate中的Akima1DInterpolator类，可完成该项工作，只需要指定参数method='makima'即可。

akima=Akima1DInterpolator(x,y,method='makima')
print(akima.c)
'''
[[  0.    0.    0.    0.    0.    -0.4412   -0.2699     -42.01      17.3187     -2.3684] 
 [  0.    0.    0.    0.    0.    0.9412    4.2111      68.8387     -27.1375    9.8684] 
 [  0.    0.    0.    0.    0.    0.        0.5588      8.1713      19.8187     17.5   ] 
 [ 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()

效果如下：

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