数值计算方法之插值方法的Python实现

本系列将不加解释地给出各自数值计算方法的Python实现，旨在方便开卷考时的快速查找。本文是插值方法部分。

1.线性插值

import numpy as np

def linear_interpolation(x_points, y_points, x_target):
    """
    线性插值算法。

    参数：
        x_points (list or np.ndarray): 已知点的x坐标列表
        y_points (list or np.ndarray): 已知点的y坐标列表
        x_target (float or np.ndarray): 需要插值的目标x坐标

    返回：
        float or np.ndarray: 对应的插值结果
    """
    if len(x_points) != len(y_points):
        raise ValueError("x_points和y_points的长度必须相等！")

    # 转为numpy数组，方便操作
    x_points = np.array(x_points)
    y_points = np.array(y_points)

    # 如果目标x是单个值
    if np.isscalar(x_target):
        if x_target < x_points[0] or x_target > x_points[-1]:
            raise ValueError("目标x在已知点范围之外！")
        for i in range(len(x_points) - 1):
            if x_points[i] <= x_target <= x_points[i + 1]:
                # 线性插值公式
                return y_points[i] + (y_points[i + 1] - y_points[i]) * (x_target - x_points[i]) / (x_points[i + 1] - x_points[i])

    # 如果目标x是多个值
    else:
        x_target = np.array(x_target)
        results = []
        for xt in x_target:
            if xt < x_points[0] or xt > x_points[-1]:
                results.append(None)  # 如果超出范围，返回None
            else:
                for i in range(len(x_points) - 1):
                    if x_points[i] <= xt <= x_points[i + 1]:
                        # 线性插值公式
                        results.append(y_points[i] + (y_points[i + 1] - y_points[i]) * (xt - x_points[i]) / (x_points[i + 1] - x_points[i]))
                        break
        return np.array(results)

# 测试代码
if __name__ == "__main__":
    x = [1, 2, 3, 4, 5]
    y = [2, 4, 6, 8, 10]
    target_x = 2.5

    # 单点插值
    interpolated_value = linear_interpolation(x, y, target_x)
    print(f"插值结果: x = {target_x}, y = {interpolated_value}")

    # 多点插值
    target_x_array = [1.5, 3.5, 4.5]
    interpolated_values = linear_interpolation(x, y, target_x_array)
    print(f"插值结果: x = {target_x_array}, y = {interpolated_values}")


    #调用scipy包进行插值对比

    from scipy import interpolate
    f = interpolate.interp1d(x, y, kind='linear')
    print(f(2.5))
    print(f([1.5, 3.5, 4.5]))

2.Nevile插值

def nevilles_interpolation(x_points, y_points, x_target):
    """
    使用Neville插值法进行插值。

    参数：
        x_points (list or np.ndarray): 已知点的x坐标列表
        y_points (list or np.ndarray): 已知点的y坐标列表
        x_target (float): 需要插值的目标x坐标

    返回：
        float: 插值结果
    """
    if len(x_points) != len(y_points):
        raise ValueError("x_points和y_points的长度必须相等！")

    n = len(x_points)
    Q = [[0 for _ in range(n)] for _ in range(n)]

    # 初始化Q表的第一列
    for i in range(n):
        Q[i][0] = y_points[i]

    # 填充Neville表
    for j in range(1, n):
        for i in range(n - j):
            Q[i][j] = ((x_target - x_points[i + j]) * Q[i][j - 1] - (x_target - x_points[i]) * Q[i + 1][j - 1]) / (x_points[i] - x_points[i + j])

    return Q[0][n - 1]

# 测试代码
if __name__ == "__main__":
    x = [1, 2, 3, 4, 5]
    y = [2, 4, 6, 8, 10]
    target_x = 2.5

    interpolated_value = nevilles_interpolation(x, y, target_x)
    print(f"Neville插值结果: x = {target_x}, y = {interpolated_value}")

3.Newton插值

import numpy as np

def divided_differences(x, y):
    """
    计算并返回差商表。
    :param x: 自变量数组
    :param y: 因变量数组
    :return: 差商表
    """
    n = len(y)
    diff_table = np.zeros((n, n))
    diff_table[:, 0] = y  # 第一列为 y 值

    for j in range(1, n):
        for i in range(n - j):
            diff_table[i, j] = (diff_table[i + 1, j - 1] - diff_table[i, j - 1]) / (x[i + j] - x[i])

    return diff_table[0]  # 返回第一行的差商

def newton_interpolation(x, y, x_val):
    """
    使用 Newton 插值法计算插值值。
    :param x: 自变量数组
    :param y: 因变量数组
    :param x_val: 要插值的点
    :return: 插值结果
    """
    coeffs = divided_differences(x, y)  # 计算差商表的第一行

    n = len(coeffs)
    result = coeffs[0]
    product = 1.0

    for i in range(1, n):
        product *= (x_val - x[i - 1])  # 计算 (x - x0)(x - x1)...
        result += coeffs[i] * product  # 依次累加各项

    return result

# 示例用法
if __name__ == "__main__":
    # 已知数据点
    x_points = np.array([1.0, 2.0, 3.0, 4.0])
    y_points = np.array([1.0, 4.0, 9.0, 16.0])  # 对应 y = x^2

    # 计算插值点
    x_to_interpolate = 2.5
    interpolated_value = newton_interpolation(x_points, y_points, x_to_interpolate)

    print(f"插值点 x = {x_to_interpolate} 的插值结果为: {interpolated_value}")

4.Hermite插值

import numpy as np

def divided_differences_hermite(x, y, y_der):
    """
    计算 Hermite 插值的差商表。
    :param x: 自变量数组（包含重复点）
    :param y: 因变量数组（与 x 对应）
    :param y_der: 因变量的一阶导数数组（与 x 对应）
    :return: 差商表
    """
    n = len(x)
    diff_table = np.zeros((n, n))
    diff_table[:, 0] = y  # 第一列为 y 值

    for i in range(1, n):
        for j in range(n - i):
            if x[j] == x[j + i]:
                diff_table[j, i] = y_der[j] / np.math.factorial(i)  # 使用导数计算
            else:
                diff_table[j, i] = (diff_table[j + 1, i - 1] - diff_table[j, i - 1]) / (x[j + i] - x[j])

    return diff_table

def hermite_interpolation(x, y, y_der, x_val):
    """
    使用 Hermite 插值法计算插值值。
    :param x: 自变量数组（包含重复点）
    :param y: 因变量数组（与 x 对应）
    :param y_der: 因变量的一阶导数数组（与 x 对应）
    :param x_val: 要插值的点
    :return: 插值结果
    """
    # 构造 Hermite 差商表
    diff_table = divided_differences_hermite(x, y, y_der)
    
    # Hermite 插值多项式计算
    n = len(x)
    result = diff_table[0, 0]
    product = 1.0

    for i in range(1, n):
        product *= (x_val - x[i - 1])
        result += diff_table[0, i] * product

    return result

# 示例用法
if __name__ == "__main__":
    # 已知数据点及导数值
    x_points = np.array([1.0, 1.0, 2.0, 2.0])  # 包含重复点
    y_points = np.array([1.0, 1.0, 4.0, 4.0])  # 对应 y = x^2
    y_derivatives = np.array([0.0, 0.0, 4.0, 4.0])  # 对应 y' = 2x

    # 计算插值点
    x_to_interpolate = 1.5
    interpolated_value = hermite_interpolation(x_points, y_points, y_derivatives, x_to_interpolate)

    print(f"插值点 x = {x_to_interpolate} 的插值结果为: {interpolated_value}")