插值方法—数值分析实验

目的与要求：
（一）目的通过设计、编制、调试2~3个多项式插值、拟合曲线的程序，加深对其数值计算方法及有关的基础理论知识的理解。
（二）要求 用编程语言实现拉格朗日（Lagrange）插值多项式、牛顿（Newton）插值、用线性函数拟合给定数据的程序。
示例
1、问题已知插值节点序列，用Lagrange插值多项式计算的函数在点的近似值。
2、算法描述（略）
3、程序中变量说明 （略）
4、源程序清单及运行结果（略）
5、按以上4点要求编写上机实验报告。
做法：
查看老师的matlab实验代码，给出了几个插值节点(0.4,-0.9163) (0.5,-0.6931) (0.6,-0.5108) (0.7,-0.3567)

1. 拉格朗日插值公式
2. 牛顿插值公式：
   ①构造差商表
   
   ②插值多项式代码：

import numpy as np
import matplotlib.pyplot as plt
def main(x):
    # 模拟 p(x)=2x+1 给定
    b,a = 2,1
    # 给出4个插值节点
    xi = np.array([0.4,0.5,0.6,0.7])
    yi = np.array([-0.9163,-0.6931,-0.5108,-0.3567])
    yLangage = 0
    # 拉格朗日插值节点值为yLangage
    for k in range(4):
        t = yi[k]
        for j in range(4):
            if j!=k:
                t*=(x-xi[j])/(xi[k]-xi[j])
        yLangage += t
        
    # 牛顿插值节点值为yNewton
    # [[0,0,0,0],[0,0,0,0],[0,0,0,0],[0,0,0,0]]
    # 首先构造差商表
    yNewton = 0
    Q = [[0] * 4 for _ in range(4)]
    for n in range(4):
        Q[n][0]=yi[n]
    for i in range(1,4):
        for j in range(i,4):
            Q[j][i]=(Q[j][i-1]-Q[j-1][i-1])/(xi[j]-xi[j-i])
    # 得到差商表对角线值
    qDiagonal=[]
    for i in range(4):
        qDiagonal.append(Q[i][i])
    yNewton += qDiagonal[0]
    # 代入②插值多项式
    for i in range(1,4):
        p = qDiagonal[i]
        for j in range(i):
            p*=(x-xi[j])
        yNewton+=p
        
    # 线性拟合
    z = np.polyfit(xi, yi, 1)
    p = np.poly1d(z)
    y = []
    for i in range(4):
        y.append(p(xi[i]))
    # 拉格朗日插值节点值：  yLangage
    # 牛顿插值节点值：      yNewton
    # 通过matplotlib模块进行画图
    # 图分为三块 序号1位置为拉格朗日插值图 序号2位置为牛顿插值图
    # 序号3为线性拟合图
    plt.subplot(221)
    plt.plot(xi,yi,'k')
    plt.scatter(x,yLangage,color='red')
    plt.title('Lagrange Example')
    
    plt.subplot(222)
    plt.plot(xi,yi,'k')
    plt.scatter(x, yNewton, color='red')
    plt.title('Newton Example')
    
    plt.subplot(223)
    plt.plot(xi,y,'k')
    plt.scatter(xi,yi,color='red')
    plt.title('Linar Example')
    plt.savefig('DifValue test',dpi=600)
    plt.show()
    print("拉格朗日插值结果为：",yLangage)
    print("牛顿插值结果为：",yNewton)
if __name__ == '__main__':
    x = 0.54
    main(x)