数值分析——Newton插值的Python实现

python代码

import matplotlib.pyplot as plt
from pylab import mpl
import numpy as np
import math
x0=0.596
x=[0.4,0.55,0.65,0.80,0.90,1.05]
y=[0.41075,0.57815,0.69675,0.88811,1.02652,1.25382]
n=len(x)
num_df=5 #4次插值，所以要计算5次差商的值
A=np.zeros((n,n))
A[:,0]=np.transpose(y)

def difference_quotient(x, y):
    A=np.zeros((n,n))
    A[:,0]=np.transpose(y)
    k=0
    j=0
    i=n
    for j in range(1,n):
        for i in range(j,n):
            A[i,j]=(A[i,j-1]-A[i-1,j-1])/(x[i]-x[i-1-k]) #求差商的值
        k+=1
    print(A)
    return A
def function(data):
    y0=0
    for k in range(1,num_df):
        d=A[k,k]
        for i in range(0,k):
            d*=(data-x[i])
        y0+=d
    y0+=x[0]
    return y0
def error(data):  
    d=A[n-1,n-1]
    for i in range(0,n):
        d*=(data-x[i])
    print(d)
    return d
def calculate_data(x,A):
    returnData=[]
    for data in x:
        returnData.append(function(data))
    return returnData
def draw(newData):
    plt.scatter(x,y,label="离散数据",color="red")
    plt.plot(x,newData,label="牛顿插值拟合曲线",color="black")
    plt.scatter(0.596,function(0.596),label="预测函数点",color="blue")
    plt.title("牛顿插值法")
    mpl.rcParams['font.sans-serif'] = ['SimHei']
    mpl.rcParams['axes.unicode_minus'] = False
    plt.legend(loc="upper left")
    plt.text(0.596,function(0.596), (0.596,function(0.596)))
    plt.show()
A=difference_quotient(x,y)
yy=calculate_data(x,A)
draw(yy)
error(0.596)

结果如图