三次样条插值算法详解

《三次样条插值算法详解》由会员分享，可在线阅读，更多相关《三次样条插值算法详解（59页珍藏版）》请在装配图网上搜索。

1、1 三次样条插值三次样条插值鉴于高次插值不收敛又不稳定的特点，低次插值既具有收敛鉴于高次插值不收敛又不稳定的特点，低次插值既具有收敛性又具有稳定性，因此低次值更具有实用价值，但是低次插性又具有稳定性，因此低次值更具有实用价值，但是低次插值的光滑性较差，比如分段线性插值多项式在插值区间中仅值的光滑性较差，比如分段线性插值多项式在插值区间中仅具有连续性，在插值节点处有棱角，一阶导数不存在；分段具有连续性，在插值节点处有棱角，一阶导数不存在；分段三次三次Hermite插值多项式在插值区间中仅具有一阶导数即一插值多项式在插值区间中仅具有一阶导数即一阶光滑性但不具备二阶光滑性，不能满足某些实际应用比如阶

2、光滑性但不具备二阶光滑性，不能满足某些实际应用比如汽车、轮船、飞机等的外形中流线形设计。样条是在二十世汽车、轮船、飞机等的外形中流线形设计。样条是在二十世纪初期经常用于图样设计的一种富有弹性的细长条，多个样纪初期经常用于图样设计的一种富有弹性的细长条，多个样条互相弯曲连接后沿其边缘画出的曲线就是三次样条曲线。条互相弯曲连接后沿其边缘画出的曲线就是三次样条曲线。后来数学上对其进行了抽象，定义了后来数学上对其进行了抽象，定义了m次样条函数，并成为次样条函数，并成为数值逼近的重要研究分枝，进一步扩大了样条函数的应用范数值逼近的重要研究分枝，进一步扩大了样条函数的应用范围。围。2样条函数的定义样条函数

