实验名称插值法

1、探索实验5 插值法一、 实验目的 了解插值问题及其适用的场合，理解并掌握常用的插值算法的构造和计算，了解差商概念、Runge现象及样条插值方法，学习用计算机求近似函数的一些科学计算方法和简单的编程技术。 二、概念与结论1. 插值问题与插值函数： 由实验或测量的方法得到所求函数 y=f(x) 在互异点x0 , x1, . , xn 处的值 y0 , y1 , , yn ,构造一个简单函数 j(x) 作为函数 y=f(x) 的近似表达式 y= f(x) » j(x)使 j(x0)=y0 , j(x1)=y1 , ¼, j(xn)=yn ,(1)这类问题称为插值问题。 f(x)

2、称为被插值函数，j(x) 称为插值函数， x0 , x1, . , xn 称为插值节点。(1)式称为插值条件。常用的插值函数是多项式函数。且当n=1时是称为线性插值，n=2时称为Simpson插值或抛物线插值。2插值定理： 假设x0 ,x1,xn 是n+1个互异节点,函数f(x)在这组节点的值f(xk)(k=0,1,n)是给定的，那么存在唯一的n 次次多项式pn (x)满足pn (xk)=f(xk), k=0,1,n3.插值的截断误差 设jn(x)是过点x0 ，x1 ，x2 ，xn的 n 次插值多项式， f(n+1)(x)在（a，b）上存在，其中a，b是包含点x0 ，x1 ，x2 ，,xn的任

3、一区间，则对任意给定的xÎa，b，总存在一点xÎ（a，b）（依赖于x）使 其中wn+1(x)=(x x0) (x - x1)(x-xn ) ，f(n+1)(x) 是f(x)的n+1阶微商在 x 的值。 4. 差商：给定一个函数表 x | x0 x1 . xn - - y | y0 ,y1 . yn其中当i ¹j时 ,xi¹xj 记 fxi=f(xi) ,定义f(x)关于xi,xj的一 阶差商 一般的, f(x)关于xi,xi+1,xi+k的k 阶差商定义为:5. 分段线性插值： 在区间a,b上给定一组节点: a=x0<x1<x2<<

4、;xn=b且在节点上的函数值f(xi )=yi, i=0,1,2,n, 如果函数j(x) 满足：1) j (x)在a , b上连续2) j (xi )=yi, i=0,1,2,n3) j (x)在每个子区间xi , xi+1(i=0,1,2,n-1)上是次数为1的多项式则称j (x)是定义在a ,b上的分段线性插值函数。6. m次样条插值函数： 设f(x)是区间a,b上的一个连续可微函数,在区间a,b上给定一组节点: a=x0<x1<x2<<xn=b且y = f(x)在点 x0，x1，x2, xn的值为y0，y1，y2, yn，如果函数s(x)满足条件 (1) s(x)

5、在每个子区间xi , xi+1(i=0,1,2,n-1)上是次数不超过m的多项式;(2) s(x)在区间a , b上有m-1阶连续导数;(3) s(xi)=f(xi) =yi , i=0,1,2,n则称s(x) 为函数f(x)的在a ,b上的m次样条插值函数。x0，x1，x2, 称为样条结点,其中x1,xn-1称为内结点, x0 , xn 称为边界结点。当m=3时,则称s(x)为函数f(x)的三次样条插值函数。三、程序中Mathematica语句解释：1. Product f(i) , i ，imin，imax，h 计算乘积 f(imin)´f(imin +h) ´f(im

6、in +2h) ´´f(imin +nh) imax h £ imin + nh £ imax , h>02Simplifyexpr 对expr进行化简3. Appendlist, elem 在表list的后面加一个元素 elem 4. Prependlist, elem 在表list的前面加一个元素 elem四、方法与程序 插值法是函数逼近的重要方法之一, 它是求近似函数的一种方法，有着广泛的应用。插值法有很多种,其中以拉格朗日(Lagrange)插值和牛顿(Newton)插值为代表的多项式插值最有特点,常用的插值还有Hermit插值,分段插值和

