三次样条插值算法详解PPT学习教案

1、会计学1三次样条插值算法详解三次样条插值算法详解2样条函数的定义样条函数的定义定义定义4.1设区间设区间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次样条函数。次样条函数。第1页/共59页3处的函数值为在节点如果函数nxxxxf,)(10n

2、jyxfjj, 1 , 0,)(njyxSjj, 1 , 0,)(三次样条插值函数的定义三次样条插值函数的定义并且关于这个节点集的三次样条函数并且关于这个节点集的三次样条函数s(x)满足满足插值条插值条件：件：则称这个三次样条函数则称这个三次样条函数s(x)为为三次样条插值函数三次样条插值函数。第2页/共59页4则其必满足的三次样条插值函数是如果,)()(xfxSnjyxSjj, 1 , 0,)(1, 1,)()(limnjmxSxSjjxxj1, 1),()(lim njxSxSjxxj1, 1,)()(limnjyxSxSjjxxj三次样条插值函数的边界条件三次样条插值函数的边界条件插值条

3、件：插值条件：连续性条件：连续性条件：一阶导数连续条件：一阶导数连续条件：二阶导数连续条件：二阶导数连续条件：第3页/共59页5（1）因为）因为s(x)在每个小区间上是一个次小于三次的多在每个小区间上是一个次小于三次的多项式，故有四个未知系数；项式，故有四个未知系数；（2）因为）因为s(x)有有n分段，从而共有分段，从而共有4n个未知系数！个未知系数！（3）但）但插值条件与样条条件插值条件与样条条件仅给出仅给出4n-2个条件，无法个条件，无法定出定出4n个未知系数，还差个未知系数，还差2个条件！这个条件！这2个条件我们用个条件我们用边界条件边界条件给出！给出！ 第4页/共59页6通常我们对插值

4、多项式在两端点的状态加以要求也就是所谓的边界条件:第一边界条件：由区间端点处的一阶导数给出即30003()(),()(),nnns xmfxs xmfx第5页/共59页7第二边界条件：由区间端点处的二阶导数给出即30003()(),()(),nnns xMfxs xMfx 特殊情况为自然边界条件：由区间端点处的二阶导数恒为0给出即3003()0()0nns xMs xM 第6页/共59页8这样三次样条插值问题就分成三类！其实不止这三类!第三类又称周期边界条件：由区间端点处的函数值或导数值满足周期条件给出 303303303(0)(0)(0)(0)(0)(0)nnns xs xs xs xs x

5、s x 第7页/共59页9样条函数的例子样条函数的例子容易验证：容易验证：3215)18218(2115)1427163(1015)152611()(232323xxxxxxxxxxxxxS是满足如下数据的第一类边界样条插值问题解：是满足如下数据的第一类边界样条插值问题解：x0123y0000y10第8页/共59页10样条函数的例子样条函数的例子3215)18218(2115)1427163(1015)152611()(232323xxxxxxxxxxxxxS第9页/共59页11通常有通常有三转角法、三弯矩法、三转角法、三弯矩法、B样条基函数法样条基函数法。三次样条插值函数的求法这三种方法的基

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

7、虑到带一阶导数的分段三次Hermite插值多项式)()()()()(11111010)(3iiiiiiiihxxiihxxiihxxihxxiihmhmyyxH1, 1 , 0,1nixxhiiinnnxxxxHxxxxHxxxxHxH1)1(321)1(310)0(33)()()()(第11页/共59页13njmxSjj, 1 , 0,)(我们采用待定一阶导数的方法即设我们采用待定一阶导数的方法即设因为分段三次Hermite插值多项式已经至少是一阶连续可导了，为了让它成为三次样条函数只需确定节点处的一阶导数使这些节点处的二阶导数连续即可！1, 1),0()0( nixSxSii)()()()

