数值分析之拉格朗日插值、牛顿插值、分段线性插值实现

最新推荐文章于 2023-11-27 20:51:00 发布

BUAA-XX

最新推荐文章于 2023-11-27 20:51:00 发布

阅读量8.5k

点赞数 2

分类专栏：数值分析文章标签：编程算法导论

本文链接：https://blog.csdn.net/sinat_33829806/article/details/78387843

版权

数值分析专栏收录该内容

7 篇文章

订阅专栏

本文介绍了三种常用的插值方法：拉格朗日插值法、牛顿插值法及分段线性插值法，并通过实例展示了每种方法的计算过程及特点。拉格朗日插值法适用于全局信息的插值，而牛顿插值法在数据更新时更为高效。分段线性插值法则因其简单性和实用性，在工业界应用广泛。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

1、拉格朗日插值法

考虑全局信息的比较经典的插值方法，编程简单，计算量大。

#coding=utf-8
from matplotlib import pyplot as plt

def 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 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=[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、牛顿插值

相比于拉格朗日的优势在于牛顿插值法在重新引入数据时可以仅仅更新F矩阵进行。（手动维护差商表的话，效率提高很多，但对计算机而言，觉得差不多）。复杂度和拉格朗日一样，O(n2)。

#coding=utf-8
from matplotlib import pyplot as plt

def 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 F

def 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 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]
    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、分段线性插值

最简单，工业使用最多。

#coding=utf-8
from matplotlib import pyplot as plt
def 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]<testdata and  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)


![结果](http://img.blog.csdn.net/20171029205027016?watermark/2/text/aHR0cDovL2Jsb2cuY3Nkbi5uZXQvc2luYXRfMzM4Mjk4MDY=/font/5a6L5L2T/fontsize/400/fill/I0JBQkFCMA==/dissolve/70/gravity/SouthEast)