牛顿均差插值多项式

1.均差定义

称$f[x_0,x_k]=\frac{f(x_k)-f(x_0)}{x_k-x_0}为函数f(x)关于点x_0,x_k的\text{一阶均差}$
$f[x_0,x_1,x_k]=\frac{f[x_0,x_k]-f[x_0,x_1]}{x_k-x_1}称为f(x)的\text{二阶均差.}$
一般地，称
$$f[x_0,x_1,\dots ,x_k]=\frac{f[x_0,\dots ,x_{k-2},x_k]-f[x_0,x_1,\dots ,x_{k-1}]}{x_k-x_{k-1}}$$
为$f(x)的\text{k阶均差}(均差也称为\text{差商.})$

2.均差表

3.牛顿插值多项式

$一般情形已知f在插值点x_i(i=0,1,\dots ,n)$上的值为$f(x_i)(i=0,1,\dots ,n),要求n次插值多项式P_n(x)满足条件$
$$P_n(x_i)=f(x_i),i=0,1,\dots ,n,$$
则$P_n(x)$可表示为
$$P_n(x)=f(x_0)+f[x_0,x_1](x-x_0)+f[x_0,x_1,x_2](x-x_0)(x-x_1)+\dots +f[x_0,x_1,\dots ,x_n](x-x_0)\dots (x-x_{n-1})$$
与拉格朗日插值不同，这里$P_n(x)$是由基函数{1,x−x0,…,(x−x0)…(x−xn−1)}\{1,x-x_0,\ldots,(x-x_0)\ldots(x-x_{n-1})\}{1,x−x0​,…,(x−x0​)…(x−xn−1​)}逐次递推得到.

4.Python实现牛顿插值多项式均差系数

import numpy as np
from sympy import expand
from sympy.abc import x


# 自己原创
def difference_quotient_coefficient(x_array: np.ndarray, y_array: np.ndarray):
    """
    实现牛顿均差插值多项式各阶均差系数
    :param x_array: 横坐标数组
    :param y_array: 纵坐标数组
    :return: 各阶均差系数数组
    """
    nth = x_array.shape[0]
    # 初始化均差系数数组
    difference_quotient_array = np.zeros(nth)
    difference_quotient_array[0] = y_array[0]
    # 求各阶均差
    for i in range(1, nth):
        nth_difference_quotient = np.zeros(nth - i)
        # 求每阶均差
        for j in range(nth - i):
            nth_difference_quotient[j] = (y_array[j + 1] - y_array[j]) / (x_array[i + j] - x_array[j])
        difference_quotient_array[i] = nth_difference_quotient[0]
        y_array = nth_difference_quotient
    return difference_quotient_array

5.Python实现牛顿插值多项式基函数

# 自己原创
def newtown_base_function(symbol_x, x_array: np.ndarray):
    """
    实现牛顿均差插值基函数
    :param symbol_x: 符号变量x
    :param x_array: 横坐标数组
    :return: 牛顿均差插值基函数数组
    """
    nth = x_array.shape[0]
    base_function_array = np.zeros(nth, dtype=np.ndarray)
    # base_function_array = list(np.zeros(nth))
    base_function_array[0] = 1
    for i in range(1, nth):
        base_function_array[i] = np.prod(symbol_x - x_array[:i], dtype=np.ndarray)
    return base_function_array

6.Python实现牛顿插值多项式

# 自己原创
def newtown(x_array: np.ndarray, y_array: np.ndarray, symbol_x):
    """
    实现牛顿均差插值多项式
    :param x_array: 插值结点横坐标数组
    :param y_array: 插值结点纵坐标数组
    :param symbol_x: 符号变量x
    :return: 牛顿均差插值多项式
    """
    return expand(np.sum(difference_quotient_coefficient(x_array, y_array) * newtown_base_function(symbol_x, x_array)))

7.秦九韶算法

# 自己原创
def qin_jiu_shao(coefficient_array: np.ndarray, value, is_ascending_order: bool = False):
    """
    多项式求值的秦九韶算法,即Horner算法
    :param coefficient_array: 多项式系数数组
    :param value: 自变量的值
    :param is_ascending_order: 多项式排列顺序
    :return: the value of polynomial
    example:
    2*x**4 - 3*x**2 + 3*x - 4=((((2x+0)x-3)x+3)x-4)
    3*x**5 - 2*x**3 + x + 7=(((((3x+0)x-2)x+0)x+1)x+7)
    """
    if is_ascending_order:
        coefficient_array = coefficient_array[::-1]
    degree = coefficient_array.shape[0]
    b = coefficient_array[0]
    for i in range(1, degree):
        b = b * value + coefficient_array[i]
    return b

8.测试

if __name__ == '__main__':
	# 牛顿均差插值多项式测试成功，来源详见李庆扬数值分析第5版P48，e.x.1
	x_coefficient_list = [1, -1, 2]
    x_coefficient_array = np.array(x_coefficient_list)
    y = np.array([0, -3, 4])
    # 牛顿均差插值多项式均差系数测试成功
    print("牛顿均差插值多项式均差系数:{}".format(difference_quotient_coefficient(x_coefficient_array, y)))
    # 牛顿均差插值多项式基函数测试成功
    print("牛顿均差插值多项式基函数:{}".format(newtown_base_function(x, x_coefficient_array)))
    newtown_polynomial = newtown(x_coefficient_array, y, x)
    # 计算牛顿插值多项式(降幂排列)在x=2处的值
    print("牛顿插值多项式(降幂排列)在x=2处的值为:{}".format(newtown_polynomial.subs(x, 2)))
    # 创建牛顿插值多项式对象
    newtown_interpolation_polynomial = Poly(newtown_polynomial, x)
    # 获得牛顿插值多项式降幂排列系数
    newtown_polynomial_coefficient = np.array(newtown_interpolation_polynomial.coeffs())  # list 类型
    # 使用秦九韶算法计算牛顿插值多项式(降幂排列)在x=2处的值
    print("秦九韶算法计算牛顿插值多项式(降幂排列)在x=2处的值为:{}".format(qin_jiu_shao(newtown_polynomial_coefficient, 2)))
    # 牛顿插值多项式次数
    print("牛顿插值多项式次数为:{}".format(newtown_interpolation_polynomial.degree()))
    # 牛顿插值多项式子项次数
    print("牛顿插值多项式子项次数为:{}".format(newtown_interpolation_polynomial.monoms()))

9.运行结果