《数值分析》教案6

数值分析教案6 1.7.2三次样条插值的基本原理三次样条插值也是一种分段插值方法，用分段的三次多项式构造成一个整体上具有函数、一阶和二阶导函数连续的函数，近似地替代已知函数)(x f，“样条”一词源于过去绘图员使用的一种绘图工具样条，它是用于富于弹性、能弯曲的木条（或塑料）制成的软尺，把它弯折靠近所有的基点用画笔沿着样条就可以画出连续基点的光滑曲线。 假设已知函数)(x f在区间,b a上的)1(?n个节点b x x x x x an n?1210?及其对应的函数值i iy x f?)(,),2,1,0(n i?，即给出)1(?n组样本点数据),(,),(),(1100n ny x yx yx?，可以构造一个定义在,b a上的函数)(x S，满足下述条件。 i iy x S?)(,),2,1,0(n i?)(x S在每个小区间,1?i ix x)1,2,1,0(?n i?上，都是一个三次多项式332210)(x a x a x a a x Si i i i i?（1-42）)(),(),(x S x S x S?在,b a上连续。 可见，)(x S是一个光滑的分段函数，这样的函数称为三次样条（Spline）插值函数。 构造的函数)(x S是由n个小区间上的分段函数组成，根据条件，每个小区间上构造出一个三次多项式，第i个小区间上的三次多项式为332210)(x ax axaax Si i i i i?，共有n个多项式，每个多项式有4个待定系数。 要确定这n个多项式，就需要确定4n个系数a3210,i i i iaaa)1,2,1,0(?n i?。 为此，应该找到包含这些系数的4n个独立方程。 根据)(x S满足的条件，在所有的节点上可得出)1(?n个条件方程i iy x S?)(,),2,1,0(n i?(1-43)根据)(x S满足的条件，除两端点外在所有节点上，又可得出)1(3?n个方程?)()()()()()(111i i i ii i i ii i i ix S x Sx S x Sx Sx S)1,2,1(?n i?(1-44)由式(1-43)和(1-44)式可知共有)24(?n个独立方程，还差两个。 通常的办法是在区间,b a的两个端点上各加一个条件，即称之为边界条件。 常用的边界条件有以下三种 （1）给定两端点处的导数值n nyx Syx S?)(,)(00，特别地，当00?ny y时，样条曲线在端点处呈水平状态。 （2）给定两端处的二阶导数值n nyx Syx S?)(,)(00，特别地，当0)()(0?nx Sx S时，称为自然边界条件。 （3）如果)(x f是以a b?为周期的周期函数，则)(x S也应是具有同样的周期的周期函数，在端点处需要满足)0()0(),0()0(?b Sa Sb Sa S这样，在已有的)24(?n个条件方程基础上，再加上任何一种边界条件，即可求出这4n个系数，从而就可以求得三次样条插值函数)(x S。 1.7.3以节点处的二阶导数值为参数的三次样条插值函数设i iMx S?)(),1,0(n i?，因为在小区间,1?i ix x上)()(x Sx Si?是三次多项式，故)(x Si?为线性函数，由Lagrange插值公式得iiiiiii iiii iii ihx xMhx xMx xx xMx xx xM x S111111)(?，,1?i ix x x（1-45）式中，i i ix x h?1对式（1-45）两边进行一次积分，得出iiiiiii iAhx xMhx xMx S?2) (2)()(2121，,1?i ix x x（1-46）再对式（1-46）两边进行一次积分，得出i i iiiiiii iBx xAhx xMhx xMx S?) (6) (6)()(3131，,1?i ix x x（1-47）式中iA，iB都是积分常数。 代入插值条件i i iy x S?)(,11)(?i i iyx S，得i iii i iyBhM x S?6)(2，12116)(?i i i iii i iyB hAhM xS i iyx,是已知的，解这两个方程，得出iA与iB，代入式（1-47）得出S)1,1,0 (666) (6)()(122113131?n ihx xh Myhx xhMyhx xMhx xMxii i iiii i iiiiiiii i?（1-48）求导得6) (2) (2)()(112121iiiii iiiiiii iM Mhhy yhxxMhxxMxS?，,1?i ixxx（1-49）)(x Si?在,b a上连续，所以在相邻两个小区间的分界点ix（节点）上取值相等)0()0(1?iii ixSxS，1,2,1?n i?(1-50）由式（1-49）和式（1-50）便可得出1111111163)0 (63)0(?iiiiiiiiiiiii iiiiiMhhy yhM xSMhhy yhM xS点,)0(1iiixxx?，点,)0(1?iiixxx，因此，在代入式（1-49）时区间长度分别用11?iiixxh和iiixxh?1。 移项可得11111111162?iiii iii iii iiiii iih h hy yhy yMhhhM Mhhh若令111111161?iiiiiii iiiiiiii iiihh hy yhy ychhhh hh?(1-51)它们均为常数，于是上式变成iiiii icM M M?112?，)1,2,1(?n i?（1-52）式（1-52）给出了含有)1(?n个参数nM M M,10?的)1(?n个方程，称为三次样条的M关系式，或按其动力学意义，称为三弯矩方程。 这)1(?n个方程中共有)1(?n个量nM MM,10?，要求出它们尚缺两个条件方程，这就需要借助边界条件。 将这三种边界条件中的任一种与三弯矩方程（1-51）联立，即可求出参数),1,0(n iMi?。 如果问题要求)(xS满足边界条件 （1），由式（1-49）得) (62)(010001000M Mhhy yMhy aS?) (62)(11111?n nnnn nnnnM Mhhy yMhy bS化简2?00010106yhy yhMM?111162nn nnnn nhy yyhMM（1-53）式(1-52)和式(1-53)联立，即得关于参数),1,0(n iMi?的)1(?n阶线性方程组，其矩阵形式为?nnnnnn ncMMMMMa12101210112211022222?(1-54)其中?11100010006,6,1,1nnnnnnnhy yyhcyhy yhc?。 这就是三对角方程组，可用追赶法求解。 如果问题要求)(xS满足边界条件 （2），即n nyM yM?,00此时方程组(1-52)中实际只有)1(?n个数，其矩阵形式的?nnnnnnnn nycy cMMMM112xx12211222212222?(1-55)这仍是三对角方程组。 如果问题要求)(xS满足周期边界条件 （3），将边界条件)0()0(0?nx SxS，nM M?0与方程组(1-52)联立，得到n阶方程组?nnnnn nnnMMMM a1211212222112222?(1-56)其中?11001101001016,nn nnnnnnnnhyyhy yhhch hhh hha?。 其求解过程类似于求解三对角方程组的追赶法。 容易验证，这三种情形所得到的方程组的系数矩阵均为严格对角占优阵，因此方程组有唯一解，将解代入(1-48)式，即得三次样条插值函数。 2.以节点处的导数值为参数的三次样条插值函数。 类似地，也可以节点处的导数值为参数来导出三次样条插值函数的表达式。 设i imxS?)(),1,0(n i?。 按分段三次Hermite插值公式，)(xS可表成,)() (2121)(11212112121?iiiiii iiiiiiiiiiiiiixxx mhx xx xmhx xx xyhxxhx xyhxxhx xxS(1-57)由二阶导数在节点处连续S)0()0(1?iiiixSx)1,1(?n i?可导出关于参数im的方程组iiiii idm m m?112?，)1,2,1(?n i?(1-58)其中?)()(311111111iiiiiiiiii iiiii iiiyyhay yhdhhhhhh?(1-59)式(1-58)称为三次样条的m关系式，或按其动力学意义，称为三转角方程。 方程组(1-57)含有)1(?n个量，)1(?n个方程，与对M关系式的讨论类似，增加边界条件 （1）， （2）后，可得关于参数im的三对角方程组；增加周期边界条件 （3）后，得n阶方程组。 这些方程组的系数矩阵均非奇异，故方程组有唯一解。 解出参数im)1,1,0(?n i?后代入式（1-56），即得三次样条插值函数。 【例1-20】已知函数表i0123ix0123iy02316求满足边界条件1)0(?y，0)3(?y的三次样条插值函数。 解如果用三弯矩方程求解，由已知，按式(1-51)和式(1-53)计算方程组的系数及右端项，结果如下ii?i?ic01611/21/2-321/21/23631-78将上述数据代入式(1-52)和式(1-53)得方程组?783636212/122/102/122/100123210MMMM其解为3170,3106,338,3283210?MMMM。 将此解代入式(1-48)，得1,0)31411 (3)1 (914937)1 (914919)1 (63280)0 (633826)1 (3286)0 (338)(233330?xx xxx xxxx xxxxS1,0)351081924 (31)(231?xxxxxSS1,0)52573232946 (31)(232?xxxxx如果用三转角方程求解，由式(1-59)得212121?29) (21)(2131xx?yyyyd21) (21)(21323122?yyyyd将上述数据及0,130?m m代入式(1-58)，得三转角方程21221,29212212121?m m mm其解为332,3221?mm。 将此解代入式(1-57)，即得三次样条插值函数的分段表达式?3,2)52573232946(312,1)351081924(311,0)31411 (3)(23232xxxxxxxxxxxxx S比较上面两种解法，对第一种边界条件，用三转角方程计算较简单。 小结计算三次样条插值函数的步骤为 （1）根据给定的点),(i iyx及相应的边界条件计算iiic a,?或id。 一般地，对第一种边界条件用三转角方程，对第二种边界条件用三弯矩方程较为简单； （2）解方程组(1-51)或(1-58)，求出参数?iM或?im； （3）将求得的参数代入式(1-48)或式(1-57)，即得三次样条插值函数)(xS的分段表示式。 三次样条插值函数有较好的收敛性。 定理设)(x f在,b a上四次连续可微，)(xS为)(x f在,b a上的三次样条插值函数，则)3,2,1,0(,)(max)()()4 (4