8、()(11111010iiiiiiiihxxiihxxiihxxihxxihmhmyyxS1, 1 , 0,11nixxhxxxiiiii212) 1()(,) 1)(12()(xxxxxx第12页/共59页14并整理后得求二阶导数对,)(xSi)()2(6)(131iiiiiiyyhxxxxS 11122624642iiiiiiiixxxxxxmmhh112642lim( )(0)()iiiiiixxiiiSxSxyymmhhh112111624lim( )(0)()iiiiiixxiiiSxSxyymmhhh第13页/共59页15(0)(0),1,2,1iiSxSxin由于在内部节点处二阶

9、导数连续条件：11111)11(21iiiiiiimhmhhmh)(321121iiiiiihyyhyy1, 1ni111112iiiiiiiiimhhhmmhhh)(3111111iiiiiiiiiiiihyyhhhhyyhhh整理化简后得：整理化简后得：第14页/共59页16112,1,1iiiiiimmmg in1iiiihhh111iiiiihhh 1113()iiiiiiiiiyyyyghh个未知量个方程共个1,1nn称为三转角法基本方程组以上推导还没有考虑边界条件！针对不同类型的三次样条问题，就可以导出不同的方程组！第15页/共59页17第一类三次样条插值问题方程组第一类三次样条插

10、值问题方程组00)(mxSnnmxS)(基本方程组化为基本方程组化为n-1阶方程组阶方程组0112112fgmmkkkkkkgmmm1122, 3 ,2nknnnnnnfgmm111212由于已知：由于已知：化为矩阵形式化为矩阵形式第16页/共59页181111022223333422221111222222nnnnnnnnnmgmmgmgmgmgm这是一个严格对角占优的三对角方程组，这是一个严格对角占优的三对角方程组，用追赶法可以求解！用追赶法可以求解！第17页/共59页19第二类三次样条插值问题的方程组第二类三次样条插值问题的方程组00)(MxS nnMxS )(时，称为自然边界条件00n

11、MM由于已知：由于已知：)(6)(01200yyhxS 004mh102mh0M)(6)(121 nnnnyyhxS112nnmhnnmh14nM故得：故得：第18页/共59页200000110232MhhyymmnnnnnnnMhhyymm23211110gng稍加整理得212222121132211nnnnmmmmm1210nnggggg1210联合基本方程组得一个联合基本方程组得一个n+1阶三对角方程组，阶三对角方程组，化成矩阵形式为：仍然是严格对角占优化成矩阵形式为：仍然是严格对角占优第19页/共59页21第三类样条插值问题的方程组立即可得下式：nnmmxSxS00)0()0(由于：)

12、0()0(0 nxSxS)2(2)(6)0(01001200mmhyyhxS )2(2)(6)0(11121nnnnnnnmmhyyhxS 第20页/共59页22nnnnngmmm211110013nnnnnnhyyhyygnnnnnnhhhhhh1,101100其中：联合基本方程得一个广义三对角或周期三对角方程组：联合基本方程得一个广义三对角或周期三对角方程组：1111112222311222222nnnnnnmgmgmg这个方程组的系数矩阵仍然是严格对角占优阵！第21页/共59页23求解这些不同类型的样条插值问题的方程组，我们可得求解这些不同类型的样条插值问题的方程组，我们可得所要待定的一

13、阶导数：所要待定的一阶导数：称为三次样条插值问题三转角公式！nmmmm,210再代入再代入S(x)的每一段表达式，就求得三次样条函数的表达式的每一段表达式，就求得三次样条函数的表达式!)()()()()(11111010iiiiiiiihxxiihxxiihxxihxxihmhmyyxS1,.,2 , 1 , 0ni第22页/共59页24例1. 对于给定的节点及函数值1 2 4 5() 1 3 4 2kkxf x的近似值并求插值函数的三次样条求满足自然边界条件)3(),(0)()(0fxSxSxSn 解:这是自然边界条件下的样条问题。3213123113221kkkkhhh11kkkkhhhk

14、 1第23页/共59页25291g272g60g63g我们可以将上述计算列于表中：)(3111kkkkkkkkkhyyhyyg00001023MhhyygnnnnnnMhhyyg23111k0123xk1245yk1342mk?Mk0?0hk121*k*2/31/3*k*1/32/3*gk69/2-7/2-6第24页/共59页26213/223/13/123/2123210mmmm3210gggg由些得如下方程组：819,45,47,8173210mmmm利用三转角公式：)()()()()(11111010iiiiiiiihxxiihxxiihxxihxxihmhmyyxS第25页/共59页2

15、7)()()()()(0100010101100010000hxxhxxhxxhxxhmhmyyxS)2() 1()2(3) 1(147181700 xxxx21147838123xxxx)()()()()(1211121111211102011hxxhxxhxxhxxhmhmyyxS)()()(4)(32412522127240220 xxxx42147838123xxxx同样可以求得第三段表达式！第26页/共59页283232321371 12884137( )1 2488434510333 45884xxxxS xxxxxxxxx 417)3()3( Sf第27页/共59页29第28页/

16、共59页30三弯矩法三弯矩法:待定二阶导数待定二阶导数() (0,1,2, )kkMSxkn选择二阶导数作为待定参数选择二阶导数作为待定参数:(1)1111()( )() , ,()kkkkkkkkkkkkkkSxMMMSxxxMxxxSxMh由于三次样条由于三次样条S(x)是三次多项式是三次多项式,故它的二阶导数是一次多项故它的二阶导数是一次多项式式,从而从而1(0,1,1)kkkhxx kn思考思考:(1)的原因的原因?第29页/共59页31(2)21( )()()2kkkkkkkkMMSxxxMxxph(3)321( )()()()62kkkkkkkkkkMMMSxxxxxpxxqh()

17、kkkkkSxyqy111111112()(),6()=)kkkkkkkkkkkkkkMMSxSxypf xxhSxy利用第30页/共59页32321234( )()()() (2.7)kkkkkkkS xsxxsxxsxxs1121314622,6kkkkkkkkkkkkkkkkMMshMsMMspf xxhsqy从而推导出了三次样条从而推导出了三次样条S(x)在第在第k个小区间个小区间xk,xk+1上的上的表达式为表达式为:0,1,2,1kn它的系数都是它的系数都是用二阶导数与用二阶导数与函数值表示函数值表示!第31页/共59页3311()() (1,2,1)()k-1k()kk+1kkk

18、kkkkkSxS xknSxS x（是在区间， ；是在区间 ，）1126 (2.16)kkkkkkMMMd对所有中间节点对所有中间节点xk,k=1,2,n-1,左边小区间与右边小区间上左边小区间与右边小区间上的三次多项式的一阶导数应当连续的三次多项式的一阶导数应当连续!确定二阶导数确定二阶导数kx1kx1kx1(0,1,1)kkkhxx kn第32页/共59页34三弯矩法基本方程三弯矩法基本方程1126 (2.16)1,2, -1kkkkkkMMMdkn11(1) ,1,kkkkkkhhh 110011, (1,2,1)(2.17),(2.18), , ,kkkknnndf xxxkndf x

19、 xdf xx注意到这个基本方程只注意到这个基本方程只包括了包括了n-1个方程个方程!但却但却有有n个二阶导数需要待定个二阶导数需要待定,这是一个欠定方程组这是一个欠定方程组,还还需要根据边界条件再确定两个方程需要根据边界条件再确定两个方程!第33页/共59页35曲率调整样条曲率调整样条00(),()nnMfxMfx这种样条的边界条件是已知两端点的二阶导数值这种样条的边界条件是已知两端点的二阶导数值!这样从三弯矩基本方程可以导数确定其它这样从三弯矩基本方程可以导数确定其它n-2个待定参数个待定参数的方程组的方程组:1111022222222111126262626nnnnnnnnnMdMMdM

20、dMdM第34页/共59页36自然样条自然样条00,0nMM这种样条的边界条件是这种样条的边界条件是:已知两端点的二阶导数值为已知两端点的二阶导数值为0!这样从三弯矩基本方程可以导数确定其它这样从三弯矩基本方程可以导数确定其它n-2个待定参数个待定参数的方程组的方程组:1112222222211126262626nnnnnnnMdMdMdMd第35页/共59页37固支样条固支样条00(),()nnmfxmfx这种样条的边界条件是这种样条的边界条件是:已知两端点的一阶导数值已知两端点的一阶导数值!根据前面推导过程中得到的样条函数根据前面推导过程中得到的样条函数S(x)的一阶导数的表的一阶导数的表

21、达式达式(2.11),得方程得方程00001000111()26()/()26()/nnnnnnnnSxmMMdmhSxmMMmdh第36页/共59页38固支样条固支样条这样从三弯矩基本方程可以导数确定这样从三弯矩基本方程可以导数确定n个待定个待定参数的方程组参数的方程组:00001111111116()/2162 (2.20)626()/12nnnnnnnnMdmhMdMdMmdh 第37页/共59页39非扭结样条非扭结样条0121( )( ),( )( )nnSxS x SxSx这种样条的边界条件是这种样条的边界条件是:要求样条要求样条S(x)在开始的两个小在开始的两个小区间区间x0,x1

22、,x1,x2上的三阶导数相同上的三阶导数相同,在最后两个小区间在最后两个小区间xn-2,xn-1,xn-1,xn上的三阶导数相同上的三阶导数相同.对表达式对表达式(2.9)再求一次导数得方程再求一次导数得方程1( )() (2.9)kkkkkkMMSxxxMh0112111211() (2.21)1()nnnnnMMMMMM第38页/共59页40非扭结样条非扭结样条再由三弯矩基本方程再由三弯矩基本方程,可得可得11222222221162162 (2.22)62612nnnnnnMdaaMdMdMdbb第39页/共59页41周期样条周期样条01126nnnnnnMMMMMd这种样条的边界条件是

23、这种样条的边界条件是:要求样条要求样条S(x)及其导数是以区及其导数是以区间长度间长度xn-x0为周期的函数即为周期的函数即000()(),()(),()()nnnS xS xS xS xSxSx这些条件可以确定如下两个方程这些条件可以确定如下两个方程:第40页/共59页42再由三弯矩基本方程再由三弯矩基本方程,可得可得1111222211116262 (2.23)6226nnnnnnnndMdMdMMd周期样条周期样条第41页/共59页431.1.在在MatlabMatlab中数据点称之为断点。如果三次样条插值没有边中数据点称之为断点。如果三次样条插值没有边界条件，最常用的方法，就是采用非扭

24、结（界条件，最常用的方法，就是采用非扭结（not-a-knotnot-a-knot）条件。这个条件强迫条件。这个条件强迫第第1 1个和第个和第2 2个三次多项式的三阶导数个三次多项式的三阶导数相等，对最后一个和倒数第相等，对最后一个和倒数第2 2个三次多项式也做同样地处个三次多项式也做同样地处理。理。2.2.MatlabMatlab中三次样条插值也有现成的函数：中三次样条插值也有现成的函数：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,co

25、nds),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是插值点的是插值点的函数值。函数值。3.3.对于三次样条插值，我们提倡使用函数对于三次样条插值，我们提倡使用函数csapecsape，csapecsape的返的返回值是回值是pppp形式，要求插值点的近似函数值，必须调用函数形式，要求插值点的近似函数值，必须调用函数ppvalpp

