基于Python实现Hermite插值求解函数近似值的习题指导

问题说明 #

已知函数 $f(x)=3xe^{x}-e{2x}$，完成以下两项任务：

• a) 基于节点 $x_0=1$、$x_1=1.05$ 构造三次Hermite插值多项式，近似计算 $f(1.03)$，并对比实际误差与理论误差界；
• b) 基于节点 $x_0=1$、$x_1=1.05$、$x_2=1.07$ 构造五次Hermite插值多项式，重复(a)的操作。

实现步骤与代码 #

1. 定义函数及各阶导数 #

首先定义原函数、各阶导数（用于插值计算和误差界推导）：

import numpy as np
from scipy.optimize import minimize_scalar

# 原函数
def f(x):
    return 3 * x * np.exp(x) - np.exp(2 * x)

# 一阶导数
def f1(x):
    return 3 * np.exp(x) + 3 * x * np.exp(x) - 2 * np.exp(2 * x)

# 二阶导数
def f2(x):
    return 6 * np.exp(x) + 3 * x * np.exp(x) - 4 * np.exp(2 * x)

# 四阶导数（用于三次Hermite误差界）
def f4(x):
    return 12 * np.exp(x) + 3 * x * np.exp(x) - 16 * np.exp(2 * x)

# 六阶导数（用于五次Hermite误差界）
def f6(x):
    return 18 * np.exp(x) + 3 * x * np.exp(x) - 64 * np.exp(2 * x)

2. 三次Hermite插值（两点） #

插值计算实现

# 三次Hermite插值函数
def cubic_hermite(x, x0, x1, f0, f1, f0_prime, f1_prime):
    h = x1 - x0
    # 计算Hermite基函数
    alpha0 = (1 - 2*(x - x0)/(-h)) * ((x - x1)/(-h))**2
    alpha1 = (1 - 2*(x - x1)/h) * ((x - x0)/h)**2
    beta0 = (x - x0) * ((x - x1)/(-h))**2
    beta1 = (x - x1) * ((x - x0)/h)**2
    # 插值结果
    return f0*alpha0 + f1*alpha1 + f0_prime*beta0 + f1_prime*beta1

# 节点与目标计算点
x0, x1 = 1.0, 1.05
x_target = 1.03

# 计算节点的函数值和导数值
f0, f1 = f(x0), f(x1)
f0_prime, f1_prime = f1(x0), f1(x1)

# 执行插值计算
h3_val = cubic_hermite(x_target, x0, x1, f0, f1, f0_prime, f1_prime)
true_val = f(x_target)
actual_error_cubic = abs(true_val - h3_val)

误差界计算

需要找到四阶导数在区间$[1,1.05]$上的绝对值最大值，代入误差界公式：

# 定义四阶导数绝对值函数，用于求区间最大值
def abs_f4(x):
    return abs(f4(x))

# 求解区间内的最大值
res = minimize_scalar(lambda x: -abs_f4(x), bounds=(x0, x1), method='bounded')
M4 = -res.fun

# 计算理论误差界
error_bound_cubic = M4 / 24 * abs((x_target - x0)**2 * (x_target - x1)**2)

结果输出

print("=== 三次Hermite插值结果 ===")
print(f"f(1.03)实际值: {true_val:.8f}")
print(f"插值近似值: {h3_val:.8f}")
print(f"实际误差: {actual_error_cubic:.8e}")
print(f"理论误差界: {error_bound_cubic:.8e}")

3. 五次Hermite插值（三点） #

五次Hermite插值需要每个节点的函数值、一阶导、二阶导（共6个条件），通过构造方程组求解五次多项式系数：

# 构造五次Hermite插值的系数方程组
def quintic_hermite_coeffs(x_nodes, f_vals, f1_vals, f2_vals):
    # 五次多项式有6个系数a0~a5
    A = np.zeros((6, 6))
    b = np.zeros(6)
    idx = 0
    for i in range(len(x_nodes)):
        x = x_nodes[i]
        # 函数值条件
        A[idx] = [x**k for k in range(6)]
        b[idx] = f_vals[i]
        idx +=1
        # 一阶导数条件
        A[idx] = [k * x**(k-1) if k>0 else 0 for k in range(6)]
        b[idx] = f1_vals[i]
        idx +=1
        # 二阶导数条件
        A[idx] = [k*(k-1)*x**(k-2) if k>=2 else 0 for k in range(6)]
        b[idx] = f2_vals[i]
        idx +=1
    # 求解系数
    coeffs = np.linalg.solve(A, b)
    return coeffs

# 节点与目标计算点
x_nodes = [1.0, 1.05, 1.07]
x_target = 1.03

# 计算节点的函数值、一阶导、二阶导
f_vals = [f(x) for x in x_nodes]
f1_vals = [f1(x) for x in x_nodes]
f2_vals = [f2(x) for x in x_nodes]

# 求解五次多项式系数
coeffs = quintic_hermite_coeffs(x_nodes, f_vals, f1_vals, f2_vals)

# 计算插值值
h5_val = sum(coeffs[k] * x_target**k for k in range(6))
actual_error_quintic = abs(true_val - h5_val)

误差界计算

找到六阶导数在区间$[1,1.07]$上的绝对值最大值，代入误差界公式：

# 定义六阶导数绝对值函数
def abs_f6(x):
    return abs(f6(x))

# 求解区间内的最大值
res_quintic = minimize_scalar(lambda x: -abs_f6(x), bounds=(min(x_nodes), max(x_nodes)), method='bounded')
M6 = -res_quintic.fun

# 计算理论误差界
error_bound_quintic = M6 / 720 * abs((x_target - x_nodes[0])**2 * (x_target - x_nodes[1])**2 * (x_target - x_nodes[2])**2)

结果输出

print("\n=== 五次Hermite插值结果 ===")
print(f"f(1.03)实际值: {true_val:.8f}")
print(f"插值近似值: {h5_val:.8f}")
print(f"实际误差: {actual_error_quintic:.8e}")
print(f"理论误差界: {error_bound_quintic:.8e}")

结果说明 #

• 三次Hermite插值的实际误差远小于理论误差界，因为误差界是基于区间导数最大值的保守估计；
• 五次Hermite插值引入了更多节点和高阶导数信息，误差会比三次插值更小，符合数值插值的收敛特性。

内容的提问来源于stack exchange，提问作者angelabayona