第 3 章插值法 实际问题中的函数f(x)是各种各样的,有的数学表达式过于复杂,有的 只给出某些点上的函数值、导数值等离散数据。为了研究或使用f(x), 往往 构造一个简单函数p(x)作为近似函数,通过处理p(x)而获得使用f(x)的 效果。如果只需要p(x)处理或给出离散数据,则称为f(x)的插值函数。 选取不同的插值函数p(x), 近似f(x)的效果就会有所不同。由于代数 多项式结构简单且计算与分析都比较方便,故常用作科研或工程上的插值函 数,这就是代数插值。拉格朗日插值、牛顿插值、埃尔米特插值都可用于寻求 整个区间上的代数插值多项式;但当多项式次数太高(涉及数据点太多)时插 值效果可能会变差,这就需要使用分段拉格朗日插值、分段埃尔米特插值等 方法。如果要求插值函数p(x)是光滑曲线(插值点处的一阶二阶导数连 续), 还可以使用三次样条函数插值多项式。 3.插值及代数插值 1 插值法就是构造函数的近似表达式p(近似地表达表格函数或者复杂函 数f(x)的一种数学方法。如果构造p(时 x, ), ix与f(i) x) 能够满足x处的p(i) x相 等,而在别处则以p(x), 为f(的插值函数。 近似地代替f(则称p( 插值法的第一步是根据待解问题选择恰当的函数类作为原来函数的近 似表达式;第二步是具体构造p(x)表达式。 x) x) x) 3.1.1 插值的概念 假定通过实验观测,得到如表31所示的一批数据(xi,yi),1, 2,…,n,需要从中找到自变量 x 与因变量 y 之间的函数关系,一般可用一个 近似函数y=f(x)来表示。 -i=0, 表3- 1 一般插值数据表 x x0 x1 x2 … xn y y0 y1 y2 … yn 6 　数值计算方法——人工智能、大数据分析的数学基础 依据一批给定的数据点(输入、输出变量的数据)来确定函数,就是确定满足特定要求的 曲线或者曲面。函数y=f(的产生办法常因观测数据及要求的不同而不同, x) 一般可采用 函数插值或数据拟合两种办法来实现。 如果只有一个输入变量、一个输出变量,则属于一元函数的拟合和插值;如果有多个输 入变量,则为多元函数的拟合和插值。 如果要求曲线或曲面通过给定的所有数据点,就是插值问题;如果只是希望反映对象整 体的变化趋势,而不强求曲线或曲面通过所有数据点,这就是数据拟合,又称为曲线拟合或 曲面拟合。 注:在人工智能、大数据分析的数据挖掘过程中,原始数据中往往存在许多不完整或者 偏离的数据。轻则影响执行效率,重则影响执行效果,一般都需要进行数据预处理,较为常 用的是数据插补方法。插值法和拟合是两种主要的数据插补方法。 例3- 1 已知sin(35°10')、sin(35°20')、sin(35°30')三个函数值,如表3-2所示。求另外 两个函数值sin(35°16')和sin(35°27')。 表3- 2 已知点函数值与未知点待求函数值 序号 1 2 3 4 5 x 35°10' 35°16' 35°20' 35°27' 35°30' sinx 0.5759568 ? 0.5783323 ? 0.5807030 解:构造yi=sin(xi)的近似函数(已知点函数值与原函数相同): p(xi)=sin(xi-1) + [sin(xi+1)-sin(xi-1)] × (xi-xi-1) 代入求解两个函数值: si°'0.0.5759568)×0.0. n3516=5759568+(5783323-0.6=5773821 si°'=5783323+(5807030-0.7=0. n35270.0.5783323)×0.5799918 调用Windows计算器求得(精确值): si°'sn(=5773827… n3516=i35+16÷60)0. sin35°27'=in(35+27÷60)=5799922… s0. 考虑到代入p(x)求未知函数值时,使用的已知函数值是近似数,因此,由近似函数 p(x)求得的函数值的精度是比较高的。也就是说,p(x)可以作为sin(x)的插值函数 y=sin(x)≈p(x) 3.1.2 代数插值 由于代数多项式结构简单,计算及理论分析都很方便,故科研与工程上经常选取代数多 项式作为插值函数,这就是代数插值。代数插值方法很多,有拉格朗日插值、分段插值、牛顿 插值、等距结点插值等。 1. 插值多项式 y=x), x1,xn , 对于给定的函数(如表格函数)f(已知n+1 个互异点x0,x2,…,且x0 这些点上的函数取值为y0,y2,…,即 = n xn ì . í .. .. 这样,插值结点的选择便可由选择结构来控制而自动实现了。 注:分段线性插值也有缺点,因其连线为折线,整体的曲线不够光滑。 例3-6 已知一批y=x2 的函数值,如表3-5所示。 表3-5 一批y=x2 的函数值 xi -30 -20 -10 0 10 20 30 40 50 60 70 yi 900 400 100 0 100 400 900 1600 2500 3600 4900 通过分段线性插值法,构造函数y=x2 在区间x∈[-30,70]的近似曲线(多段折线)。 解:本例通过执行Python程序完成任务。程序运行结果如图3-4所示。 附:分段线性插值的Python程序。 第3 章　插值法　 67 图3-4 分段线性插值构造的y=x2 近似曲线 import numpy as np import matplotlib.pyplot as plt def getLine(xn, yn): "分段线性插值闭包" def line(x): index=-1 for i in range(1,len(xn)): "寻找x 所属区间" if x>xn[i]: i=i+1 else: index=i-1 break if index==-1: return -100 '''a0=(x-x_i+1)/(x_i-x_i+1)、a1=(x-x_i)/(x_i+1-x_i) 插值: y=a0*y(i)+a1*y(i+1)''' a0=(x-xn[index+1])/float((xn[index]-xn[index+1])) a1=(x-xn[index])/float((xn[index+1]-xn[index])) return a0*yn[index]+a1*yn[index+1] return line #生成并输出已知函数值列表 xn,yn=[],[] for i in range(-30,80,10): print((i,i**2),end=' ') xn.append(i) yn.append(i**2) #分段线性插值 interpolat=getLine(xn, yn) x=[i for i in range(-30, 70)] 68 　数值计算方法——人工智能、大数据分析的数学基础 y=[interpolat(i) for i in x] #画函数图像 plt.plot(xn, yn, 'ro') plt.plot(x, y, 'b-') plt.show() 3.分段抛物插值 如果有较高的计算精度要求,一般采用分段抛物插值。将抛物插值公式中的x0、x1、x2 换成xi-1、xi、xi+1,则成为 y = (x -xi)(x -xi+1) (xi-1 -xi)(xi-1 -xi+1)yi-1 + (x -xi-1)(x -xi+1) (xi -xi-1)(xi -xi+1)yi + (x -xi-1)(x -xi) (xi+1 -xi-1)(xi+1 -xi)yi+1 这就是分段抛物插值公式。 对于给定的插值点x,究竟选用哪3个结点来进行插值计算呢? 换句话说,该式中下标 i 应取何值,插值效果才最好呢? (1)如果插值点x 位于某两个结点xk-1与xk 之间,则第三个结点可取xk-1,也可取 xk+1,为了提高计算精度,当然应该选靠近x 的那一点。 (2)进一步判断x 究竟偏向区间(xk-1,xk)的哪一侧。 . 当x 靠近xk-1,即|x-xk-1|≤|x-xk|时,在xk-2和xk+1中xk-2靠近x,这时应取 xk-2为第三个插值结点,即令i=k-1; . 反之,当x 靠近xk ,即当|x-xk-1|>|x-xk|时,应取xk+1为第三个插值结点,这 时令i=k。 与分段线性插值相类似,当x≤x1 时,要进行外推,这时令i=1,即取x0、x1、x2 为插 值结点;x>xn-1时,令i=n-1,即取xn-2、xn-1、xn 为插值结点。 按上述分析,下标i 的取值应为 i= 1, x ≤x1 k -1, xk-1 |x -xk | (k =2,3,…,n -1) n -1, x >xn-1 ì . í ... ... 3.4 差商、差分与牛顿插值 应用拉格朗日插值法对y=f(x)插值时,如果阶数不同,则每一项都必须重新计算。 这样,进行高阶插值时,就不能利用低阶插值的结果而要重复计算。牛顿插值法解决了这个 问题。牛顿插值多项式来自差商,其意义在于具有“承袭性”。即在插值过程中,增加的一项 可以从上一项推算出来。 3.4.1 差商与拉格朗日插值公式 差商即均差。k 阶差商可表示为f(x0),f(x1),…,f(xk)的线性组合。一次和二次插 值的拉格朗日公式可以用差商形式表示。为了建立具有承袭性的插值公式,有必要引进差 商的一般概念并研究其基本性质。 第 3 章　插值法　 69 1. 差商的概念与性质 一阶差商是函数值之差与自变量之差之比 xk ]xk -x0 f[x0,= f(xk )-f(x0) 可以作为一阶导数的近似值。二阶差商是一阶差商的差商 xk ]x0,f[x0,xk ]= f[x1,-f[xk ] x1,xk -x0 一般地,有了k-1阶差商,即可递推地定义 k 阶差商 f[x0,x1,…,xk ]f[x1,…,1,xk ]-f[x0,x1,…,1] xk-xk = xk -x0 利用差商的递推定义,可以构造如表3-6所示的差商表来计算差商。 表3- 6 差商表 xi f(xi) 一阶差商二阶差商三阶差商四阶差商 x0 f(x0) x1 f(x1) f[x0,x1] x2 f(x2) f[x1,x2] f[x0,x1,x2] x3 f(x3) f[x2,x3] f[x1,x2,x3] f[x0,x1,x2,x3] x4 f(x4) f[x3,x4] f[x2,x3,x4] f[x1,x2,x3,x4]x0,x1,x2,x3,x4]f[ . . . . . . 下面是将要用到的两个 n 阶差商的基本性质。 ( n 阶差商f[x1,…,]是由函数值f(x0),x1),…,)线性组合而成的 1)x0,xn f(f(xn x1,…,] f[x0,xn = Σ(n) f(xk ) k=0 Π(xk -xj) j=0 j≠ k n (2)差商具有对称性,即在 n 阶差商f[x0,x1,…,xn ]中任意调换两个结点的顺序,其 值不变。例如 f[x0,x1]=f[x1,x0] f[x0,x2]f[x0,=x0,x1] x1,=x1,x2]f[x2, 2. 线性插值公式的差商形式 如果将 p0(x)=f(x0) 当作零次插值多项式,考察线性插值公式 p1(x)f(x0)+ f(x1)-f(x0)(x-x0) = x1-x0 可以理解为:p1(x)是由p0(x)修正得到的。其中修正项(x-x0)的系数 f(x1)-f(x0)c1= x1-x0 0 　数值计算方法——人工智能、大数据分析的数学基础 它实际上是函数增量与自变量增量之比,也就是函数y=f(x)在相应区间(x0,x1)上的平 均变化率,称这为f(x)的一阶差商,并记为 f[x0,x1]= f(x1)-f(x0) x1-x0 于是线性插值公式成为 p1(x)=f(x0) + (x-x0)f[x0,x1] 3. 抛物插值公式的差商形式 可以修正p1(x), 得到抛物插值公式。 由于p2(x)通过点(x0,x0))、(f(x1)) 及(x2,f(x2)), 故可将p2(x)写成 f(x1, p2(x)=p1(x)+c2(x-x0)(x-x1) 也就是 p2(x)=f(x0)+ f(x1)-f(x0)(x-x0)+c2(x-x0)(x-x1) x1-x0 该式中,因为 c2 是修正项的系数。显然 , p2(x0)=f(x0),p2(x1)=f(x1),p2(x2)=f(x2) 并将x=x2 代入上式,可得 p2(x2)f(x0)+ f(x1)-f(x0)(-x0)+c2(-x0)(-x1)f(x2) =x1x2x2= x1-x0 改写为 f(x2)-f(x0) f(x1)-f(x0)- x2-x0 x1-x0 c2= x2-x1 这个系数实际上是一阶差商的差商,用 x1] x2-x1 表示,称为二阶差商。 f[x0,x1,x2]= f[x0,x2]-f[x0, 这样,抛物插值公式就可以表示为 x0) + (x0,-x0)(x0,x2] p2(x)=f(x-x0)f[x1] + ( x x-x1)f[x1, 3.4.2 牛顿插值 根据差商定义,可将一次与二次差商形式的拉格朗日插值公式推广到n+1 个插值结 点的情形,同时还可得到插值多项式的余项。 1. 牛顿插值公式 i=1,n)之外, i=1,2,…, 在n+1 个插值结点xi(2,…,再给一个结点x≠xi(n), 则由 差商定义得 x)x, f(=f(x0) + (x-x0)f[x0] f[x,=f[x1] + (xf[x0, x0]x0,-x1)x,x1] x0,f[x1,] + (-x2)x,x1, . f[x,x1]=x0,x2xf[x0,x2] 第 3 章　插值法　 71 f[x,x1,…,f[x1,…,xf[x0,] x0,1]x0,] + ()x,x1,…, 将其中第二式f[x,x0]代入第一式f(x), 得到 f(x)x0) + (-x0)f[x0,x-x0)(xf[x0, xn-=xn -xn xn =f(xx1] + (-x1)x,x1] 也就是 f(x)=p1(x)+R1(x) 其中,p1(x)就是差商形式的一次(线性)插值多项式,而R1(x)为线性插值余项,可表 示为 R1(x)f(-p1(=(-x0)(xf[x0,x1] =x)x) x -x1)x, 类似地,将第三式代入,则 有 f(x)=f(x0) + (x-x0)f[x0,x1] + (-x0)(-x1)f[x1, x0,x2] + )((x) -x1)(-x2)x,x1, (x(x) -x0x xf[x0,x2] 也就是 f(x)=p2(x)+R2(x) 其中,p2(x)就是差商形式的二次(抛物)插值多项式,而R2(x)为线性插值余项,可表示为 R2(x)=x)-p2(=(-x0)(-x1)(-x2)x,x1, f(x) x x xf[x0,x2] 至此,可以推断出一般规律:每增加一个插值点,只要将高一阶差商代入其前一公式。 以此类推,即可得到 =x0) + (-x0)x0,f(x)f(xf[x1] + (-x0)(-x1)f[x1, x0,x2] + … + f[x1,…, (-x0)(-x1)…(-xn )x,x1,…, 记该式中最后一项为Rn (x(x) ), 即 x xf[x0,xn ] Rn (=(-x0)(-x1)…(-xn f[x0,xn (x(x) -x0)(x(x) -x1)…(x-xn-1)x0,xn ] + x) x )x,x1,…,] 再记该式中前n+1 项为pn (x), 即((x) 差商形式的插(x) 值公式) x)f(x0, pn (=x0) + (-x0)f[x1] + f[x1, (x-x0)(x(x) -x1)x0,x2] + … + (-x0)(-x1)…(x1,…,] x x x -xn-1)f[x0,xn 则有 f(x)=pn (x)+Rn (x) 注意到插值余项Rn (i=1,n), ( xi)=0(2,…,故构造而成 n 次多项pn x)式必然 满足 (xi)f(xi),i=1,2,…, 可见,pn (x)就是符合要求的插值多项式。这种差商形式的插值公式称为牛顿插值公式。 pn= n 2. 牛顿插值的计算 牛顿插值公式具有承袭性,其优点是计算工作量小,特别是当提高插值阶数时,拉格朗 日插值公式中每项都必须重新计算,而牛顿插值公式中却只需计算最后一项,将此项值累加 到前一阶插值结果上即可。牛顿插值公式可用于求非等距结点的函数值。 应用牛顿插值法,大体上按以下步骤操作。 72 　数值计算方法——人工智能、大数据分析的数学基础 S1 计算各阶差商值,即构造差商表。 S2 按牛顿插值公式计算插值结果。 S3 比较插值结果,判断:需要高一阶插值吗? 是则转S2。 S4 输出插值结果。 S5 结束。 例3-7 列表函数f(x)=lg(x)的值如表3-7所示,用牛顿插值法求lg4.01的值。 表3-7 lgx 函数表 x 4.0002 4.0104 4.0233 4.0294 f(x) 0.6020817 0.6031877 0.6045824 0.6052404 解:已知x0=4.0002,x1=4.0104,x2=4.0233,x3=4.0294。 根据给定的函数表作差商表,如表3-8所示。 表3-8 求lg4.01的差商表 x f(x) 一阶差商二阶差商三阶差商 4.0002 0.6020817 4.0104 0.6031877 0.108431 4.0233 0.6045824 0.108116 -0.013636 4.0294 0.6052404 0.107869 -0.013000 0.021781 代入牛顿插值公式,求得 lg4.01=f(x0)+ (x -x0)f[x0,x1]+ (x -x0)(x -x1)f[x0,x1,x2]+ (x -x0)(x -x1)(x -x2)f[x0,x1,x2,x3] =0.6020817+ (4.01-4.0002)×0.108431+ (4.01-4.0002)× (4.01-4.0104)× (-0.013636)+ (4.01-4.0002)× (4.01-4.0104)× (4.01-4.0294)×0.021781 =0.6031443 注:插值多项式pn (x)的系数就是差商表斜线上的各阶差商。如果需要更高阶的差 商,则可按差商定义进行计算,将表向下向右延伸,当给定x 而要计算pn(x)的值时,可将x -xi 的值列在表右端,更便于计算。 附 牛顿插值的Python程序。 import numpy as np def newton(xi,fi,x): #牛顿插值: 计算各阶差商 n=len(xi) c=np.zeros((n,n)) for i in range(n): 第3 章　插值法　 73 c[i,0]=fi[i] for i in range(1,n): for j in range(i,n): c[j,i]=(c[j,i-1]-c[j-1,i-1])/(xi[j]-xi[j-i]) #牛顿插值: 计算插值结果 s=fi[0] for i in range(1,n): t=1 for j in range(0,i): t*=(x-xi[j]) s+=c[i,i]*t return s if _ _name _ _=='_ _main_ _': #已知结点及其函数值 xKnown=np.array([4.0002, 4.0104, 4.0233, 4.0294]) fKnown=np.array([0.6020817, 0.6031877, 0.6045824, 0.6052404]) #新结点求值 xNew=4.01 fNew=newton(xKnown, fKnown, xNew) print(fNew) 程序的运行结果如下: 0.6031443812538274 3.牛顿插值公式的余项 对于差商表示的插值余项 Rn(x)=(x -x0)(x -x1)…(x -xn)f[x,x0,x1,…,xn] 可通过差商与导数的关系,推导出用导数表示的插值余项。 根据微分中值定理,在(x0,x1)区间内存在一点ξ,使得 f'(ξ)= f(x1)-f(x0) x1 -x0 =f[x0,x1] 从该式得到差商与一阶导数的关系,再推广到n 阶导数,则有 f[x0,x1,…,xn]=f(n)(ξ) n! 其中,ξ 为插值区间内一点。这个公式的正确性可用洛尔定理证明。 类似地,只要再增加一个点x≠xi(i=1,2,…,n),即可将插值余项Rn (x)中的n+1 阶差商,用n+1阶导数表示为 f[x,x0,x1,…,xn]=f(n+1)(ξ) (n +1)! 其中,ξ、x 均为插值区间内的点,且ξ 依赖于x,即当x 变化时ξ 也随之变化。将此结果代 入差商表示的插值余项式,即可得到用导数表示的插值余项 Rn(x)=(x -x0)(x -x1)…(x -xn)f(n+1)(ξ) (n +1)! 改写为 74　数值计算方法——人工智能、大数据分析的数学基础 ξ) (x)(x-xi)fn( +1)! Rn = Π(n) ( n+1)( = 这就是牛顿插值的余项公式,与拉格朗日(i) 插(0) 值余项公式完全一样。 由插值多项式的唯一性可知,牛顿插值多项式与拉格朗日插值多项式是相等的,它们的 余项也是相等的。故可得到这样的一个等式关系。 与拉格朗日插值多项式相比,牛顿插值多项式更具有一般性。其误差项对于仅由离散 点给出的f(x)或者导数不存在的f(x)仍然适用,从而应用更为广泛。 3.4.3 差分、差商及导数的关系 在实际问题中,列表函数的自变量往往按等距离的点给出。这种情况下,函数插值就要 用差分而不是差商了。也就是说,当所有插值结点都是等距离时,可将插值公式表示为简单 的差分形式。 1. 向前差分 设函数f(x)在等距结点x=h(=2,…,上的值为y0=xi), x)在每个 x0+ii1,n) f(则f( xi+1] i+1-yi, x) 小区间[xi,上的改变量y称为f(在点xi 上的一阶差分或向前差分。记作 Δyi=yi+1-yi 对一阶差分再取一次差分就是二阶差分,记为 Δ2yi=Δyi+1-Δyi =(yi+2-(-yi) -yi+1)yi+1=yi+2-2yi+1+ y 一般地, n 阶差分为n-1阶差分的差分,即 Δnyi=Δn-1yi+1-Δn-1yi 依据差分逐层递推的特点,计算向前差分时,用如表3-9所示向前差分表(给出四阶差 分)较为方便。 表3- 9 向前差分表 xi yi Δyi Δ2yi Δ3yi Δ4yi x0 y0 Δy0 Δ2y0 Δ3y0 Δ4y0 x1 y1 Δy1 Δ2y1 Δ3y1 x2 y2 Δy2 Δ2y2 x3 y3 Δy3 x4 y4 Δ 2. 向后差分 与向前差分相仿,对于函数f(x), 将每个等长小区间[xi,xi+1]上的改变量f(xi) 1)称为f(x) 处的一阶向后差分。记作 f(xi 在点xi - yi=yi-yi-1 ΔΔ 类似地,二阶差分为一阶差分的差分 2 yi=yi Δ yi-1