7、样条插值等。这里只给出Lagrange插值、Newton插值 、分段线性插值和样条插值的构造过程及程序。1Lagrange插值 Lagrange插值是将待求的n次多项式插值函数Pn(x)改写成另一种表示方式，再利用插值条件确定其中的待定插值基函数，从而求出插值多项式。Lagrange插值是多项式插值，它成功地用构造插值基函数的方法解决了求多项式插值函数出现的病态问题。1.1 Lagrange插值的构造过程： 令插值多项式为 Ln(x)=l0n(x)y0 + l1n(x)y1 + l2n(x)y2 + lnn(x)yn 式中 l0n(x), l1n(x), l2n(x), ,lnn(x) 都是次

8、数£ n 的待定多项式，利用插值条件 Pn(x0)=y0 , Pn(x1)=y1 , . , Pn(xn)=yn可以得到 l0n(x0) = 1, l1n(x0) = 0, l2n(x0) = 0, , lnn(x0) = 0 l0n(x1) = 0, l1n(x1) = 1, l2n(x1) = 0 , , lnn(x1) = 0 . l0n(xn) = 0, l1n(xn) = 0, l2n(xn) = 0, , lnn(xn) =1考虑函数 l0n(x),由其在 x1 , x2 . xn 的零值，可以知道它含有 x-x1 , x-x2 ,. ,x-xn n个因子，于是可令 l0

9、n(x) = a(x-x1)(x-x2). (x-xn)利用 l0n(x0)=1 确定其中的系数a，得 (x-x1)(x-x2)(x-xn )l0n(x) = - (x0-x1)(x0-x2 )(x0-xn ) 类似的可以得出 l1n(x) , l2n(x), lnn(x) (x-x0)(x-x1) (x-xi-1 ) (x-xi+1) (x-xn )lin(x)= - (xi-x0)(xi-x1 )(xi-xni-1 ) (xi-xi+1) (xi-xn ) l0n(x), l1n(x) , lnn(x)称为以x0 , x1 ,. ,xn为节点的插值基函数。于是，将这些插值基函数代入前面的多

10、项式Ln(x)即得到如下所求的n次Lagrange插值多项式： 1.2 求n次Lagrange插值多项式算法： 1 输入n+1个插值点： (xi, yi),i=0,1,n2 计算插值基函数l0n(x), l1n(x) , lnn(x)3 给出n次Lagrange插值多项式：Ln(x)= y0 l0n(x)+ y1 l1n(x) +yn lnn(x)1.3 求Lagrange插值多项式程序：Clearlag,xi,x,yi;xi=Input"xi="yi=Input"yi="n=Lengthxi-1;p=Sumyii*(Product(x-xij)/(xi

11、i-xij),j,1,i-1*Product(x-xij)/(xii-xij),j,i+1,n+1),i,1,n+1;lagx_=Simplifyp说明：本程序用于求n次Lagrange插值多项式。程序执行后，按要求通过键盘输入插值基点xi:x0 , x1, . , xn 和对应函数值yi: y0 , y1 , , yn 后，程序即可给出对应的n次Lagrange插值多项式lagx。程序中变量说明：xi：存放插值基点x0 , x1, . , xn yi: 存放对应函数值y0 , y1 , , ynlagx: 存放求出的n次Lagrange插值多项式Ln(x)注：语句lagx_=Simplify

12、p用简化形式给出对应的n次Lagrange插值多项式。 1.4 例题与实验 例1.给定数据表x 0 1 2 3y=f(x) 1 3 5 12用Lagrange插值法求三次插值多项式,并给出函数f(x)在x =1.4的近似值。解： 执行Lagrange插值程序后，在输入的两个窗口中按提示分别输入0, 1, 2, 3、1, 3, 5, 12，每次输入后用鼠标点击窗口的“OK”按扭，得如下插值函数： 6 + 22 x - 15 x2 + 5 x3- 6所以得到三次插值多项式L3(x)=1+11 x/3-5 x2/2+5 x3/6接着键入： lag1.4得输出:3.52，因此f(x)在x =1.4的近

