三次样条插值算法详解

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

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

3、编辑课件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次样条函数。次样条函数。编辑课件编辑课件3处的函数值为在节点如果函数nxxxxf,)(10njyxfjj, 1 , 0,)(njyx

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

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

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

7、界条件：第三类又称周期边界条件：由区间端点处的函数值或导数值满足周期条件给出由区间端点处的函数值或导数值满足周期条件给出 303303303(0)(0)(0)(0)(0)(0)nnns xs xs xs xs xs x 编辑课件编辑课件9样条函数的例子样条函数的例子容易验证：容易验证：3215)18218(2115)1427163(1015)152611()(232323xxxxxxxxxxxxxS是满足如下数据的第一类边界样条插值问题解：是满足如下数据的第一类边界样条插值问题解：x0123y0000y10编辑课件编辑课件10样条函数的例子样条函数的例子3215)18218(2115)1427

8、163(1015)152611()(232323xxxxxxxxxxxxxS编辑课件编辑课件11通常有通常有三转角法、三弯矩法、三转角法、三弯矩法、B样条基函数法样条基函数法。三次样条插值函数的求法三次样条插值函数的求法这三种方法的基本思想是类似的，都是这三种方法的基本思想是类似的，都是通通过待定某些参数过待定某些参数来确定插值函数，但肯定来确定插值函数，但肯定不是待定不是待定4n个参数。而是利用已知条件将个参数。而是利用已知条件将待定参数减小到最少。待定参数减小到最少。 比如：比如：待定一阶导数、待定二阶导数待定一阶导数、待定二阶导数、采采用基函数方法用基函数方法来确定插值函数。来确定插值函

9、数。编辑课件编辑课件12三转角法三转角法:待定一阶数待定一阶数为了确定三次样条插值函数的表达式为了确定三次样条插值函数的表达式 S(x)，我们采用待定系数法来求解，我们待定什么系数呢？我们采用待定系数法来求解，我们待定什么系数呢？考虑到带一阶导数的分段三次考虑到带一阶导数的分段三次Hermite插值多项式插值多项式)()()()()(11111010)(3iiiiiiiihxxiihxxiihxxihxxiihmhmyyxH1, 1 , 0,1nixxhiiinnnxxxxHxxxxHxxxxHxH1)1(321)1(310)0(33)()()()(编辑课件编辑课件13njmxSjj, 1 ,

10、 0,)(我们采用待定一阶导数的方法即设我们采用待定一阶导数的方法即设因为因为分段三次分段三次Hermite插值多项式已经至少是一阶连续插值多项式已经至少是一阶连续可导可导了，为了让它成为三次样条函数只需确定节点处了，为了让它成为三次样条函数只需确定节点处的一阶导数使这些节点处的二阶导数连续即可！的一阶导数使这些节点处的二阶导数连续即可！1, 1),0()0( nixSxSii)()()()()(11111010iiiiiiiihxxiihxxiihxxihxxihmhmyyxS1, 1 , 0,11nixxhxxxiiiii212) 1()(,) 1)(12()(xxxxxx编辑课件编辑课件

11、14并整理后得求二阶导数对,)(xSi)()2(6)(131iiiiiiyyhxxxxS 11122624642iiiiiiiixxxxxxmmhh112642lim( )(0)()iiiiiixxiiiSxSxyymmhhh112111624lim( )(0)()iiiiiixxiiiSxSxyymmhhh编辑课件编辑课件15(0)(0),1,2,1iiSxSxin由于在内部节点处二阶导数连续条件：11111)11(21iiiiiiimhmhhmh)(321121iiiiiihyyhyy1, 1ni111112iiiiiiiiimhhhmmhhh)(3111111iiiiiiiiiiiihy

12、yhhhhyyhhh整理化简后得：整理化简后得：编辑课件编辑课件16112,1,1iiiiiimmmg in1iiiihhh111iiiiihhh 1113()iiiiiiiiiyyyyghh个未知量个方程共个1,1nn称为三转称为三转角法基本角法基本方程组方程组以上推导还没有考虑边界以上推导还没有考虑边界条件！针对不同类型的三条件！针对不同类型的三次样条问题，就可以导出次样条问题，就可以导出不同的方程组！不同的方程组！编辑课件编辑课件17第一类三次样条插值问题方程组第一类三次样条插值问题方程组00)(mxSnnmxS)(基本方程组化为基本方程组化为n-1阶方程组阶方程组0112112fgmm

13、kkkkkkgmmm1122, 3 ,2nknnnnnnfgmm111212由于已知：由于已知：化为矩阵形式化为矩阵形式编辑课件编辑课件181111022223333422221111222222nnnnnnnnnmgmmgmgmgmgm这是一个严格对角占优的三对角方程组，这是一个严格对角占优的三对角方程组，用追赶法可以求解！用追赶法可以求解！编辑课件编辑课件19第二类三次样条插值问题的方程组第二类三次样条插值问题的方程组00)(MxS nnMxS )(时，称为自然边界条件00nMM由于已知：由于已知：)(6)(01200yyhxS 004mh102mh0M)(6)(121 nnnnyyhxS

14、112nnmhnnmh14nM故得：故得：编辑课件编辑课件200000110232MhhyymmnnnnnnnMhhyymm23211110gng稍加整理得212222121132211nnnnmmmmm1210nnggggg1210联合基本方程组得一个联合基本方程组得一个n+1阶三对角方程组，阶三对角方程组，化成矩阵形式为：仍然是严格对角占优化成矩阵形式为：仍然是严格对角占优编辑课件编辑课件21第三类样条插值问题的方程组立即可得下式：nnmmxSxS00)0()0(由于：)0()0(0 nxSxS)2(2)(6)0(01001200mmhyyhxS )2(2)(6)0(11121nnnnnn

15、nmmhyyhxS 编辑课件编辑课件22nnnnngmmm211110013nnnnnnhyyhyygnnnnnnhhhhhh1,101100其中：联合基本方程得一个广义三对角或周期三对角方程组：联合基本方程得一个广义三对角或周期三对角方程组：1111112222311222222nnnnnnmgmgmg这个方程这个方程组的系数组的系数矩阵仍然矩阵仍然是严格对是严格对角占优阵！角占优阵！编辑课件编辑课件23求解这些不同类型的样条插值问题的方程组，我们可得求解这些不同类型的样条插值问题的方程组，我们可得所要待定的一阶导数：所要待定的一阶导数：称为三次样条插称为三次样条插值问题三转角公值问题三转角

16、公式！式！nmmmm,210再代入再代入S(x)的每一段表达式，就求得三次样条函数的表达式的每一段表达式，就求得三次样条函数的表达式!)()()()()(11111010iiiiiiiihxxiihxxiihxxihxxihmhmyyxS1,.,2 , 1 , 0ni编辑课件编辑课件24例1. 对于给定的节点及函数值1 2 4 5() 1 3 4 2kkxf x的近似值并求插值函数的三次样条求满足自然边界条件)3(),(0)()(0fxSxSxSn 解:这是自然边界条件下的样条问题。3213123113221kkkkhhh11kkkkhhhk 1编辑课件编辑课件25291g272g60g63g

17、我们可以将上述计算列于表中：)(3111kkkkkkkkkhyyhyyg00001023MhhyygnnnnnnMhhyyg23111k0123xk1245yk1342mk?Mk0?0hk121*k*2/31/3*k*1/32/3*gk69/2-7/2-6编辑课件编辑课件26213/223/13/123/2123210mmmm3210gggg由些得如下方程组：819,45,47,8173210mmmm利用三转角公式：)()()()()(11111010iiiiiiiihxxiihxxiihxxihxxihmhmyyxS编辑课件编辑课件27)()()()()(010001010110001000

18、0hxxhxxhxxhxxhmhmyyxS)2() 1()2(3) 1(147181700 xxxx21147838123xxxx)()()()()(1211121111211102011hxxhxxhxxhxxhmhmyyxS)()()(4)(32412522127240220 xxxx42147838123xxxx同样可以求得第三段表达式！编辑课件编辑课件283232321371 12884137( )1 2488434510333 45884xxxxS xxxxxxxxx417)3()3( Sf编辑课件编辑课件29编辑课件编辑课件30三弯矩法三弯矩法:待定二阶导数待定二阶导数() (0,

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

20、kkkkkkkkkMMSxSxypf 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

21、 (2.16)kkkkkkMMMd对所有中间节点对所有中间节点xk,k=1,2,n-1,左边小区间与右边小区间上左边小区间与右边小区间上的三次多项式的一阶导数应当连续的三次多项式的一阶导数应当连续!确定二阶导数确定二阶导数kx1kx1kx1(0,1,1)kkkhxx kn编辑课件编辑课件34三弯矩法基本方程三弯矩法基本方程1126 (2.16)1,2, -1kkkkkkMMMdkn11(1) ,1,kkkkkkhhh 110011, (1,2,1)(2.17),(2.18), , ,kkkknnndf xxxkndf x xdf xx注意到这个基本方程只注意到这个基本方程只包括了包括了n-1个

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

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

24、)/nnnnnnnnSxmMMdmhSxmMMmdh编辑课件编辑课件38固支样条固支样条这样从三弯矩基本方程可以导数确定这样从三弯矩基本方程可以导数确定n个待定个待定参数的方程组参数的方程组:00001111111116()/2162 (2.20)626()/12nnnnnnnnMdmhMdMdMmdh编辑课件编辑课件39非扭结样条非扭结样条0121( )( ),( )( )nnSxS x SxSx这种样条的边界条件是这种样条的边界条件是:要求样条要求样条S(x)在开始的两个小在开始的两个小区间区间x0,x1,x1,x2上的三阶导数相同上的三阶导数相同,在最后两个小区间在最后两个小区间xn-2

25、,xn-1,xn-1,xn上的三阶导数相同上的三阶导数相同.对表达式对表达式(2.9)再求一次导数得方程再求一次导数得方程1( )() (2.9)kkkkkkMMSxxxMh0112111211() (2.21)1()nnnnnMMMMMM编辑课件编辑课件40非扭结样条非扭结样条再由三弯矩基本方程再由三弯矩基本方程,可得可得11222222221162162 (2.22)62612nnnnnnMdaaMdMdMdbb编辑课件编辑课件41周期样条周期样条01126nnnnnnMMMMMd这种样条的边界条件是这种样条的边界条件是:要求样条要求样条S(x)及其导数是以区及其导数是以区间长度间长度xn

26、-x0为周期的函数即为周期的函数即000()(),()(),()()nnnS xS xS xS xSxSx这些条件可以确定如下两个方程这些条件可以确定如下两个方程:编辑课件编辑课件42再由三弯矩基本方程再由三弯矩基本方程,可得可得1111222211116262 (2.23)6226nnnnnnnndMdMdMMd 周期样条周期样条编辑课件编辑课件43在在MatlabMatlab中数据点称之为断点。如果三次样条插值没有边中数据点称之为断点。如果三次样条插值没有边界条件，最常用的方法，就是采用非扭结（界条件，最常用的方法，就是采用非扭结（not-a-knotnot-a-knot）条件。这个条件强

27、迫条件。这个条件强迫第第1 1个和第个和第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,conds),pp=csape(x0,y0,conds),pp=csape(x0,y0,conds,valcon

28、ds)pp=csape(x0,y0,conds,valconds)，y=ppval(pp,x)y=ppval(pp,x)。其中其中x0,y0 x0,y0是已知数据点，是已知数据点，x x是插值点，是插值点，y y是插值点的函是插值点的函数值。数值。对于三次样条插值，我们提倡使用函数对于三次样条插值，我们提倡使用函数csapecsape，csapecsape的返的返回值是回值是pppp形式，要求插值点的近似函数值，必须调用函数形式，要求插值点的近似函数值，必须调用函数ppvalppval。MATLAB中三次样条函数法实现中三次样条函数法实现编辑课件编辑课件44pp=csape(x0,y0,con

29、ds,valconds)pp=csape(x0,y0,conds,valconds)condsconds指定插值的边界条件，其值可为：指定插值的边界条件，其值可为：complete complete 边界为一阶导数，一阶导数的值在边界为一阶导数，一阶导数的值在valcondsvalconds参数中给出。参数中给出。not-a-knot not-a-knot 非扭结条件非扭结条件periodic periodic 周期条件周期条件second second 边界为二阶导数，二阶导数的值在边界为二阶导数，二阶导数的值在valcondsvalconds参数中给出，若忽略参数中给出，若忽略valcon

30、dsvalconds参数，二阶导数的参数，二阶导数的缺省值为缺省值为0, 00, 0。MATLAB中三次样条函数法实现中三次样条函数法实现编辑课件编辑课件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代码代码编辑课件编辑课件46第一边界条件的例题第一边界条件的例题pp = form: pp breaks: 1 2 4 5 coefs: 3x4 double pieces: 3 order: 4 dim:

31、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 分段多项式结构与系数矩阵分段多项式结构与系数矩阵编辑课件编辑课件47如下代码求解上述样条问题：如下代码求解上述样条问题：x-4-3-2-101234y00.151.122.362.361.460.490.060y00例例3：自然边界条件例题

32、：自然边界条件例题编辑课件编辑课件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分段多项式结构：分段多项式结构：编辑课件编辑课件50 pp

33、.coefs= 0.161 0.000 -0.161 0 -0.803 0.54256811487482 0.51171207658321 0.000 -0.39373527245950 0.28972754050074 1.344 1.120 0.14922128129602 -0.896 0.74225699558174 2.360 0.540 -0.44381443298969 -0.593 2.360 0.236 -0.348 -1.888 1.460 -0.483 0.36687039764359 -0.73650773195876 0.490 -0.303 0.18578240058910 -0.607 0.000分段多项式的系数矩阵：分段多项式的系数矩阵：编辑课件编辑课件51图像图像编辑课件编辑课件5221 1,1125yxx 例例4：对如下：对如下Runge现象中的函数，求用现象中的函数，求用n分点作节分点作节点的三次样条插值多项式点的三次样条插值多项式s3(x)的图象。的图象。取取n=5,10,15,20等，将区间等分成等，将区间等分成n份，取分点作为插值份，取分点作为插值节点，利用