数值分析之数值积分

import numpy as np


# 线性插值
def linear_interpolation(x_data, y_data, x):
    n = len(x_data)  # 获取数据点的数量
    for i in range(n - 1):  # 遍历每个区间
        if x_data[i] <= x <= x_data[i + 1]:  # 检查 x 是否在当前区间内
            # 计算斜率
            slope = (y_data[i + 1] - y_data[i]) / (x_data[i + 1] - x_data[i])
            # 计算插值结果
            y = y_data[i] + slope * (x - x_data[i])

            # 打印详细计算步骤
            print(f"线性插值计算细节:")
            print(f"使用点: ({x_data[i]}, {y_data[i]}), ({x_data[i + 1]}, {y_data[i + 1]})")
            print(
                f"斜率: (y_data[i+1] - y_data[i]) / (x_data[i+1] - x_data[i]) = ({y_data[i + 1]} - {y_data[i]}) / ({x_data[i + 1]} - {x_data[i]}) = {slope}")
            print(f"结果: y = y_data[i] + 斜率 * (x - x_data[i]) = {y_data[i]} + {slope} * ({x} - {x_data[i]}) = {y}")

            return y  # 返回插值结果
    return None  # 如果 x 不在任何区间内，返回 None


# 二次插值=抛物线插值
def quadratic_interpolation(x_data, y_data, x):
    n = len(x_data)  # 获取数据点的数量
    for i in range(n - 2):  # 遍历每个区间
        if x_data[i] <= x <= x_data[i + 2]:  # 检查 x 是否在当前区间内
            # 使用三个点 (x_data[i], y_data[i]), (x_data[i+1], y_data[i+1]), (x_data[i+2], y_data[i+2]) 进行二次插值
            A = np.array([
                [1, x_data[i], x_data[i] ** 2],
                [1, x_data[i + 1], x_data[i + 1] ** 2],
                [1, x_data[i + 2], x_data[i + 2] ** 2]
            ])  # 构建系数矩阵 A
            B = np.array([y_data[i], y_data[i + 1], y_data[i + 2]])  # 构建常数矩阵 B
            coeffs = np.linalg.solve(A, B)  # 解线性方程组求解系数 a, b, c
            a, b, c = coeffs  # 提取系数 a, b, c

            # 打印详细计算步骤
            print(f"二次插值计算细节:")
            print(
                f"使用点: ({x_data[i]}, {y_data[i]}), ({x_data[i + 1]}, {y_data[i + 1]}), ({x_data[i + 2]}, {y_data[i + 2]})")
            print(f"系数矩阵 A:\n{A}")
            print(f"常数矩阵 B:\n{B}")
            print(f"解得系数 a, b, c: {a}, {b}, {c}")

            y = a + b * x + c * x ** 2  # 计算插值结果
            print(f"结果: y = {a} + {b} * {x} + {c} * {x}^2 = {y}")
            return y  # 返回插值结果
    return None  # 如果 x 不在任何区间内，返回 None


# Lagrange插值
def lagrange_interpolation(x_data, y_data, x):
    n = len(x_data)  # 获取数据点的数量
    y = 0  # 初始化插值结果
    for i in range(n):  # 遍历每个数据点
        l = 1  # 初始化基函数 L_i(x)
        for j in range(n):  # 计算基函数 L_i(x)
            if i != j:
                l *= (x - x_data[j]) / (x_data[i] - x_data[j])  # 计算 L_i(x)
        y += y_data[i] * l  # 累加每个基函数的结果

        # 打印详细计算步骤
        print(f"Lagrange插值计算细节:")
        print(f"使用点: ({x_data[i]}, {y_data[i]})")
        print(f"L_i(x) = {l}")
        print(f"y += y_data[i] * L_i(x) = {y_data[i]} * {l} = {y_data[i] * l}")

    print(f"最终结果: y = {y}")
    return y  # 返回插值结果


# 牛顿法
def divided_difference(x_data, y_data):
    n = len(x_data)  # 获取数据点的数量
    coef = y_data.copy()  # 复制 y_data 作为初始系数
    for j in range(1, n):  # 计算差分系数
        for i in range(n - 1, j - 1, -1):
            coef[i] = (coef[i] - coef[i - 1]) / (x_data[i] - x_data[i - j])  # 更新系数
    return coef  # 返回差分系数


