Lagrange、Newton、分段插值法及Python实现

点击上方“机器学习与统计学”,关注并设为“星标” 重磅干货,第一时间送达 数据分析中,经常需要根据已

重磅干货,第一时间送达

打开凤凰新闻,查看更多高清图片

数据分析中,经常需要根据已知的函数点进行数据、模型的处理和分析,而通常情况下现有的数据是极少的,不足以支撑分析的进行,这里就需要使用差值法模拟新的数值来满足需求。

插值法又称“内插法”,是利用函数f(x)在某区间中已知的若干点的函数值,作出适当的特定函数,在区间的其他点上用这特定函数的值作为函数f (x)的近似值。

常用的插值方法有Lagrange插值、Newton插值、分段插值、Hermite插值、样条插值等等。这里我们就介绍一下最常用到的Lagrange、Newton、分段插值法及Python实现。

1、拉格朗日插值法Lagrange插值基本思想是将待求的n次多项式插值函数pn(x)改写成另一种表示方式,再利用插值条件确定其中的待定函数,从而求出插值多项式。它是n次多项式插值,成功地用构造插值基函数的方法解决了求n次多项式插值函数问题。

一般地,若已知 

 在互不相同 n+1 个点 

 处的函数值 

 ( 即该函数过 

 这n+1个点),则可以考虑构造一个过这n+1 个点的、次数不超过n的多项式 

,使其满足:

要估计任一点ξ,ξ≠xi,i=0,1,2,...,n,则可以用Pn(ξ)的值作为准确值f(ξ)的近似值。称式(*)为插值条件(准则),含xi(i=0,1,...,n)的最小区间[a,b],其中a=min{x0,x1,...,xn},b=max{x0,x1,...,xn}。

设集合 

 是关于点 

 的角标的集合, 

 ,作n个多项式 

 。对于任意 

 ,都有 

使得

 是n-1次多项式,且满足 

 并且 

 。最后可得

 

 。

形如上式的插值多项式 

称为拉格朗日(Lagrange)插值多项式。

#coding=utf-8from matplotlib import pyplot as pltdef Lg(data,testdata):    predict=0    data_x=[data[i][0] for i in range(len(data))]    data_y=[data[i][1] for i in range(len(data))]    if testdata in data_x:        #print "testdata is already known"        return data_y[data_x.index(testdata)]    for i in range(len(data_x)):        af=1        for j in range(len(data_x)):            if j!=i:                af*=(1.0*(testdata-data_x[j])/(data_x[i]-data_x[j]))        predict+=data_y[i]*af    return predictdef plot(data,nums):    data_x=[data[i][0] for i in range(len(data))]    data_y=[data[i][1] for i in range(len(data))]    Area=[min(data_x),max(data_x)]    X=[Area[0]+1.0*i*(Area[1]-Area[0])/nums for i in range(nums)]    X[len(X)-1]=Area[1]    Y=[Lg(data,x) for x in X]    plt.plot(X,Y,label=result)    for i in range(len(data_x)):        plt.plot(data_x[i],data_y[i],ro,label="point")    plt.savefig(Lg.jpg)    plt.show()data=[[0,0],[1,2],[2,3],[3,8],[4,2]]print Lg(data,1.5)plot(data,100)

2、牛顿插值Newton插值基本思想是将待求的n次插值多项式Pn(x)改写为具有承袭性的形式,然后利用插值条件⑴确定Pn(x)的待定系数,以求出所要的插值函数。

牛顿差值引入了差商的概念,使其在差值节点增加时便于计算。

设函数 

 ,已知其n+1个插值节点为 

 , 

,我们定义:

 在 

 的零阶差商为 

 在点 

 与 

 的一阶差商为

 

 在点 

 ,

 ,

 的二阶插商为

 

一般的, 

 在点 

 的k 阶差商为

可将k阶差商表示为函数值 

 的组合:

经过分别变形,依次代入,可得牛顿差值公式:

可记为:

其中, 

 为牛顿差值公式的余项或截断误差,当n趋于无穷大时为零。

#coding=utf-8#coding=utf-8from matplotlib import pyplot as pltdef calF(data):    #差商计算  n个数据 0-(n-1)阶个差商 n个数据    data_x=[data[i][0] for i in range(len(data))]    data_y=[data[i][1] for i in range(len(data))]    F= [1 for i in range(len(data))]       FM=[]    for i in range(len(data)):        FME=[]        if i==0:            FME=data_y        else:            for j in range(len(FM[len(FM)-1])-1):                delta=data_x[i+j]-data_x[j]                value=1.0*(FM[len(FM)-1][j+1]-FM[len(FM)-1][j])/delta                FME.append(value)        FM.append(FME)    F=[fme[0] for fme in FM]    print FM    return Fdef NT(data,testdata,F):    #差商之类的计算    predict=0    data_x=[data[i][0] for i in range(len(data))]    data_y=[data[i][1] for i in range(len(data))]    if testdata in data_x:        return data_y[data_x.index(testdata)]    else:        for i in range(len(data_x)):            Eq=1            if i!=0:                for j in range(i):                    Eq=Eq*(testdata-data_x[j])            predict+=(F[i]*Eq)        return predictdef plot(data,nums):    data_x=[data[i][0] for i in range(len(data))]    data_y=[data[i][1] for i in range(len(data))]    Area=[min(data_x),max(data_x)]    X=[Area[0]+1.0*i*(Area[1]-Area[0])/nums for i in range(nums)]    X[len(X)-1]=Area[1]    F=calF(data)    Y=[NT(data,x,F) for x in X]    plt.plot(X,Y,label=result)    for i in range(len(data_x)):        plt.plot(data_x[i],data_y[i],ro,label="point")    plt.savefig(Newton.jpg)    plt.show()data=[[0,0],[1,2],[2,3],[3,8],[4,2]]plot(data,100)

3、分段线性插值

对每一个分段区间(xi,xi+1)分别进行插值,将被插值函数f(x)的插值节点由小到大排序,然后每对相邻的两个节点为端点的区间上用m次多项式去近似f(x)。

将每两个相邻的节点用直线连起来,如此形成的一条折线就是分段线性插值函数。计算x点的插值时,只用到x左右的两个节点,计算量与节点个数n(初始值x0,y0的长度,n=length(x0))无关,而拉格朗日插值与n值有关。分段线性插值中n越大,分段越多,插值误差越小。

#coding=utf-8from matplotlib import pyplot as pltdef DivideLine(data,testdata):    #找到最邻近的    data_x=[data[i][0] for i in range(len(data))]    data_y=[data[i][1] for i in range(len(data))]    if testdata in data_x:        return data_y[data_x.index(testdata)]    else:        index=0        for j in range(len(data_x)):            if data_x[j] 1 ]>testdata:                 index=j                  break         predict= 1.0 *(testdata-data_x[j])*(data_y[j+ 1 ]-data_y[j])/(data_x[j+ 1 ]-data_x[j])+data_y[j]          return  predict def plot(data,nums):     data_x=[data[i][ 0 ]  for  i  in  range(len(data))]     data_y=[data[i][ 1 ]  for  i  in  range(len(data))]     Area=[min(data_x),max(data_x)]     X=[Area[ 0 ]+ 1.0 *i*(Area[ 1 ]-Area[ 0 ])/nums  for  i  in  range(nums)]     X[len(X) -1 ]=Area[ 1 ]     Y=[DivideLine(data,x)  for  x  in  X]     plt.plot(X,Y,label= result )      for  i  in  range(len(data_x)):         plt.plot(data_x[i],data_y[i], ro ,label= "point" )     plt.savefig( DivLine.jpg )     plt.show() data=[[ 0 , 0 ],[ 1 , 2 ],[ 2 , 3 ],[ 3 , 8 ],[ 4 , 2 ]] print DivideLine(data, 1.5 ) plot(data, 100 )