import numpy as np
import matplotlib.pyplot as plt

# 构造数据
T = np.array([0, 20, 40, 60, 80])
sigma = np.array([78.2, 73.4, 70.2, 66.7, 63.2])

# 计算差商


def divided_difference(x, y):
    n = len(x)
    F = np.zeros([n, n])
    F[:, 0] = y.astype(float)
    for j in range(1, n):
        for i in range(n-j):
            F[i][j] = (F[i+1][j-1] - F[i][j-1])/(x[i+j]-x[i])
    return F[0]

# 求解插值多项式


def newton_interpolation(x, y, X):
    a = divided_difference(x, y)
    n = len(x) - 1
    p = a[n]
    for k in range(1, n+1):
        p = a[n-k] + (X - x[n-k])*p
    return p


# 生成插值点
T_interp = np.linspace(0, 80, 100)
sigma_interp = newton_interpolation(T, sigma, T_interp)

# 绘制图形
plt.plot(T, sigma, 'o', label='data')
plt.plot(T_interp, sigma_interp, label='Newton Interpolation')
plt.xlabel('T(℃)')
plt.ylabel(r'\sigma \times 10 ^ 3 (N/m)')
plt.legend()
plt.ylim(0, 100)
plt.savefig('sigma_vs_T.png')

在以上代码中补充代码使结果输出利用牛顿插值法求出的穿过上述数据点的多项式