13、似值为3.52，即f(1.4)»3.52.2. Newton插值 Newton插值是多项式插值的另一种表示形式，它在增加插值节点时具有灵活性，只要在原有的插值多项式中增加一项就得到新的插值函数。2.1 Newton插值构造过程： 为了使插值Newton多项式具有承袭性，令插值函数具有下列形式：式中称为Newton插值基函数。为求出Nn(x),利用插值条件和差商，得出如下n次Newton插值多项式：计算 Newton插值多项式Nn(x)常用如下差商表：2.2 求Newton插值多项式算法： 1. 输入n+1个插值点： (xi, yi),i=0,1,n2. 计算差商表3.给出n次Newt

14、on插值多项式。2.3 求Newton插值多项式程序：Clearnewt,s,x;xi=Input"xi="yi=Input"yi="n=Lengthxi;(*计算差商表*)f=Table0,n,n;Dofi,1=yii,i,1,nDofi,j+1=(fi,j-fi+1,j)/(xii-xii+j), j,1,n-1,i,1,n-jPrint"差商表"DoPrintxii," ",fi,i,1,n(*求Newton插值多项式*)fa=1;s=f1,1;Dofa=(x-xik)*fa;s=s+fa*f1,k+1,k,

15、1,n-1newtx_=sSimplify%说明：本程序用于求n次Newton插值多项式。程序执行后，按要求通过键盘输入插值基点xi:x0 , x1, . , xn 和对应函数值yi: y0 , y1 , , yn 后，程序依次给出输入的数据表、计算出的差商表、Newton插值多项式、Newton插值多项式的简化形式。程序中变量说明：xi：存放插值基点x0 , x1, . , xn yi: 存放对应函数值y0 , y1 , , ynf:存放函数值y0 , y1 , , yn及所有差商newtx: 存放求出的n次newton插值多项式Nn(x)注：1)语句f=Table0,n,n用于产生一个n&

16、#180;n的矩阵变元用于存放函数值y0 , y1 , , yn及所有差商。2)在Mathematica中有一个求n次插值多项式的命令，命令形式为： InterpolatingPolynomialx01,y0，x1,y1,x2,y2,xn,yn, x 它可以求过n+1个插值点x01,y0，x1,y1,x2,y2,xn,yn的n次插值多项式Pn(x)。2.4例题与实验例2.给定数据表 x 4.0002 4.0104 4.0233 4.0294 y 0.6020817 0.6031877 0.6045824 0.60524041) 计算差商表2) 用Newton插值法求三次插值多项式Nn(x)3)

17、 求f(4.011)的近似值解： 执行Newton插值程序后，在输入的两个窗口中按提示分别输入4.0002, 4.0104, 4.0233, 4.0294、0.6020817, 0.6031877, 0.6045824, 0.6052404,每次输入后用鼠标点击窗口的“OK”按扭，得如下插值函数：4.0002, 4.0104, 4.0233, 4.02940.6020817, 0.6031877, 0.6045824, 0.6052404 差商表4.0002 0.6020817, 0.108431, -0.0136404, 0.02116294.0104 0.6031877, 0.108116

18、, -0.0130225, 04.0233 0.6045824, 0.107869, 0, 04.0294 0.6052404, 0, 0, 00.6020817 + 0.108431 (-4.0002 + x) - 0.0136404 (-4.0104 + x) (-4.0002 + x) + 0.0211629 (-4.0233 + x) (-4.0104 + x) (-4.0002 + x)-1.41642 + 1.23926 x - 0.268313 x2 + 0.0211629 x3于是我们得到本题的差商表：xi yi f, f , , f , , ,4.0002 0.6020817

19、, 0.108431, -0.0136404, 0.02116294.0104 0.6031877, 0.108116, -0.01302254.0233 0.6045824, 0.107869 4.0294 0.6052404 和三次插值多项式N3(x)= -1.41642 + 1.23926 x - 0.268313 x2 + 0.0211629 x3接着键入： newt4.011得输出0.603253，因此f(4.011)» 0.603253例3. 多项式插值的误差估计式中可以看到，当插值节点越多时误差会越小，这个结论正确吗？通过实验说明该结论的正确性。 解: 考虑函数f(x)

