133 lines
4.1 KiB
Python
133 lines
4.1 KiB
Python
import math
|
||
import numpy as np
|
||
|
||
# 定义插值点
|
||
list_x = [0.0, 0.2, 0.4, 0.6, 0.8] # 自变量列表
|
||
list_y = [1.0, 1.2214, 1.4918, 1.8221, 2.2255] # 因变量列表
|
||
|
||
def GetDyList(list_y):
|
||
"""
|
||
计算前向差分表。
|
||
|
||
参数:
|
||
list_y (list of float): 因变量列表。
|
||
|
||
返回:
|
||
list of list of float: 前向差分表，每一行是对应阶数的差分值。
|
||
"""
|
||
result = []
|
||
result.append(list_y)
|
||
for i in range(len(list_y) - 1, 0, -1):
|
||
tmp = []
|
||
for j in range(i):
|
||
tmp.append(result[len(list_y) - 1 - i][j + 1] - result[len(list_y) - 1 - i][j])
|
||
result.append(tmp)
|
||
return result
|
||
|
||
def FxDiff_n(x, n):
|
||
"""
|
||
计算函数 e^x 的 n 阶导数值。
|
||
|
||
参数:
|
||
x (float): 自变量值。
|
||
n (int): 导数的阶数，n >= 0。
|
||
|
||
返回:
|
||
float: e^x 的 n 阶导数值。
|
||
"""
|
||
return math.exp(x)
|
||
|
||
def NewtonForwardRegression(x, n, list_x):
|
||
"""
|
||
计算牛顿前向插值的余项估计值。
|
||
|
||
参数:
|
||
x (float): 待插值的自变量值。
|
||
n (int): 插值多项式的阶数。
|
||
list_x (list of float): 自变量列表，要求等间距。
|
||
|
||
返回:
|
||
float: 牛顿前向插值的余项估计值。
|
||
"""
|
||
h = list_x[1] - list_x[0] # 等间距
|
||
t = (x - list_x[0]) / h # 归一化变量
|
||
ks = max([abs(FxDiff_n(list_x[i], n + 1)) for i in range(n + 1)]) # n+1 阶导数的最大值
|
||
result = 1
|
||
for i in range(n + 1):
|
||
result *= (t - i) * h / (i + 1) # 余项公式中的乘积项
|
||
return result * ks
|
||
|
||
def NewtonForwardInterpolation(x, n, list_x, list_y):
|
||
"""
|
||
使用牛顿前向插值法计算插值结果及其余项估计值。
|
||
|
||
参数:
|
||
x (float): 待插值的自变量值。
|
||
n (int): 插值多项式的阶数。
|
||
list_x (list of float): 自变量列表，要求等间距。
|
||
list_y (list of float): 因变量列表，与 list_x 一一对应。
|
||
|
||
返回:
|
||
tuple: 插值结果和余项估计值 (result, r)。
|
||
"""
|
||
list_dyk = GetDyList(list_y) # 差分表
|
||
h = list_x[1] - list_x[0] # 等间距
|
||
t = (x - list_x[0]) / h # 归一化变量
|
||
result = list_y[0] # 插值结果初值
|
||
mul = 1 # 乘积项
|
||
for i in range(1, n):
|
||
mul *= ((t - i + 1) / i) # 计算乘积项
|
||
result += list_dyk[i][0] * mul # 累加差分项
|
||
r = abs(NewtonForwardRegression(x, n, list_x)) # 余项估计
|
||
return (result, r)
|
||
|
||
def GetDQyList(list_x, list_y):
|
||
"""
|
||
计算牛顿基函数插值的差商表。
|
||
|
||
参数:
|
||
list_x (list of float): 自变量列表。
|
||
list_y (list of float): 因变量列表。
|
||
|
||
返回:
|
||
list of list of float: 差商表，每一行是对应阶数的差商值。
|
||
"""
|
||
result = []
|
||
result.append(list_y)
|
||
for i in range(len(list_y) - 1, 0, -1):
|
||
tmp = []
|
||
for j in range(i):
|
||
tmp.append((result[len(list_y) - 1 - i][j + 1] - result[len(list_y) - 1 - i][j]) /
|
||
(list_x[j + len(list_y) - i] - list_x[j]))
|
||
result.append(tmp)
|
||
return result
|
||
|
||
def NewtonBaseInterpolation(x, n, list_x, list_y):
|
||
"""
|
||
使用牛顿基函数插值法计算插值结果。
|
||
|
||
参数:
|
||
x (float): 待插值的自变量值。
|
||
n (int): 插值多项式的阶数。
|
||
list_x (list of float): 自变量列表。
|
||
list_y (list of float): 因变量列表，与 list_x 一一对应。
|
||
|
||
返回:
|
||
float: 插值结果。
|
||
"""
|
||
list_dqy = GetDQyList(list_x, list_y) # 差商表
|
||
result = list_dqy[0][0] # 插值结果初值
|
||
mul = 1 # 乘积项
|
||
for i in range(1, n):
|
||
mul *= (x - list_x[i - 1]) # 计算乘积项
|
||
result += list_dqy[i][0] * mul # 累加差商项
|
||
return result
|
||
|
||
if __name__ == '__main__':
|
||
# 测试牛顿前向插值
|
||
print(NewtonForwardInterpolation(0.12, 2, list_x, list_y))
|
||
print(NewtonForwardInterpolation(0.12, 3, list_x, list_y))
|
||
# 测试牛顿基函数插值
|
||
print(NewtonBaseInterpolation(0.12, 2, list_x, list_y))
|
||
print(NewtonBaseInterpolation(0.12, 3, list_x, list_y))
|