26、val。MATLAB中三次样条函数法实现中三次样条函数法实现第42页/共59页44pp=csape(x0,y0,conds,valconds)pp=csape(x0,y0,conds,valconds)condsconds指定插值的边界条件，其值可为：指定插值的边界条件，其值可为：complete complete 边界为一阶导数，一阶导数的值在边界为一阶导数，一阶导数的值在valcondsvalconds参数中给出。参数中给出。not-a-knot not-a-knot 非扭结条件非扭结条件periodic periodic 周期条件周期条件second second 边界为二阶导数，二阶导

27、数的值在边界为二阶导数，二阶导数的值在valcondsvalconds参数中给出，若忽略参数中给出，若忽略valcondsvalconds参数，二阶导数的缺省值为参数，二阶导数的缺省值为0, 00, 0。MATLAB中三次样条函数法实现中三次样条函数法实现第43页/共59页45例例2：第一边界条件的例题：第一边界条件的例题x1245y1342y17/8-19/8clear;x=1,2,4,5;y=1,3,4,2;pp=csape(x,y,complete,17/8,-19/8);pp.coefsMATLAB代码代码第44页/共59页46第一边界条件的例题第一边界条件的例题pp = form:

28、pp breaks: 1 2 4 5 coefs: 3x4 double pieces: 3 order: 4 dim: 1pp.coefs -0.1250 0 2.1250 1.0000 -0.1250 -0.3750 1.7500 3.0000 0.3750 -1.1250 -1.2500 4.00003323320.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 分段多项式结构与系数矩阵分段多项式结构与系数矩阵第45页/共59页47如下代码求解上述样条问题：如下代码求解上述样条问

29、题：x-4-3-2-101234y00.151.122.362.361.460.490.060y00例例3：自然边界条件例题：自然边界条件例题第46页/共59页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程序程序第47页/共59页49pp = form: pp breaks: -4 -3 -2 -1 0 1 2 3 4 coefs: 8x

30、4 double pieces: 8 order: 4 dim: 1分段多项式结构：分段多项式结构：第48页/共59页50 pp.coefs= 0.18085603829161 0.00000000000000 -0.03085603829161 0 -0.08428019145803 0.54256811487482 0.51171207658321 0.15000000000000 -0.39373527245950 0.28972754050074 1.34400773195876 1.12000000000000 0.14922128129602 -0.89147827687776 0

31、.74225699558174 2.36000000000000 0.13685014727540 -0.44381443298969 -0.59303571428571 2.36000000000000 0.13337812960236 -0.03326399116348 -1.07011413843888 1.46000000000000 -0.06036266568483 0.36687039764359 -0.73650773195876 0.49000000000000 -0.06192746686303 0.18578240058910 -0.18385493372607 0.06000000000000分段多项式的系数矩阵：分段多项式的系数矩阵：第49页/共59页51图像图像第50页/共59页5221 1,1125y