20、 = (1+x2)-1 在区间-4,4内选取不同个数的等距插值节点做观察，这里分别选-4,4 内的9个和11个的等距节点来做实验，将对应的插值函数图与被插函数f(x) = (1+x2)-1画在一起做观察，为简单起见，这里用Mathematica 命令做实验，对应命令为： In1:= u=Tablex,(1+x2)-1,x,-4,4 ; （*采取f(x) 在-4,4 内的9个插值点 In2:= g=ListPlotu, PlotStyle->PointSize0.04 （*将散点图图形文件存放在变量g中In3:= s=InterpolatingPolynomialu , x ; （*将插值

21、函数存放在变量s中In4:= t= Plots, (1+x2)-1, x,-4,4, PlotStyle->Thickness0.005,Thickness0.006 （*将插值函数s与f(x)画在一起的图形文件存放在变量t中In5:= Showt, g （*将散点图, 插值函数s, f(x)画在一个坐标系中图 1 在-4,4中选9个等距节点的插值函数与被插函数图，粗线为被插函数图In6:= u=Tablex,(1+x2)-1,x,-4,4，0.8 ; （*采取f(x) 在-4,4 内的11插值点 In7:= g=ListPlotu, PlotStyle->PointSize0.0

22、4 In8:= s=InterpolatingPolynomialu , x ; In9:= t= Plots, (1+x2)-1, x,-4,4, PlotStyle->Thickness0.005,Thickness0.006 In10:= Showt, g 图 2 在-4,4中选11个等距节点的插值函数与被插函数图 图中细的曲线是插值函数, 粗的曲线是原来函数(被插函数)。从图形上观察, 说明增加插值节点不能保证使插值误差变小，本题图象可以看到，节点的增多得到的插值函数产生的误差在靠近端点的函数值不但没有减少反而增大了。因此当插值节点越多时误差会越小，这个结论不正确。这个实验也说明

23、高次插值函数效果不好。插值点个数增加，出现误差增大的现象称为Runge现象。Runge现象的发现使人们知道了高次插值函数误差较大, 这也是为什么实用中一般不使用高次多项式插值的原因。3. 分段线性插值 分段线性插值是最简单的分段插值，它在每个子区间xi , xi+1(i=0,1,2,n-1)上用直线段连接数据点构成的折线函数作为插值函数的。分段线性插值函数具有计算简单，且有每个插值子区间趋于零时，分段线性插值函数可以克服Runge现象，且具有随着插值节点的增多，它可以是插值误差不断减小、一致逼近被插函数的特点。 3.1 分段线性插值函数的构造过程： 在每个子区间xi , xi+1(i=0,1,

24、2,n-1)上，用对应的两点：x xi xi+1y yi yi+1做过此两点的直线： L(x)=k(x -xi)+ yi这里k=(yi+1 -yi)/(xi+1 -xi)，然后将所有直线连起来即可。显然这样得到的折线满足插值条件。 分段线性插值函数是一个分段函数。3.2 用分段线性插值函数做近似计算的算法：1. 输入n个插值点： (xi, yi),i=1,n2. 输入要近似计算的自变量点xa3.寻找包含xa的小区间xi , xi+14.用 xi , xi+1 上的线性插值函数在xa处的值作为f(xa)的函数值3.3 用分段线性插值函数做近似计算的程序：Clearx,a,b;lia_,b_,x_

25、:=(b2-a2)/(b1-a1)*(x-a1)+a2xi=Input"xi="yi=Input"yi="xa=Input"xa="n=Lengthxi;Ifxa<xi1|xa>xin,Print"超限"Break, DoIfxa>=xik,m=k,k,1,n-1; Print"m=",m,"x=",xim; lixim,yim,xim+1,yim+1,xa说明：本程序用分段线性插值做近似计算。程序执行后，按要求通过键盘输入插值基点xi: x1, . , x

26、n 和对应函数值yi: y1 , , yn 和要做近似计算的点xa后，程序将给出包含xa小区间xi , xi+1的下标i和xi ,和函数f(xa)的近似值。如果xa不在x1, xn 内程序给出超限提示。程序中变量说明：xi：存放插值基点 x1, . , xn yi: 存放对应函数值 y1 , , ynxa: 存放要做近似计算的自变量值m:存放包含xa小区间xi , xi+1的下标i注：在Mathematica中有一个求分段线性插值函数的命令，命令形式为： Interpolation x1,y1,x2,y2,.,xn,yn，InterpolationOrder -> 1用如上Mathema

