第三十六篇 正向差分的插值求解

差分法

另一种求插值多项式的方法是通过np个数据点(xi, yi), i = 0,1,2,…,n(其中n = np−1),首先将插值多项式写成另外一种形式:

python diff 差分 python求解差分方程_python diff 差分


其中常数Ci, i = 0,1,2,…,n可以按照要求来确定

python diff 差分 python求解差分方程_python_02


重新计算之后得到

python diff 差分 python求解差分方程_python diff 差分_03


计算实例

本文重复使用上篇中给出的例子。使用差分法去推导通过这些点的多项式

python diff 差分 python求解差分方程_差分_04


然后考虑额外的点去修正多项式

python diff 差分 python求解差分方程_python diff 差分_05


计算x=4.5时的y值

这个方程的第一部分牵涉到3个数据点,因此插值多项式将是二次型的。首先计算方程的常量

python diff 差分 python求解差分方程_python_06


然后带入方程得到

python diff 差分 python求解差分方程_差分_07


这和上篇中的结果相同。

第四个数据点x3 = 5,y3 = 9代入将得到一个三次插值多项式。额外的常量将得到

python diff 差分 python求解差分方程_差分_08


然后将这个三次项加入之前的二次项中

python diff 差分 python求解差分方程_差分_09



python diff 差分 python求解差分方程_差分_10


三次项展示在下图然后得到x=4.5和y=8.175。

等间隔的差分法

如果数据以等距的x值提供,使xi−xi−1 = h,则对系数Ci的求导得到了很大程度的简化。

进行求解之前,介绍一个‘从前差分’的概念

python diff 差分 python求解差分方程_python_11


和一个‘从后差分’的概念

python diff 差分 python求解差分方程_差分_12


而且,有一个‘二次差分’,也就是‘差分中的差分’,第二次从前差分能写成

python diff 差分 python求解差分方程_python diff 差分_13


一般形式为

python diff 差分 python求解差分方程_线性代数_14


在这次程序中,主要使用’从前差分’的方法

对于之前的方程,系数为

python diff 差分 python求解差分方程_python_15


它能展现成一个通用形式

python diff 差分 python求解差分方程_python_16


python diff 差分 python求解差分方程_差分_17

从之前可以得出Δjy0项可以通过递归来求值。为了便于手算,上表已经给出了计算步骤,其中有六个数据点

上表描述的特殊格式称为(牛顿)正向差分,其特征是沿从左上角到右下角向下倾斜的对角线下标保持相同。其他的布局是可能的,如逆向差分和高斯方法,但这些将不在这里讨论。当把表中的数据代入方程,“零”差项,也就是刚开始的差值项,是y值本身,因此

python diff 差分 python求解差分方程_差分_18


计算实例

给与下面的数据是函数y = cosx 的数据,求cons27

首先计算正向差分的数据,写入表格

python diff 差分 python求解差分方程_差分_19


一般方程的插值多项式,x0可以选择任意的初始值x。然而,如果我们希望包括所有五个数据点,生成一个四阶插值多项式, x0应设置等于表中的最高值,也就是说,x0 = 20。x值之间的常量间隔为h = 5,系数可以用以表格形式计算如下

python diff 差分 python求解差分方程_python diff 差分_20


保留小数点后五位,方程的插值多项式可以写成

python diff 差分 python求解差分方程_线性代数_21


扩展得到

python diff 差分 python求解差分方程_线性代数_22


替换需要求解的值,保留到小数点后五位得到

python diff 差分 python求解差分方程_python_23


从上面看出,从上面的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))

终端输出结果如下

python diff 差分 python求解差分方程_多项式_24