数值计算-利用Python实现牛顿插值法(不调用函数库)

Xenia ·

更新时间:2024-11-13

· 849 次阅读

数值计算-利用Python实现牛顿插值法（不调用函数库）牛顿插值法介绍实现代码使用实例 牛顿插值法介绍

牛顿插值公式:

已知对于给定节点x0,x1,...,xnx_0,x_1,...,x_nx0,x1,...,xn，对应函数值f(x0),f(x1),...,f(xn)f(x_0),f(x_1),...,f(x_n)f(x0),f(x1),...,f(xn)，

则有牛顿插值公式:

Pn(x)=f(x0)+f[x0,x1]⋅(x−x0)+...+f[x0,x1,...,xn]⋅(x−x0)⋅(x−x1)⋅...⋅(x−xn)P_n(x)=f(x_0)+f[x_0,x_1]\cdot(x-x_0)+...+f[x_0,x_1,...,x_n]\cdot(x-x_0)\cdot(x-x_1)\cdot...\cdot(x-x_n)Pn(x)=f(x0)+f[x0,x1]⋅(x−x0)+...+f[x0,x1,...,xn]⋅(x−x0)⋅(x−x1)⋅...⋅(x−xn)

其中，f[x0,x1,...xn]f[x_0,x_1,...x_n]f[x0,x1,...xn]为nnn阶均差。

实现代码

# newton.py
import matplotlib.pyplot as plt
from pylab import mpl
import numpy as np
def newton_dd(x: list, f: list, n: int):  # 此处 n表示多项式阶数，为节点数-1
    dd = [0 for _ in range(0, n)]  # divided difference
    for i in range(0, n):
        for j in range(n-1, i-1, -1):  # n-1,n-2,...,i 必须从后往前
            if i == 0:
                dd[j] = (f[j + 1] - f[j]) / (x[j + 1] - x[j])
            else:
                dd[j] = (dd[j] - dd[j - 1]) / (x[j + 1] - x[j - i])
    return dd
def get_function(data: np.ndarray, x: list, f: list, n: int):
    parm = newton_dd(x, f, n)
    # print(parm)
    res = f[0]
    w = 1
    for i in range(0, n):
        w *= data - x[i]
        res += parm[i] * w
    return res
def draw(data: np.ndarray, x: list, f: list, n: int):
    plt.plot(data, get_function(data, x, f, n), label="牛顿插值拟合曲线", color="red")
    plt.scatter(x, f, label="节点数据", color="red")
    plt.title("%s次牛顿插值法结果"%n)
    mpl.rcParams['font.sans-serif'] = ['SimHei']
    mpl.rcParams['axes.unicode_minus'] = False
    plt.legend(loc="upper right")

使用实例

考虑龙格函数：
f(x)=1/(1+25x2)f(x)=1/(1+25x^2) f(x)=1/(1+25x2)
在[−1,1][-1,1][−1,1]上取11个等距节点，对这些节点进行牛顿插值。

import numpy as np
import newton as nt # 这是上文的代码文件
import matplotlib.pyplot as plt
def draw(data: np.ndarray, x: list, f:list, n: int):
    nt.draw(data, x, f, n-1)
    x_= list(np.arange(-1, 1.01, 0.01))
    f_= [1 / (1 + 25*pow(num, 2)) for num in x_]
    plt.plot(x_, f_, label="龙格函数", color="black")
    plt.show()
if __name__ == "__main__":
    n = 11
    x = list(np.arange(-1, 1 + 2 / (n - 1), 2 / (n - 1)))
    f = [1 / (1 + 25*pow(num, 2)) for num in x]
    draw(np.arange(-1, 1.01, 0.01), x, f, n)

运行得到结果:
原创文章 2获赞 2访问量 30 关注私信展开阅读全文
作者：IIDEAT