27、tica命令求出的分段线性插值函数没有给出具体的分段函数表达式, 而是用“InterpolatingFunctionx1, xn, <> ”作为所求的分段插值函数, 通常可以用 变量=Interpolation 数据点集合把所得的分段插值函数存放在变量中, 如果要计算插值函数在某一点的如p点的值,只要输入 “变量p” 即可，如果要表示这个插值函数应该用 “变量x” 。此外，命令 Interpolation x1,y1,x2,y2,.,xn,yn，InterpolationOrder -> 3可以给出更好的分段插值函数。 3.4例题与实验例4.给定数据表xi 1 2 3 4 y

28、i 0 -5 -6 31） 用分段线性插值函数求函数f(x)在x =1.4的近似值2） 用Mathematica命令画出分段线性插值图形解：执行求分段线性插值函数程序后，在输入的两个窗口中按提示分别输入1，2，3，4、0，-5，-6，3、1.4每次输入后用鼠标点击窗口的“OK”按扭，得如下插值函数：m=1 x=1-2此结果说明f(x)在x =1.4的近似值为2，它对应的小区间为x1 , x2。为完成第二个问题，键入命令：d=Interpolation 1,0,2,-5,3,-6,4,3,InterpolationOrder -> 1Plotdx,x,1,4得输出如下：Interpolat

29、ingFunction1, 4, <>图3例4的分段线性插值函数图形。例5.在例3中，用分段插值函数代替n次多项式插值函数，观察插值节点增多时误差的变化情况。 解: 对函数f(x) = (1+x2)-1 在区间-4,4内选取不同个数的等距插值节点做观察，这里分别选-4,4 内的9个和21个的等距节点来做实验，将对应的插值函数图与被插函数f(x) = (1+x2)-1画在一起做观察，为简单起见，这里用Mathematica 命令做实验，对应命令为： In1:= u=Tablex,(1+x2)-1,x,-4,4 ; （*采取f(x) 在-4,4 内的9个插值点 In2:= g=List

30、Plotu, PlotStyle->PointSize0.04 （*将散点图图形文件存放在变量g中In3:= s= Interpolationu,InterpolationOrder -> 1 ; （*将线性插值函数存放在变量s中In4:=t=Plotsx, (1+x2)-1, x,-4,4,PlotStyle->Thickness0.005,Thickness0.006 （*将线性插值函数s(x)与f(x)画在一起的图形文件存放在变量t中In5:= Showt, g （*将散点图, 线性插值函数s(x), f(x)画在一个坐标系中图 4 在-4,4中选9个等距节点的线性插值

31、函数与被插函数图，粗线为被插函数图In6:= u=Tablex,(1+x2)-1,x,-4,4，0.4 ; （*采取f(x) 在-4,4 内的21插值点 In7:= g=ListPlotu, PlotStyle->PointSize0.04 In8:= s= Interpolationu,InterpolationOrder -> 1 In9:=t=Plotsx, (1+x2)-1, x,-4,4,PlotStyle->Thickness0.005,Thickness0.006 In10:= Showt, g （*将散点图, 线性插值函数s(x), f(x)画在一个坐标系中图

32、 5 在-4,4中选21个等距节点的插值函数与被插函数图，粗线为被插函数图 从图形上观察, 说明增加插值节点能保证使插值误差变小。因此当插值节点越多时分段线性插值误差会越小。4样条插值 分段插值会导致插值函数在子区间的端点(衔接处)不光滑,即导数不连续,对于一些实际问题,不但要求一阶导数连续,而且要求二阶导数连续。三次样条插值可以满足这些要求的概念。 所谓“样条”(SPLINE)是工程绘图中的一种工具,它是有弹性的细长木条,绘图时,用细木条连接相近的几个结点,然后再进行拼接,连接全部结点,使之成为一条光滑曲线,且在结点处具有连续的曲率。样条函数是对这样的曲线进行数学模拟得到的。它除了要求给出各

