数值分析之插值法

import numpy as np


# Lagrange插值
def lagrange_interpolation(x, x_points, y_points):
    n = len(x_points)
    result = 0
    print(f"\n计算拉格朗日插值在 x = {x} 的值：")
    for i in range(n):
        # 计算第 i 个基函数 Li(x)
        Li = y_points[i]
        term = 1
        print(f"\n基函数 L{i}(x): y[{i}] = {y_points[i]}")
        for j in range(n):
            if j != i:
                term *= (x - x_points[j]) / (x_points[i] - x_points[j])
                print(
                    f"  项 (x - x[{j}]) / (x[{i}] - x[{j}]) = ({x} - {x_points[j]}) / ({x_points[i]} - {x_points[j]})")
        result += term * y_points[i]
        print(f"  L{i}(x) = {term:.4f}, 贡献 = {term * y_points[i]:.4f}")
    print(f"\n插值结果: P({x}) = {result:.4f}")
    return result


# 差商表计算
def divided_difference(x_points, y_points):
    n = len(x_points)
    table = np.zeros((n, n))
    table[:, 0] = y_points  # 第0阶差商
    print("\n差商表计算过程:")
    for j in range(1, n):
        for i in range(n - j):
            table[i, j] = (table[i + 1, j - 1] - table[i, j - 1]) / (x_points[i + j] - x_points[i])
            print(
                f"f[{i},{i + j}] = ({table[i + 1, j - 1]} - {table[i, j - 1]}) / ({x_points[i + j]} - {x_points[i]}) = {table[i, j]}")
    return table[0, :]  # 返回最高阶差商的第一行


def newton_interpolation(x, x_points, y_points):
    coeffs = divided_difference(x_points, y_points)
    n = len(x_points)
    result = coeffs[0]
    term = 1
    print(f"\n牛顿插值计算 P({x}):")
    print(f"  P0 = {coeffs[0]}")
    for i in range(1, n):
        term *= (x - x_points[i - 1])
        result += coeffs[i] * term
        print(f"  第{i}阶项: {coeffs[i]} * (x - x[{i - 1}])... = {coeffs[i]} * {term} = {coeffs[i] * term}")
    print(f"插值结果: P({x}) = {result:.4f}")
    return result


# Hermite插值计算
def hermite_interpolation(x, x_points, y_points, dy_points):
    n = len(x_points) * 2  # 2个点需要4次方程
    print(f"\nHermite插值计算 P({x}):")
    # 构造三次多项式 H(x) = a + bx + cx^2 + dx^3
    A = np.array([[1, x_points[0], x_points[0] ** 2, x_points[0] ** 3],
                  [0, 1, 2 * x_points[0], 3 * x_points[0] ** 2],
                  [1, x_points[1], x_points[1] ** 2, x_points[1] ** 3],
                  [0, 1, 2 * x_points[1], 3 * x_points[1] ** 2]])
    b = np.array([y_points[0], dy_points[0], y_points[1], dy_points[1]])
    coeffs = np.linalg.solve(A, b)
    print(f"系数解: a = {coeffs[0]}, b = {coeffs[1]}, c = {coeffs[2]}, d = {coeffs[3]}")

    # 计算 H(x)
    result = coeffs[0] + coeffs[1] * x + coeffs[2] * x ** 2 + coeffs[3] * x ** 3
    print(f"  H({x}) = {coeffs[0]} + {coeffs[1]}*{x} + {coeffs[2]}*{x}^2 + {coeffs[3]}*{x}^3")
    print(f"插值结果: H({x}) = {result:.4f}")
    return result
# 分段线性插值计算
def piecewise_linear_interpolation(x, x_points, y_points):
    print(f"\n分段线性插值计算 P({x}):")
    for i in range(len(x_points) - 1):
        if x_points[i] <= x <= x_points[i + 1]:
            slope = (y_points[i + 1] - y_points[i]) / (x_points[i + 1] - x_points[i])
            result = y_points[i] + slope * (x - x_points[i])
            print(f"  区间 [{x_points[i]}, {x_points[i+1]}]:")
            print(f"    斜率 = ({y_points[i+1]} - {y_points[i]}) / ({x_points[i+1]} - {x_points[i]}) = {slope}")
            print(f"    P({x}) = {y_points[i]} + {slope} * ({x} - {x_points[i]}) = {result:.4f}")
            return result
    return None  # x 超出范围

# 三次样条插值计算
def cubic_spline_interpolation(x, x_points, y_points):
    n = len(x_points) - 1
    h = np.diff(x_points)  # 步长
    a = y_points[:-1]  # 系数 a
    print(f"\n三次样条插值计算 P({x}):")

    # 构造三对角矩阵求解二阶导数 m
    A = np.zeros((n + 1, n + 1))
    b = np.zeros(n + 1)
    A[0, 0] = 1  # 自然边界条件 m0 = 0
    A[n, n] = 1  # mn = 0
    for i in range(1, n):
        A[i, i - 1] = h[i - 1]
        A[i, i] = 2 * (h[i - 1] + h[i])
        A[i, i + 1] = h[i]
        b[i] = 6 * ((y_points[i + 1] - y_points[i]) / h[i] - (y_points[i] - y_points[i - 1]) / h[i - 1])
    m = np.linalg.solve(A, b)
    print(f"二阶导数 m = {m}")

    # 计算插值
    for i in range(n):
        if x_points[i] <= x <= x_points[i + 1]:
            t = (x - x_points[i]) / h[i]
            result = (1 - t) * y_points[i] + t * y_points[i + 1] + \
                     h[i] ** 2 / 6 * ((1 - t) ** 3 - (1 - t)) * m[i] + \
                     h[i] ** 2 / 6 * (t ** 3 - t) * m[i + 1]
            print(f"  区间 [{x_points[i]}, {x_points[i + 1]}]: t = {t}")
            print(f"插值结果: S({x}) = {result:.4f}")
            return result
    return None
if __name__ == '__main__':
    # Lagrange测试
    # 数据点
    x_points1 = np.array([0.32, 0.34, 0.36])
    y_points1 = np.array([0.314567, 0.333487, 0.352274])
    x_test1 = 0.3367
    lagrange_interpolation(x_test1, x_points1, y_points1)
    # Newton测试
    # 数据点
    x_points2 = np.array([0, 1, 2, 3], dtype=float)
    y_points2 = np.array([1, 1.6487, 2.7183, 4.4817], dtype=float)
    x_test2 = 2.8
    newton_interpolation(x_test2, x_points2, y_points2)
    # Hermite插值计算测试
    # 数据点和导数
    x_points3 = np.array([1, 2])
    y_points3 = np.array([2, 3])
    dy_points3 = np.array([1, -1])  # 一阶导数值
    x_test3 = 0.5
    hermite_interpolation(x_test3, x_points3, y_points3, dy_points3)
    # 分段线性插值计算测试
    # 数据点
    x_points4 = np.array([-3,-1,2, 3, 9],dtype=float)
    y_points4 = np.array([12, 5, 1,6,12])
    x_test4 = 1.2
    piecewise_linear_interpolation(x_test4, x_points4, y_points4)
    x_test5 = 3.3
    piecewise_linear_interpolation(x_test5, x_points4, y_points4)
    # 三次样条插值计算
    # 数据点
    x_points5 = np.array([1, 2, 4,5])
    y_points5 = np.array([1, 3, 4,2])
    x_test6 = 3
    cubic_spline_interpolation(x_test6, x_points5, y_points5)