3、的定义定义定义设区间设区间a,b上给定一个节点划分上给定一个节点划分a=x0 x1xn-1xn=b如果存在正整数如果存在正整数k使得使得a,b上的分段函数上的分段函数s(x)满足如满足如下两条：下两条：(1)在在a,b上有直到上有直到k-1阶连续导数。阶连续导数。(2)在每个小区间在每个小区间xi,xi+1上是次数不大于上是次数不大于k的多项式。的多项式。则称则称分段函数分段函数s(x)是以是以(2.6)为节点集的为节点集的k次样条函数。次样条函数。3处的函数值为在节点如果函数nxxxxf,)(10njyxfjj,1,0,)(njyxSjj,1,0,)(三次样条插值函数的定义三次样条插值函数的

4、定义并且关于这个节点集的三次样条函数并且关于这个节点集的三次样条函数s(x)满足满足插值条插值条件：件：则称这个三次样条函数则称这个三次样条函数s(x)为为三次样条插值函数三次样条插值函数。4则其必满足的三次样条插值函数是如果,)()(xfxSnjyxSjj,1,0,)(1,1,)()(limnjmxSxSjjxxj1,1),()(lim njxSxSjxxj1,1,)()(limnjyxSxSjjxxj三次样条插值函数的边界条三次样条插值函数的边界条件件插值条件：插值条件：连续性条件：连续性条件：一阶导数连续条件：一阶导数连续条件：二阶导数连续条件：二阶导数连续条件：5（1）因为）因为s(x

5、)在每个小区间上是一个次小于三次的多在每个小区间上是一个次小于三次的多项式，故有四个未知系数；项式，故有四个未知系数；（2）因为）因为s(x)有有n分段，从而共有分段，从而共有4n个未知系数！个未知系数！（3）但）但插值条件与样条条件插值条件与样条条件仅给出仅给出4n-2个条件，无法个条件，无法定出定出4n个未知系数，还差个未知系数，还差2个条件！这个条件！这2个条件我们用个条件我们用边界条件边界条件给出！给出！6通常我们对插值多项式在两端点的状态加以要求也就是所谓的边界条件:第一边界条件：由区间端点处的一阶导数给出即第一边界条件：由区间端点处的一阶导数给出即30003()(),()(),nn

6、ns xmfxs xmfx7第二边界条件：由区间端点处的二阶导数给出即第二边界条件：由区间端点处的二阶导数给出即30003()(),()(),nnns xMfxs xMfx 特殊情况为特殊情况为自然边界条件自然边界条件：由区间端点处的二阶导数恒为由区间端点处的二阶导数恒为0给出即给出即3003()0()0nns xMs xM 8这样三次样条插这样三次样条插值问题就分成三值问题就分成三类！其实不止这类！其实不止这三类三类!第三类又称周期边界条件：第三类又称周期边界条件：由区间端点处的函数值或导数值满足周期条件给出由区间端点处的函数值或导数值满足周期条件给出 303303303(0)(0)(0)(

7、0)(0)(0)nnns xs xs xs xs xs x 9样条函数的例子样条函数的例子容易验证：容易验证：3215)18218(2115)1427163(1015)152611()(232323xxxxxxxxxxxxxS是满足如下数据的第一类边界样条插值问题解：是满足如下数据的第一类边界样条插值问题解：x0123y0000y1010样条函数的例子样条函数的例子3215)18218(2115)1427163(1015)152611()(232323xxxxxxxxxxxxxS11通常有通常有三转角法、三弯矩法、三转角法、三弯矩法、B样条基函数法样条基函数法。三次样条插值函数的求法三次样条插

8、值函数的求法这三种方法的基本思想是类似的，都是这三种方法的基本思想是类似的，都是通通过待定某些参数过待定某些参数来确定插值函数，但肯定来确定插值函数，但肯定不是待定不是待定4n个参数。而是利用已知条件将个参数。而是利用已知条件将待定参数减小到最少。待定参数减小到最少。比如：比如：待定一阶导数、待定二阶导数待定一阶导数、待定二阶导数、采采用基函数方法用基函数方法来确定插值函数。来确定插值函数。12三转角法三转角法:待定一阶数待定一阶数为了确定三次样条插值函数的表达式为了确定三次样条插值函数的表达式 S(x)，我们采用待定系数法来求解，我们待定什么系数呢？我们采用待定系数法来求解，我们待定什么系数

9、呢？考虑到带一阶导数的分段三次考虑到带一阶导数的分段三次Hermite插值多项式插值多项式)()()()()(11111010)(3iiiiiiiihxxiihxxiihxxihxxiihmhmyyxH1,1,0,1nixxhiiinnnxxxxHxxxxHxxxxHxH1)1(321)1(310)0(33)()()()(13njmxSjj,1,0,)(我们采用待定一阶导数的方法即设我们采用待定一阶导数的方法即设因为因为分段三次分段三次Hermite插值多项式已经至少是一阶连续插值多项式已经至少是一阶连续可导可导了，为了让它成为三次样条函数只需确定节点处了，为了让它成为三次样条函数只需确定节点

10、处的一阶导数使这些节点处的二阶导数连续即可！的一阶导数使这些节点处的二阶导数连续即可！1,1),0()0(nixSxSii)()()()()(11111010iiiiiiiihxxiihxxiihxxihxxihmhmyyxS1,1,0,11nixxhxxxiiiii212)1()(,)1)(12()(xxxxxx14并整理后得求二阶导数对,)(xSi)()2(6)(131iiiiiiyyhxxxxS 11122624642iiiiiiiixxxxxxmmhh112642lim()(0)()iiiiiixxiiiSxSxyymmhhh112111624lim()(0)()iiiiiixxiii

11、SxSxyymmhhh15(0)(0),1,2,1iiSxSxin由于在内部节点处二阶导数连续条件：11111)11(21iiiiiiimhmhhmh)(321121iiiiiihyyhyy1,1ni111112iiiiiiiiimhhhmmhhh)(3111111iiiiiiiiiiiihyyhhhhyyhhh整理化简后得：整理化简后得：16112,1,1iiiiiimmmg in1iiiihhh111iiiiihhh 1113()iiiiiiiiiyyyyghh个未知量个方程共个1,1nn称为三转称为三转角法基本角法基本方程组方程组以上推导还没有考虑边界以上推导还没有考虑边界条件！针对不同

12、类型的三条件！针对不同类型的三次样条问题，就可以导出次样条问题，就可以导出不同的方程组！不同的方程组！17第一类三次样条插值问题方程组第一类三次样条插值问题方程组00)(mxSnnmxS)(基本方程组化为基本方程组化为n-1阶方程组阶方程组0112112fgmmkkkkkkgmmm1122,3,2nknnnnnnfgmm111212由于已知：由于已知：化为矩阵形式化为矩阵形式181111022223333422221111222222nnnnnnnnnmgmmgmgmgmgm这是一个严格对角占优的三对角方程组，这是一个严格对角占优的三对角方程组，用追赶法可以求解！用追赶法可以求解！19第二类三

13、次样条插值问题的方程组第二类三次样条插值问题的方程组00)(MxS nnMxS)(时，称为自然边界条件00nMM由于已知：由于已知：)(6)(01200yyhxS 004mh102mh0M)(6)(121 nnnnyyhxS112nnmhnnmh14nM故得：故得：200000110232MhhyymmnnnnnnnMhhyymm23211110gng稍加整理得212222121132211nnnnmmmmm1210nnggggg1210联合基本方程组得一个联合基本方程组得一个n+1阶三对角方程组，阶三对角方程组，化成矩阵形式为：仍然是严格对角占优化成矩阵形式为：仍然是严格对角占优21第三类样

14、条插值问题的方程组立即可得下式：nnmmxSxS00)0()0(由于：)0()0(0 nxSxS)2(2)(6)0(01001200mmhyyhxS)2(2)(6)0(11121nnnnnnnmmhyyhxS 22nnnnngmmm211110013nnnnnnhyyhyygnnnnnnhhhhhh1,101100其中：联合基本方程得一个广义三对角或周期三对角方程组：联合基本方程得一个广义三对角或周期三对角方程组：1111112222311222222nnnnnnmgmgmg这个方程这个方程组的系数组的系数矩阵仍然矩阵仍然是严格对是严格对角占优阵！角占优阵！23求解这些不同类型的样条插值问题的

15、方程组，我们可得求解这些不同类型的样条插值问题的方程组，我们可得所要待定的一阶导数：所要待定的一阶导数：称为三次样条插称为三次样条插值问题三转角公值问题三转角公式！式！nmmmm,210再代入再代入S(x)的每一段表达式，就求得三次样条函数的表达式的每一段表达式，就求得三次样条函数的表达式!)()()()()(11111010iiiiiiiihxxiihxxiihxxihxxihmhmyyxS1,.,2,1,0ni24例1.对于给定的节点及函数值1 2 4 5()1 3 4 2kkxf x的近似值并求插值函数的三次样条求满足自然边界条件)3(),(0)()(0fxSxSxSn 解:这是自然边界

16、条件下的样条问题。3213123113221kkkkhhh11kkkkhhhk 125291g272g60g63g我们可以将上述计算列于表中：)(3111kkkkkkkkkhyyhyyg00001023MhhyygnnnnnnMhhyyg23111k0123xk1245yk1342mk?Mk0?0hk121*k*2/31/3*k*1/32/3*gk69/2-7/2-626213/223/13/123/2123210mmmm3210gggg由些得如下方程组：819,45,47,8173210mmmm利用三转角公式：)()()()()(11111010iiiiiiiihxxiihxxiihxxih

17、xxihmhmyyxS27)()()()()(0100010101100010000hxxhxxhxxhxxhmhmyyxS)2()1()2(3)1(147181700 xxxx21147838123xxxx)()()()()(1211121111211102011hxxhxxhxxhxxhmhmyyxS)()()(4)(32412522127240220 xxxx42147838123xxxx同样可以求得第三段表达式！283232321371 12884137()1 2488434510333 45884xxxxS xxxxxxxxx417)3()3(Sf2930三弯矩法三弯矩法:待定二阶导

18、数待定二阶导数()(0,1,2,)kkMSxkn选择二阶导数作为待定参数选择二阶导数作为待定参数:(1)1111()()(),()kkkkkkkkkkkkkkSxMMMSxxxMxxxSxMh由于三次样条由于三次样条S(x)是三次多项式是三次多项式,故它的二阶导数是一次多项故它的二阶导数是一次多项式式,从而从而1(0,1,1)kkkhxx kn思考思考:(1)的原因的原因?31(2)21()()()2kkkkkkkkMMSxxxMxxph(3)321()()()()62kkkkkkkkkkMMMSxxxxxpxxqh()kkkkkSxyqy111111112()(),6()=)kkkkkkkk

19、kkkkkkMMSxSxypf xxhSxy利用32321234()()()()(2.7)kkkkkkkS xsxxsxxsxxs1121314622,6kkkkkkkkkkkkkkkkMMshMsMMspf xxhsqy从而推导出了三次样条从而推导出了三次样条S(x)在第在第k个小区间个小区间xk,xk+1上的上的表达式为表达式为:0,1,2,1kn它的系数都是它的系数都是用二阶导数与用二阶导数与函数值表示函数值表示!3311()()(1,2,1)()k-1k()kk+1kkkkkkkkSxS xknSxS x（是在区间，；是在区间，）1126 (2.16)kkkkkkMMMd对所有中间节点

20、对所有中间节点xk,k=1,2,n-1,左边小区间与右边小区间上左边小区间与右边小区间上的三次多项式的一阶导数应当连续的三次多项式的一阶导数应当连续!确定二阶导数确定二阶导数kx1kx1kx1(0,1,1)kkkhxx kn34三弯矩法基本方程三弯矩法基本方程1126 (2.16)1,2,-1kkkkkkMMMdkn11(1),1,kkkkkkhhh 110011,(1,2,1)(2.17),(2.18),kkkknnndf xx xkndf x xdf xx注意到这个基本方程只注意到这个基本方程只包括了包括了n-1个方程个方程!但却但却有有n个二阶导数需要待定个二阶导数需要待定,这是一个欠定

21、方程组这是一个欠定方程组,还还需要根据边界条件再确定两个方程需要根据边界条件再确定两个方程!35曲率调整样条曲率调整样条00(),()nnMfxMfx这种样条的边界条件是已知两端点的二阶导数值这种样条的边界条件是已知两端点的二阶导数值!这样从三弯矩基本方程可以导数确定其它这样从三弯矩基本方程可以导数确定其它n-2个待定参数个待定参数的方程组的方程组:1111022222222111126262626nnnnnnnnnMdMMdMdMdM36自然样条自然样条00,0nMM这种样条的边界条件是这种样条的边界条件是:已知两端点的二阶导数值为已知两端点的二阶导数值为0!这样从三弯矩基本方程可以导数确定

22、其它这样从三弯矩基本方程可以导数确定其它n-2个待定参数个待定参数的方程组的方程组:1112222222211126262626nnnnnnnMdMdMdMd37固支样条固支样条00(),()nnmfxmfx这种样条的边界条件是这种样条的边界条件是:已知两端点的一阶导数值已知两端点的一阶导数值!根据前面推导过程中得到的样条函数根据前面推导过程中得到的样条函数S(x)的一阶导数的表的一阶导数的表达式达式(2.11),得方程得方程00001000111()26()/()26()/nnnnnnnnSxmMMdmhSxmMMmdh38固支样条固支样条这样从三弯矩基本方程可以导数确定这样从三弯矩基本方程

23、可以导数确定n个待定个待定参数的方程组参数的方程组:00001111111116()/2162 (2.20)626()/12nnnnnnnnMdmhMdMdMmdh 39非扭结样条非扭结样条0121()(),()()nnSxS x SxSx这种样条的边界条件是这种样条的边界条件是:要求样条要求样条S(x)在开始的两个小在开始的两个小区间区间x0,x1,x1,x2上的三阶导数相同上的三阶导数相同,在最后两个小区间在最后两个小区间xn-2,xn-1,xn-1,xn上的三阶导数相同上的三阶导数相同.对表达式对表达式(2.9)再求一次导数得方程再求一次导数得方程1()()(2.9)kkkkkkMMSx

24、xxMh0112111211()(2.21)1()nnnnnMMMMMM40非扭结样条非扭结样条再由三弯矩基本方程再由三弯矩基本方程,可得可得11222222221162162 (2.22)62612nnnnnnMdaaMdMdMdbb41周期样条周期样条01126nnnnnnMMMMMd这种样条的边界条件是这种样条的边界条件是:要求样条要求样条S(x)及其导数是以区及其导数是以区间长度间长度xn-x0为周期的函数即为周期的函数即000()(),()(),()()nnnS xS xS xS xSxSx这些条件可以确定如下两个方程这些条件可以确定如下两个方程:42再由三弯矩基本方程再由三弯矩基本

25、方程,可得可得1111222211116262 (2.23)6226nnnnnnnndMdMdMMd 周期样条周期样条43在在MatlabMatlab中数据点称之为断点。如果三次样条插值没有边中数据点称之为断点。如果三次样条插值没有边界条件，最常用的方法，就是采用非扭结（界条件，最常用的方法，就是采用非扭结（not-a-knotnot-a-knot）条件。这个条件强迫条件。这个条件强迫第第1 1个和第个和第2 2个三次多项式的三阶导数个三次多项式的三阶导数相等，对最后一个和倒数第相等，对最后一个和倒数第2 2个三次多项式也做同样地处理。个三次多项式也做同样地处理。MatlabMatlab中三次

26、样条插值也有现成的函数：中三次样条插值也有现成的函数：y=interp1(x0,y0,x,spline)y=interp1(x0,y0,x,spline)；y=spline(x0,y0,x)y=spline(x0,y0,x)；pp=csape(x0,y0,conds),pp=csape(x0,y0,conds),pp=csape(x0,y0,conds,valconds)pp=csape(x0,y0,conds,valconds)，y=ppval(pp,x)y=ppval(pp,x)。其中其中x0,y0 x0,y0是已知数据点，是已知数据点，x x是插值点，是插值点，y y是插值点的函是插值点

27、的函数值。数值。对于三次样条插值，我们提倡使用函数对于三次样条插值，我们提倡使用函数csapecsape，csapecsape的返的返回值是回值是pppp形式，要求插值点的近似函数值，必须调用函数形式，要求插值点的近似函数值，必须调用函数ppvalppval。MATLAB中三次样条函数法实现中三次样条函数法实现44pp=csape(x0,y0,conds,valconds)pp=csape(x0,y0,conds,valconds)condsconds指定插值的边界条件，其值可为：指定插值的边界条件，其值可为：complete complete 边界为一阶导数，一阶导数的值在边界为一阶导数，一

28、阶导数的值在valcondsvalconds参数中给出。参数中给出。not-a-knot not-a-knot 非扭结条件非扭结条件periodic periodic 周期条件周期条件second second 边界为二阶导数，二阶导数的值在边界为二阶导数，二阶导数的值在valcondsvalconds参数中给出，若忽略参数中给出，若忽略valcondsvalconds参数，二阶导数的参数，二阶导数的缺省值为缺省值为0,00,0。MATLAB中三次样条函数法实现中三次样条函数法实现45例例2：第一边界条件的例题：第一边界条件的例题x1245y1342y17/8-19/8clear;x=1,2,

29、4,5;y=1,3,4,2;pp=csape(x,y,complete,17/8,-19/8);MATLAB代码代码46第一边界条件的例题第一边界条件的例题pp=form:pp breaks:1 2 4 5 coefs:3x4 double pieces:3 order:4 dim:13323320.125(1)2.125(1)1()0.125(2)0.375(2)1,75(2)30.375(4)1.125(4)1.25(4)4xxs xxxxxxx 分段多项式结构与系数矩阵分段多项式结构与系数矩阵47如下代码求解上述样条问题：如下代码求解上述样条问题：x-4-3-2-101234y00.15

30、1.122.362.361.460.490.060y00例例3：自然边界条件例题：自然边界条件例题48x=-4-3-2-101234;y=00.151.122.362.361.460.490.060;pp=csape(x,y,second);xx=-4:0.01:4;yy=ppval(pp,xx);hold on;plot(x,y,ok);plot(xx,yy,k-);hold off;MATLAB程序程序49pp=form:pp breaks:-4-3-2-1 0 1 2 3 4 coefs:8x4 double pieces:8 order:4 dim:1分段多项式结构：分段多项式结构：5

31、0 pp.coefs=分段多项式的系数矩阵：分段多项式的系数矩阵：51图像图像5221 1,1125yxx 例例4：对如下：对如下Runge现象中的函数，求用现象中的函数，求用n分点作节分点作节点的三次样条插值多项式点的三次样条插值多项式s3(x)的图象。的图象。取取n=5,10,15,20等，将区间等分成等，将区间等分成n份，取分点作为插值份，取分点作为插值节点，利用节点，利用matlab函数函数csape()作出它的三次样条插值函作出它的三次样条插值函数，并作出这些插值函数与被插函数的图像。数，并作出这些插值函数与被插函数的图像。但由于插值节点数目较多，故不列出插值函数的表达式，但由于插值

32、节点数目较多，故不列出插值函数的表达式，只画出它们的图象与被插函数的图象的复合，从这些图象只画出它们的图象与被插函数的图象的复合，从这些图象可以看出，随着节点数目的增多，样条插值函数图象越来可以看出，随着节点数目的增多，样条插值函数图象越来越符合被插函数，因此避免越符合被插函数，因此避免Runge现象，也就是说三次样现象，也就是说三次样条插值函数具有收敛性。条插值函数具有收敛性。5354555657次样条插值函数，满足任意边界条件的三为节点是以设,),1,0()(,)(2nkxxSbaCxfk,min,max,10101iniiniiiihhhxxh设时则当 ch)()(,)()(xfxfba

33、xSxS和上一致收敛到在和三次样条函数的一致收敛性三次样条函数的一致收敛性,)()(lim),()(lim00baxxfxSxfxShh这个定理告诉我们：当插值节点较这个定理告诉我们：当插值节点较多时，三次样条插值函数可以十分多时，三次样条插值函数可以十分接近被插函数！不会产生接近被插函数！不会产生Runge现现象中的不收敛情况！象中的不收敛情况！58设被插函数设被插函数f(x)f(x)有直到四阶的连续导数，则对第一、二类样有直到四阶的连续导数，则对第一、二类样条插值问题的样条插值函数条插值问题的样条插值函数S(x)S(x)有如下估计：有如下估计：)(max3845)()(max)4(4xfhxSxfbxabxa三次样条函数的插值余项三次样条函数的插值余项)(max241)()(max)4(3xfhxSxfbxabxa)(max83)()(max)4(2xfhxSxfbxabxa inihh10max59在三次样条插值法中，插值节点处函数值的波动，只对这个在三次样条插值法中，插值节点处函数值的波动，只对这个节点两边的分段有影响，而对离该点较远的分段的影响会逐节点两边的分段有影响，而对离该点较远的分段的影响会逐渐变小，因此样条插值法具有较好的稳定性！渐变小，因此样条插值法具有较好的稳定性！三次样条函数法的稳定性三次样条函数法的稳定性