def newton_interpolation(x_data, y_data, x):
    coef = divided_difference(x_data, y_data)  # 获取差分系数
    n = len(x_data)  # 获取数据点的数量
    y = coef[-1]  # 初始化插值结果为最后一个系数
    print(f"牛顿插值计算细节:")
    print(f"初始系数: {coef}")
    for i in range(n - 2, -1, -1):  # 计算插值结果
        y = y * (x - x_data[i]) + coef[i]  # 更新插值结果
        print(f"y = y * (x - x_data[i]) + coef[i] = {y} * ({x} - {x_data[i]}) + {coef[i]} = {y}")
    print(f"最终结果: y = {y}")
    return y  # 返回插值结果


# Hermite插值
def L(x, i):
    l = 1  # 初始化基函数 L_i(x)
    for j in range(len(x_data)):  # 计算基函数 L_i(x)
        if i != j:
            l *= (x - x_data[j]) / (x_data[i] - x_data[j])  # 更新 L_i(x)
    return l  # 返回基函数 L_i(x)


def hermite_interpolation(x_data, y_data, y_prime_data, x):
    n = len(x_data)  # 获取数据点的数量
    H = 0  # 初始化插值结果
    for i in range(n):  # 遍历每个数据点
        L_i = L(x, i)  # 计算基函数 L_i(x)
        h1 = 1 - 2 * L_i * (x - x_data[i])  # 计算 h1
        h2 = (x - x_data[i]) * L_i  # 计算 h2
        H += y_data[i] * h1 * L_i * L_i + y_prime_data[i] * h2 * L_i * L_i  # 累加每个基函数的结果

        # 打印详细计算步骤
        print(f"Hermite插值计算细节:")
        print(f"使用点: ({x_data[i]}, {y_data[i]}), 导数: {y_prime_data[i]}")
        print(f"L_i(x) = {L_i}")
        print(f"h1 = 1 - 2 * L_i * (x - x_data[i]) = 1 - 2 * {L_i} * ({x} - {x_data[i]}) = {h1}")
        print(f"h2 = (x - x_data[i]) * L_i = ({x} - {x_data[i]}) * {L_i} = {h2}")
        print(
            f"H += y_data[i] * h1 * L_i * L_i + y_prime_data[i] * h2 * L_i * L_i = {y_data[i]} * {h1} * {L_i} * {L_i} + {y_prime_data[i]} * {h2} * {L_i} * {L_i} = {H}")

    print(f"最终结果: H = {H}")
    return H  # 返回插值结果


# 分段线性插值
def piecewise_linear_interpolation(x_data, y_data, x):
    n = len(x_data)  # 获取数据点的数量
    for i in range(n - 1):  # 遍历每个区间
        if x_data[i] <= x <= x_data[i + 1]:  # 检查 x 是否在当前区间内
            # 计算斜率
            slope = (y_data[i + 1] - y_data[i]) / (x_data[i + 1] - x_data[i])
            # 计算插值结果
            y = y_data[i] + slope * (x - x_data[i])

            # 打印详细计算步骤
            print(f"分段线性插值计算细节:")
            print(f"使用点: ({x_data[i]}, {y_data[i]}), ({x_data[i + 1]}, {y_data[i + 1]})")
            print(
                f"斜率: (y_data[i+1] - y_data[i]) / (x_data[i+1] - x_data[i]) = ({y_data[i + 1]} - {y_data[i]}) / ({x_data[i + 1]} - {x_data[i]}) = {slope}")
            print(f"结果: y = y_data[i] + 斜率 * (x - x_data[i]) = {y_data[i]} + {slope} * ({x} - {x_data[i]}) = {y}")

            return y  # 返回插值结果
    return None  # 如果 x 不在任何区间内，返回 None


