1. 数学原理
在前面的Akima插值中，计算斜率使用如下公式：
如果记：
在出现分母分子同时为零的情况时，会出现NaN的计算结果，Akima他自己也意识到这种问题，因此，在原来的算法上做了修订，也就是Modified Akima算法，有时也被简写为makima算法。修订如下：
这样的话，就解决了在出现超过2个点具有连续值时的无法被零除问题，此时。
2. 算法实现
修改原来的算法如下：
import numpy as np
from numpy.polynomial import Polynomial
import matplotlib.pyplot as pltclass makimaIntp:__coef:np.ndarray__bpx:np.ndarray__bpy:np.ndarraydef __init__(self,x:np.ndarray,y:np.ndarray):n,=x.shapem=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])/2w2=np.abs(m[i+1]-m[i])+np.abs(m[i+1]+m[i])/2if w1==0 and w2==0:d[i]=(m[i+1]+m[i+2])/2else: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.shaperet= np.zeros(n)for i in range(n):if self.__bpx[0]>=w[i]:ret[i]=self.__bpy[0]continueif self.__bpx[-1]<=w[i]:ret[i]=self.__bpy[-1]continuej=0while self.__bpx[j]<w[i]:j+=1cp=self.__coef[j-1,:]p=Polynomial(0)for t in range(len(cp)):p+=cp[t]*Polynomial([-self.__bpx[j-1],1])**tret[i]=p(w[i])return ret@propertydef c(self):return self.__coef@propertydef x(self):return self.__bpx3. 算法测试
依旧以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()效果如下：
其中蓝色点画线是我们自己实现的算法，红色虚线是软件包提供的算法，可见两者吻合的非常好。