33、个结点处的函数值外,只需提供两个边界点处导数信息,便可满足对光滑性的不同要求项式。这里给出最常用的三次样条插值的构造与算法。4.1 三次样条插值函数的构造过程： 设三次样条插值函数为S(x)。记Mi = S(xi), f(xi)= fi= yi ,考虑它在任一区间xi,xi+1上的形式.根据三次样条的定义可知 ,S(x)的二阶导数S(x)在每一个子区间xi,xi+1 ( i=0,1,2,n-1)上都是线性函数.于是在xi,xi+1 上S(x)=Si(x)的二阶导数表示成其中 hi= xi+1xi . 对S(x)连续积分两次,并利用插值条件S(xi)= yi ,得到因此，只要能求出所有的M i，

34、就能求出三次样条插值函数S(x)。 为求出Mi，利用三次样条插值函数S(x)在节点处的连续性S'(xi-)= S'(xi+) ,(i=1,2,n-1)可得iMi-1+2Mi+iMi+1= dI (i=1,2,n-1) 其中 上面的方程组有n-1个方程，但有n+1个变量Mi，故需两个方程才能求唯一解，为此引入下列常用的边界条件。第一型边界条件： 已知f(x)在两端点的导数f(a)和f(b) ，要求 S(a) = f(a) , S(b) = f(b)第二型边界条件： 已知f(x)在两端点的二阶导数f"(a)和f"(b) ，要求 S(a)=M0 = f(a) ,

35、S(b)=Mn= f(b)特别当 S(a)= S(b) =0时，S(x)称为自然三次样条第三型边界条件： 已知f(x)是以b -a为周期的周期函数 ，要求S(x)满足周期条件S (a) = S(b) , S(a+)= S(b-) , S(a+)= S(b-)三次样条插值问题加上第i型边界条件称为第i型插值问题（i，）。可以证明第i型插值问题的解是存在且唯一的。对第一型边界条件和第二型边界条件，它们的参数 Mk 对应如下的三对角方程组：其中，在第一型边界条件下：l0=1, d0=6(fx0,x1-f¢(x0)/h0mn=1, dn=6(f¢(xn)-fxn-1,xn) /hn

36、-1在第二型边界条件下：l0=0, d0=2f¢¢(x0)mn=0, dn=2f¢¢(xn)以上各组条件与方程组（*）联立，可以解出未知参数M0，M1 ,Mn，然后代入S(x) 表达式，即可求得样条函数 。上面构造方法中Mi相应于力学中细梁在xi处截面的弯矩，每一个方程中又至多出现相邻的三个Mi，通常称为三弯矩法。 4.2 求三次样条插值函数算法：1.确定边界条件，判定是第几型插值问题；2.根据所确定的条件计算各值，形成方程组(*)；3.解三对角方程组(*)，求得M0， M1 ， M2, Mn ；4.将求得的Mi值代回S(x)的表达式中4.3 求三次样条

37、插值函数程序：mar=Input"边界条件类型=1或2(1:第一型,2：第二型)"f1=Input"输入左边界条件函数值="f2=Input"输入右边界条件函数值="xi=Input"输入插值节点xi="yi=Input"输入对应函数值yi="n=Lengthxi-1;h=Tablexii+1-xii,i,1,n;mu=Tablehi/(hi+hi+1),i,1,n-1;lam=1-mu;d=Table 6(yii+1-yii)/hi-(yii-yii-1)/hi-1)/(hi-1+hi), i

38、,2,n;Ifmar=2, d=Prependd,2f1; d=Appendd,2f2; mu=AppendPrependmu,0,0; lam=AppendPrependlam,0,0, d=Prependd,6*(yi2-yi1)/h1-f1)/h1; d=Appendd,6*(f2-(yin+1-yin)/hn); mu=AppendPrependmu,0,1; lam=AppendPrependlam,1,0c=lam;b=Table2,n+1;a=mu;n=n+1;eps1=0.000000001;Do ak=ak/bk-1;bk=bk-ak*ck-1; dk=dk-ak*dk-1,

