实验名称: 插值算法的实现
实验目的:
a. 实现拉格朗日插值;
b. 验证随着插值结点的增多插值曲线的变化情况。
c. 实现牛顿插值,显示插商结果;
d. 比较拉格朗日插值与牛顿插值的插值结果是否相同。
实验内容
一、实现拉格朗日插值
为了方便起见,我选择的函数为三角函数
f
(
x
)
=
s
i
n
x
f(x) = sinx
f(x)=sinx
运用基函数实现拉格朗日插值,首先要计算基函数
基函数的一般形式:
根据基函数分子和分母的形式我们可以看出,如果想要求解基函数,可以先构造一个累乘函数,再求基函数
# 累乘函数
def T(x, i, X):
T1 = 1
for j in X:
if X[i] == j:
continue
T1 = T1 * (x-j)
return T1
# 插值基函数
def P(i, x, X, Y):
P1 = T(x, i, X)/T(X[i], i, X) * Y[i]
return P1
最后根据公式计算拉格朗日插值
选取八个点得到的拉格朗日插值
二、验证随着插值结点的增多插值曲线的变化情况
当插值节点为2个时
当插值节点为4个时
当插值节点为6个时
当插值节点为8个时
根据以上结果可以发现,随着插值节点的增多,曲线的拟合效果越好。为了更加精确的表示,将误差做成下列表格,可以看出随着插值节点的增多,误差减小。
插值节点 | 误差 |
---|---|
2个 | 0.11861563625463706 |
4个 | 0.014743131588928482 |
6个 | 0.0001773516321779134 |
8个 | 4.743257968842229e-06 |
三、实现牛顿插值,显示插商结果
根据差商的定义,容易得出高阶差商是两个低一阶差商的差商
根据上图,容易联想到或许差商可以用二位数组来实现,具体实现代码如下:
#计算差商
def chashang(X ,Y):
A = np.zeros([n,n])
for i in range(0,n):
A[i][0] = Y[i] #每一行的第一个都是f(xi)
for j in range(1,n):
for i in range(j,n):
A[i][j] = (A[i][j-1] - A[i-1][j-1])/(X[i]-X[i-j]) #根据差商公式计算
return A
根据公式实现牛顿插值:
#计算牛顿插值
def N(x, X, Y):
sum = Y[0]
t = np.zeros((len(X),len(X)))
for i in range(0,len(X)):
t[i,0] = Y[i]
sumt = 1.0
for i in range(1,len(X)):
sumt = sumt*(x-X[i-1])
for j in range(i,len(X)):
sumt[j,i] = (sumt[j,i-1] - sumt[j-1,i-1])/(X[j] - X[j-i])
sum += sumt*t[i,i]
return sum
四、比较拉格朗日插值与牛顿插值的插值结果是否相同
根据上式及以上所有结果可知,一般拉格朗日插值法和牛顿插值法的结果一致。
值得注意的是,当插值多项式从n-1次增加到n次时,拉格朗日插值必须重新计算所有的基本插值多项式;而对于牛顿插值,只需要用表格再计算一个n阶均差,然后再加上一项即可。
实验体会
1.通过用python实现拉格朗日插值方法和牛顿插值方法,加深了对插值方法的理解;
2.一般情况下在插值节点数目较多时,得到的插值公式的精度越高,插值结果误差越小,但也有特殊情况;
3.用插值方法推导已知区间之外的点得到的值往往不准确;
4.通过该实验,我学习并熟练掌握了python中的numpy和matlab库,更加熟悉了python语言的使用。