数据插值：Lagrange插值方法

数据插值：Lagrange插值方法

数学理论基础

设$y=f(x)$是区间$[a,b]$上的实值函数，$x_0,x_1,\cdots ,x_n$是$[a,b]$上的$n+1$个互异节点，$l_i(x)=\frac{{\prod }_{j=0,j\ne i}^n(x-x_j)}{{\prod }_{j=0,j\ne i}^n(x_i-x_j)}$$(i=0,1,\cdots ,n)$是拉格朗日插值基函数。

$$f(x)=\sum _{i=0}^n{y}_i{l}_i(x)+\frac{f^{(n+1)}(\xi )}{(n+1)!}\prod _{i=0}^n(x-x_i)$$

代码模板

1.插值预测值

import numpy as np


def lagrange_interpolation(x_nodes, y_nodes, x):
    """
    实现拉格朗日插值方法
    :param x_nodes: 插值节点的x坐标
    :param y_nodes: 插值节点的y坐标
    :param x: 要计算插值的点的x坐标
    :return: 在点x处的插值结果
    """
    n = len(x_nodes)
    result = 0
    for i in range(n):
        l_i = 1
        for j in range(n):
            if j != i:
                l_i *= (x - x_nodes[j]) / (x_nodes[i] - x_nodes[j])
        result += y_nodes[i] * l_i
    return result

2.插值预测函数表达式

import sympy as sp

def lagrange_polynomial(x_nodes, y_nodes):
    # 定义符号变量 x
    x = sp.symbols('x')
    n = len(x_nodes)
    # 初始化拉格朗日插值多项式为 0
    polynomial = 0
    for i in range(n):
        # 初始化拉格朗日插值基函数为 1
        l_i = 1
        for j in range(n):
            if j != i:
                # 计算拉格朗日插值基函数的每一项
                l_i *= (x - x_nodes[j]) / (x_nodes[i] - x_nodes[j])
        # 累加每一项 y_i * l_i(x) 到多项式中
        polynomial += y_nodes[i] * l_i
    # 展开多项式
    polynomial = sp.expand(polynomial)
    return polynomial

可视化插值函数与数据关系

import sympy as sp
import numpy as np
import matplotlib.pyplot as plt


def lagrange_polynomial(x_nodes, y_nodes):
    # 定义符号变量 x
    x = sp.symbols('x')
    n = len(x_nodes)
    # 初始化拉格朗日插值多项式为 0
    polynomial = 0
    for i in range(n):
        # 初始化拉格朗日插值基函数为 1
        l_i = 1
        for j in range(n):
            if j != i:
                # 计算拉格朗日插值基函数的每一项
                l_i *= (x - x_nodes[j]) / (x_nodes[i] - x_nodes[j])
        # 累加每一项 y_i * l_i(x) 到多项式中
        polynomial += y_nodes[i] * l_i
    # 展开多项式
    polynomial = sp.expand(polynomial)
    return polynomial


# 示例数据
x_nodes = [1, 2, 3, 4]
y_nodes = [2, 4, 1, 5]

# 计算拉格朗日插值多项式的解析式
poly = lagrange_polynomial(x_nodes, y_nodes)
print("拉格朗日插值多项式的解析式为:", poly)

# 将符号表达式转换为可调用的函数
x = sp.symbols('x')
poly_func = sp.lambdify(x, poly, 'numpy')

# 生成用于绘制插值函数的 x 值
x_vals = np.linspace(min(x_nodes) - 1, max(x_nodes) + 1, 400)
y_vals = poly_func(x_vals)

# 绘制原始数据点
plt.scatter(x_nodes, y_nodes, color='red', label='Original Data')

# 绘制插值函数曲线
plt.plot(x_vals, y_vals, color='blue', label='Lagrange Interpolation')

# 设置图表标题和坐标轴标签
plt.title('Lagrange Interpolation Visualization')
plt.xlabel('x')
plt.ylabel('y')

# 显示图例
plt.legend()

# 显示网格线
plt.grid(True)

# 显示图形
plt.show()

不计算基插函数情况下的Neville算法

import numpy as np
import matplotlib.pyplot as plt


def neville_interpolation_optimized(x_nodes, y_nodes, x):
    n = len(x_nodes)
    # 只使用一维数组来存储中间结果
    P = y_nodes.copy()
    for j in range(1, n):
        for i in range(n - j):
            numerator = (x - x_nodes[i + j]) * P[i] - (x - x_nodes[i]) * P[i + 1]
            denominator = x_nodes[i] - x_nodes[i + j]
            P[i] = numerator / denominator
    return P[0]


# 示例数据
x_nodes = [1, 2, 3, 4]
y_nodes = [2, 4, 1, 5]

# 生成用于绘制插值函数的 x 值
x_vals = np.linspace(min(x_nodes) - 1, max(x_nodes) + 1, 400)
y_vals = []
for x in x_vals:
    y = neville_interpolation_optimized(x_nodes, y_nodes, x)
    y_vals.append(y)

# 绘制原始数据点
plt.scatter(x_nodes, y_nodes, color='red', label='Original Data')

# 绘制插值函数曲线
plt.plot(x_vals, y_vals, color='blue', label='Neville Interpolation')

# 设置图表标题和坐标轴标签
plt.title('Neville Interpolation Visualization')
plt.xlabel('x')
plt.ylabel('y')

# 显示图例
plt.legend()

# 显示网格线
plt.grid(True)

# 显示图形
plt.show()