39、 k,2,n;x=Table0,n;xn=dn/bn;Doxk=(dk-ck*xk+1)/bk, k,n-1,1,-1;Print"Mk的值为 " , x/N;Clearx1,x2,y1,y2,m1,m2,h1,t,skskx1_,x2_,y1_,y2_,m1_,m2_,h1_,t_:=m1*(x2-t)3/(6*h1)+m2*(t-x1)3/(6*h1)+(y1- m1*h12/6)*(x2-t)/h1+(y2-m2*h12/6)*(t-x1)/h1;sp1=Tableskxii,xii+1,yii,yii+1,xi,xi+1,hi,t, i,1,n-1;Print&qu

40、ot;所求样条函数S(t)为"DoPrint"S(t)=",sp1i; Print" ",xii,"<=t<",xii+1, i,1,n-1说明：本程序用于求具有第一型边界条件和第二型边界条件的三次样条插值函数。程序执行后，按要求通过键盘输入：边界条件类型=1或2(1:第一型,2：第二型)、左边界条件函数值f¢(x0)或f¢¢(x0)、右边界条件函数值f¢(xn)或f¢¢(xn)、插值基点xi:x0 , x1, . , xn 和对应函数值yi: y0 ,

41、 y1 , , yn 后，程序即可给出M0， M1 ， M2, Mn和三次样条插值函数。程序中变量说明：xi：存放插值基点 x0, x1, . , xn yi: 存放对应函数值 y0 , y1, , ynh: 存放小区间长度hk mu: 存放数值mk lam 存放数值lk d: 存放数值dkx：存放求出的Mk值M0 , M1, . ,Mn sp1: 存放求出的三次样条插值函数s(x)skx1,x2,y1,y2,m1,m2,h1,t:在区间x1,x2上的三次样条插值函数s(t)注：1）语句sp1=Tableskxii,xii+1,yii,yii+1,xi,xi+1,hi,t, i,1,n-1是以

42、表的形式给出三次样条插值函数s(t)的，其中在sp1中，第i个分量表示三次样条插值函数s(t)的第i段.4.4例题与实验例6. 已知函数f(x）的取值如下表： x 27.7 28 29 30f(x) 4.1 4.3 4.1 3.0 f(x) 3.0 -4.0试求f(x) 在27.7,30上的三次样条插值函数并计算f(28.4)的近似值。解：这是第一类边界条件的问题，执行三次样条插值程序后，在输入的五个窗口中按提示分别输入：1、3.0 、4.0、27.7, 28, 29, 30、4.1, 4.3, 4.1, 3，每次输入后用鼠标点击窗口的“OK”按扭，得如下输出结果： Mk的值为 -23.531

43、4, 0.39604, 0.829703, -9.11485所求样条函数S(t)为S(t)=14.8432 (28 - t) - 13.073 (28 - t)3 + 14.3135 (-27.7 + t) + 0.220022 (-27.7 + t)3 27.7£t<28S(t)=4.23399 (29 - t) + 0.0660066 (29 - t)3 + 3.96172 (-28 + t) + 0.138284 (-28 + t)3 28£t<29S(t)=3.96172 (30 - t) + 0.138284 (30 - t)3 + 4.51914 (

44、-29 + t) -1.51914 (-29 + t)3 29£t<30因为28.4Î28,29接着键入： sp12/.t->28.4得输出: 4.14819，因此f(28.4)» 4.14819例7. 已知函数f(x）的取值如下表： x 0 1 2 3f(x) 0 0.5 2 1.5 试求1) 给定边界条件f¢(0)=0.2，f¢(3)=-1，求相应的三次样条插值函数s(x)并画出图形2) 给定边界条件f²(0)=-0.3，f²(3)=3.3，求相应的三次样条插值函数s(x)并画出图形 解：问题1是第一类边界条

45、件的问题，执行三次样条插值程序后，在输入的五个窗口中按提示分别输入：1、0.2 、1、0, 1, 2, 3、0, 0.5, 2, 1.5，每次输入后用鼠标点击窗口的“OK”按扭，得如下输出结果：Mk的值为 -0.36, 2.52, -3.72, 0.36所求样条函数S(x)为S(t)=0.06 (1 - t) - 0.06 (1 - t)3 + 0.08 t + 0.42 t3 0£t<1S(t)=0.08 (2 - t) + 0.42 (2 - t)3 + 2.62 (-1 + t) - 0.62 (-1 + t)3 1£t<2S(t)=2.62 (3 - t) - 0.62 (3 - t)3 +