第三十六篇 正向差分的插值求解
差分法
另一种求插值多项式的方法是通过np个数据点(xi, yi), i = 0,1,2,…,n(其中n = np−1),首先将插值多项式写成另外一种形式:
其中常数Ci, i = 0,1,2,…,n可以按照要求来确定
重新计算之后得到
计算实例
本文重复使用上篇中给出的例子。使用差分法去推导通过这些点的多项式
然后考虑额外的点去修正多项式
计算x=4.5时的y值
这个方程的第一部分牵涉到3个数据点,因此插值多项式将是二次型的。首先计算方程的常量
然后带入方程得到
这和上篇中的结果相同。
第四个数据点x3 = 5,y3 = 9代入将得到一个三次插值多项式。额外的常量将得到
然后将这个三次项加入之前的二次项中
或
三次项展示在下图然后得到x=4.5和y=8.175。
等间隔的差分法
如果数据以等距的x值提供,使xi−xi−1 = h,则对系数Ci的求导得到了很大程度的简化。
进行求解之前,介绍一个‘从前差分’的概念
和一个‘从后差分’的概念
而且,有一个‘二次差分’,也就是‘差分中的差分’,第二次从前差分能写成
一般形式为
在这次程序中,主要使用’从前差分’的方法
对于之前的方程,系数为
它能展现成一个通用形式
从之前可以得出Δjy0项可以通过递归来求值。为了便于手算,上表已经给出了计算步骤,其中有六个数据点
上表描述的特殊格式称为(牛顿)正向差分,其特征是沿从左上角到右下角向下倾斜的对角线下标保持相同。其他的布局是可能的,如逆向差分和高斯方法,但这些将不在这里讨论。当把表中的数据代入方程,“零”差项,也就是刚开始的差值项,是y值本身,因此
计算实例
给与下面的数据是函数y = cosx 的数据,求cons27
首先计算正向差分的数据,写入表格
一般方程的插值多项式,x0可以选择任意的初始值x。然而,如果我们希望包括所有五个数据点,生成一个四阶插值多项式, x0应设置等于表中的最高值,也就是说,x0 = 20。x值之间的常量间隔为h = 5,系数可以用以表格形式计算如下
保留小数点后五位,方程的插值多项式可以写成
扩展得到
替换需要求解的值,保留到小数点后五位得到
从上面看出,从上面的Q4(x)看出C3和C4项几乎对多项式没有作用。
如果系数变得足够小,截断插值多项式的能力是差分方法的一个有用特征。它不仅节省了计算的工作量,而且还提供了对数据点理论上起源的物理上的观测视角。
在拉格朗日方法中,这样的保存是不可能的,无论是否需要,完整的n阶多项式必须导出。拉格朗日方法确实有简化的优点,但是,是在只提供几个点,相对容易手算的情况下。
程序如下
#正向差分的插值求解
import numpy as np
npt=4;xi=2.8
c=np.zeros((npt))
diff=np.zeros((npt,npt))
x=np.zeros((npt))
x[0:npt]=(0.0,1.0,2.0,3.0)
diff[0:npt,0]=(1.0,1.0,15.0,61.0)
print('正向差分的插值求解')
print('数据点',' x y')
for i in range(1,npt+1):
print('{:13.4e}'.format(x[i-1]),end='')
print('{:13.4e}'.format(diff[i-1,0]))
h=x[1]-x[0]
for i in range(1,npt):
for j in range(1,npt-i+1):
diff[j-1,i]=diff[j,i-1]-diff[j-1,i-1]
c[0]=diff[0,0];yi=c[0];term=1.0;factorial=1.0
for i in range(1,npt):
factorial=factorial*i;c[i]=diff[0,i]/(factorial*h**i)
term=term*(xi-x[i-1]);yi=yi+term*c[i]
print('多项式系数C')
for i in range(1,npt+1):
print('{:13.4e}'.format(c[i-1]))
print('插值点',' x y')
print('{:13.4e}'.format(xi),end='')
print('{:13.4e}'.format(yi))
终端输出结果如下