# 三次样条插值
def cubic_spline_interpolation(x_data, y_data, x):
    n = len(x_data)  # 获取数据点的数量
    h = np.diff(x_data)  # 计算相邻数据点之间的差值
    alpha = np.zeros(n - 1)  # 初始化 alpha 数组
    for i in range(1, n - 1):  # 计算 alpha 数组
        alpha[i] = 3.0 / h[i] * (y_data[i + 1] - y_data[i]) - 3.0 / h[i - 1] * (y_data[i] - y_data[i - 1])

    l = np.zeros(n)  # 初始化 l 数组
    mu = np.zeros(n)  # 初始化 mu 数组
    z = np.zeros(n)  # 初始化 z 数组
    l[0] = 2.0 * (h[0] + h[1])  # 计算 l[0]
    mu[0] = h[1] / l[0]  # 计算 mu[0]
    z[0] = alpha[0] / l[0]  # 计算 z[0]
    for i in range(1, n - 1):  # 计算 l, mu, z 数组
        l[i] = 2.0 * (h[i - 1] + h[i]) - h[i - 1] * mu[i - 1]  # 计算 l[i]
        mu[i] = h[i] / l[i]  # 计算 mu[i]
        z[i] = (alpha[i] - h[i - 1] * z[i - 1]) / l[i]  # 计算 z[i]
    l[n - 1] = h[n - 2] * (2.0 - mu[n - 2])  # 计算 l[n-1]
    z[n - 1] = (alpha[n - 2] - h[n - 2] * z[n - 2]) / l[n - 1]  # 计算 z[n-1]
    c = np.zeros(n)  # 初始化 c 数组
    b = np.zeros(n - 1)  # 初始化 b 数组
    d = np.zeros(n - 1)  # 初始化 d 数组
    c[n - 1] = z[n - 1]  # 设置 c[n-1]
    for i in range(n - 2, -1, -1):  # 计算 c, b, d 数组
        c[i] = z[i] - mu[i] * c[i + 1]  # 计算 c[i]
        b[i] = (y_data[i + 1] - y_data[i]) / h[i] - h[i] * (c[i + 1] + 2.0 * c[i]) / 3.0  # 计算 b[i]
        d[i] = (c[i + 1] - c[i]) / (3.0 * h[i])  # 计算 d[i]

    for i in range(n - 1):  # 遍历每个区间
        if x_data[i] <= x <= x_data[i + 1]:  # 检查 x 是否在当前区间内
            dx = x - x_data[i]  # 计算 dx
            # 计算插值结果
            y = y_data[i] + b[i] * dx + c[i] * dx * dx + d[i] * dx * dx * dx

            # 打印详细计算步骤
            print(f"三次样条插值计算细节:")
            print(f"使用区间: [{x_data[i]}, {x_data[i + 1]}]")
            print(f"系数: b[i] = {b[i]}, c[i] = {c[i]}, d[i] = {d[i]}")
            print(
                f"结果: y = y_data[i] + b[i] * dx + c[i] * dx * dx + d[i] * dx * dx * dx = {y_data[i]} + {b[i]} * {dx} + {c[i]} * {dx} * {dx} + {d[i]} * {dx} * {dx} * {dx} = {y}")

            return y  # 返回插值结果
    return None  # 如果 x 不在任何区间内，返回 None


# 示例数据
x_data = np.array([1, 2, 3])  # 定义 x 数据点
y_data = np.array([2, 3, 5])  # 定义 y 数据点
x = 1.5  # 定义需要插值的 x 值
print(f"线性插值法在 x = {x} 处的结果: {linear_interpolation(x_data, y_data, x)}\n")  # 打印线性插值结果

print(f"抛物线插值=二次插值法在 x = {x} 处的结果: {quadratic_interpolation(x_data, y_data, x)}\n")  # 打印二次插值结果

print(f"拉格朗日插值法在 x = {x} 处的结果: {lagrange_interpolation(x_data, y_data, x)}\n")  # 打印拉格朗日插值结果

print(f"牛顿插值法在 x = {x} 处的结果: {newton_interpolation(x_data, y_data, x)}\n")  # 打印牛顿插值结果

print(f"分段线性插值法在 x = {x} 处的结果: {piecewise_linear_interpolation(x_data, y_data, x)}\n")  # 打印分段线性插值结果

print(f"三次样条插值法在 x = {x} 处的结果: {cubic_spline_interpolation(x_data, y_data, x)}\n")  # 打印三次样条插值结果

# 埃尔米特插值法示例数据
x_data1 = np.array([1, 2])  # 定义 x 数据点
y_data1 = np.array([2, 3])  # 定义 y 数据点
y_prime_data1 = np.array([1, 1])  # 定义导数数据点
x = 1.5  # 定义需要插值的 x 值
print(f"埃尔米特插值法在 x = {x} 处的结果: {hermite_interpolation(x_data1, y_data1, y_prime_data1, x)}\n")  # 打印埃尔米特插值结果