插值法和数据拟合 (1)

1、 2 2 3 3问题的提出问题的提出 010011( )1,( ),(0,1, ),( )( ),(0,1, ).( )1,( )( )( )niiiinnyf xnx xxyf xinyP xP xyinyf xnxyx yxyxa bR xf xP x插值问题的数学提法：已知函数在个点上的函数值求一个多项式，使其满足即要求该多项式的函数曲线要经过上已知的个点同时在其他上要估计误差 4 4( )( )( )Rxf xP x 同时在其它同时在其它上要估计误差上要估计误差。当当插值问题插值问题n 1 ( )( )( )Rxf xP x n 1 个点个点同时在其它同时在其它上要估计误差上要估计误差

2、。当当 5 5一次插值一次插值 100111(),nPxxyxy当时 ， 求 一 次 多 项 式要 求 通 过两 点 6 6 20011222( ),nP xxyx yxy当时，求二次多项式要求通过三点 7 7 设函数设函数 y=f(x)在区间在区间a,b连续连续, , 给定给定n+1个点个点已知已知 f(xk)=yk (k=0,1,n) ,在函数类中寻找一函数在函数类中寻找一函数g(x) 作为作为 f(x) 的近似表达式的近似表达式, ,使满足使满足这时称这时称 y=f(x)为为被插值函数被插值函数, , g(x) 称为称为插值函数插值函数, xk 称称为为插值节点插值节点, ,式式(4.2

3、)(4.2)称为称为插值条件插值条件, ,寻求插值函数寻求插值函数g(x) 的方法称为的方法称为插值方法插值方法. .一、一、插值问题插值问题a x0 x1 xnb (4.1)g(xk)=f(xk)=yk , k=0,1, ,n (4.2) 8 8 nnnnnnnnnnyxaxaxaayxaxaxaayxaxaxaa22101121211000202010(4.4) 设设 pn(x)=a0+ a1x+an xn ，当满足如下的插值条件当满足如下的插值条件时，即时，即得到关于系数得到关于系数 a0、 a1、 、an 的线性方程组的线性方程组pn(xk) = f(xk) = yk , k = 0,

4、1, ,n (4.3) 9 9其系数行列式为其系数行列式为Vandermonde(范德蒙范德蒙)行列式行列式 njiijnnnnnnnxxxxxxxxxxxxxxV021211020010)(111),(只要只要 ij 就有就有xi xj , ,因而因而 V(x0,x1,xn) 0 ,于是方程组于是方程组(4.4)(4.4)有唯一解有唯一解, ,也就是说也就是说, ,当节点不重合时当节点不重合时, , n+1+1个节点的个节点的插值条件能唯一确定一个插值条件能唯一确定一个n次插值多项式次插值多项式, ,从而有从而有: : 定理定理4.1 给定给定n+1 个互异节点个互异节点 x0 , x1 ,

5、 , xn 上的函数上的函数值值y0 , y1 , , yn ,则满足插值条件则满足插值条件 pn(xk )=f(xk )= yk , k=0, 1, ,n的的n 次插值多项式次插值多项式 pn(x) 存在且唯一。存在且唯一。 1010这样只要求解方程组这样只要求解方程组1. 拉格朗日插值多项式拉格朗日插值多项式2. 牛顿插值多项式牛顿插值多项式 3. 分段线性插值多项式分段线性插值多项式 4. 三次样条插值多项式三次样条插值多项式便可确定插值多项式便可确定插值多项式 pn(x) ，但相对来讲计算较复杂。，但相对来讲计算较复杂。而插值多项式的唯一性保证了无论用什么方法获得满足而插值多项式的唯一

6、性保证了无论用什么方法获得满足插值条件的多项式都是同一个多项式插值条件的多项式都是同一个多项式 pn(x) ，因此可以，因此可以采用其它更简便的方法来确定多项式。下面就介绍几种采用其它更简便的方法来确定多项式。下面就介绍几种常用的方法：常用的方法： nnnnnnnnnnyxaxaxaayxaxaxaayxaxaxaa22101121211000202010(4.4) 1111的函数值的函数值 nyyy,10已知已知 y=f(x) 在在n+1 个点个点bxxxan 10nkyxfxpkkkn, 1 , 0,)()( 构造构造n次多项式次多项式 pn(x) ,使得使得从而得到从而得到 f(x) 的

7、近似计算式的近似计算式 ,),()(baxxpxfn 1212求解求解 p1(x)=a1 x+a0ixiy0 x1x0y1y已知已知使得使得 f(x) p1(x), x x0 , x1.10100101yxxxxyxxxx 根据点斜式得到根据点斜式得到如果令如果令1010)(xxxxxl 0101)(xxxxxl 则称则称 l0(x) , l1(x)为一次插值多项式的为一次插值多项式的基函数基函数。这时。这时：并称其为一次并称其为一次Lagrange插值多项式。插值多项式。 0)(, 1)(1000 xlxl1)(, 0)(1101 xlxlf(x) p1(x)=y0l0(x)+y1 l1(x

8、)1010010( )()yyP xyxxxx 1313例例 已知已知 的两个采样点（的两个采样点（100,10）和（）和（121,11）,求求 。115xy x 100 121y=f(x) 10 11)()(0010101xxxxyyyxp 101001011)(yxxxxyxxxxxp 1115 121115 100( )1151011 10.714285100 121121 100p xg10 1 , lg20 1.3010lg12( )lg( )lg(10)1(20)1.3010102011.3010201101( )(20)( )(10)10 201020

9、1010f xx f xx ffxxyyxxlxxl xx例：已知，利用插值一次多项式求的近似值。解：，设，则插值基本多项式为：， 151510 01 1111.3010( )( )( )(20)(10)101011.3010(12)(12 20)(12 10) 1.06021010lg12 lg10 lg20 lg12 1.0792).P xy l xyl xxxP于是，拉格朗日型一次插值多项式为：故即由和两个值的线性插值得到，且具有两位有效数字（精确解 1616求解求解 p2(x)=a2x2+a1 x+a0使得使得 f(x) p2(x), x x0 , x2.关于二次多项式的构造采用如下方

10、法：令关于二次多项式的构造采用如下方法：令已知已知ixiy0 x1x0y1y2x2y并由插值条件并由插值条件得到得到)(20100 xxxxyA )(21011xxxxyB )(12022xxxxyC p2(x)=A(x-x1)(x-x2)+B(x-x0)(x-x2)+C(x-x0)(x-x1)p2(x0)=y0 , p2(x1)=y1 , p2(x2)=y2 1717于是得到于是得到则有则有 jijixlijij, 0, 1)( f(x) p2(x)=y0l0(x)+y1 l1(x)+ y2 l2(x)()()(2010210 xxxxxxxxxl 如果令如果令)()()(2101201xx

11、xxxxxxxl )()()(1202102xxxxxxxxxl 并称其为二次并称其为二次Lagrange 插值多项式。插值多项式。 则称则称 l0(x) , l1(x)，l(x)为二次插值多项式的为二次插值多项式的基函数基函数。这时。这时：0201122012010210122021()()()()()()( )()()()()()()x xx xx xx xx x x xp xyyyxx xxxxxxxxxx 1818例例已知已知 , , ,10100 11121 115y解解: : 这里这里x0 0=100=100，y0 0=10=10，x1 1=121=121，y1 1=11, =11

12、, x2 2=144=144，y2 2=12=12，利用抛物线插值公式，利用抛物线插值公式 利用抛物线插值求利用抛物线插值求144122115(121)(144)(100)(144)(115)1011(100 121)(100 144)(121 100)(121 144)(100)(121)1210.72275(144 100)(144 121)xxxxxLxx 19192例 ：已知ixlgiiyxixlgiiyx10152011.17611.3010lg12利用此三值的二次插值多项式求的近似值012012101520(15)(20)1( )1520(10 15)(1020)50(10)(20

13、)1( )1020(15 10)(1520)25(10)(15)1( )1015(20 10)(20 15)50 xxxxxlxxxxxl xxxxxlxxx 解：设，则 2020ixlgiiyx20 01 12 221( )( )( )( )2015501.17611.301010201015255011.1761(12)1220 12 1512 10 122050251.301012 10 12 151.076650lg12lg121.07923P xy lxy l xy lxxxxxxxP故所以：利用三个点进行抛物线插值得到的的值，与精确值相比，具有 位有效数字。 2121已知已知n+1

14、组离散数据组离散数据ixiy0 x1x0y1ynxny按照二次按照二次Lagrange插值多项式的构造方法，令：插值多项式的构造方法，令：将插值条件将插值条件 pn( x0 )= y0 代入，得到：代入，得到：)()(0201000nxxxxxxyA 同理，由插值条件同理，由插值条件 pn( x1 )= y1 ，得到：得到：)()(1210111nxxxxxxyA 012102011( )()()()()()()()()()nnnnnpxA xxxxxxA xxxxxxA xxxxxx 2222由插值条件由插值条件 pn( xn )= yn ,得得)()(110 nnnnnnxxxxxxyA将

15、将)()(0201000nxxxxxxyA )()(1210111nxxxxxxyA )()(,110 nnnnnnxxxxxxyA代入下式：代入下式：得到：得到：012102011( )()()()()()()()()()nnnnnp xA xxxxxxA xxxxxxA xxxxxx1202010102010121011011()()()()()()( )()()()()()()()()()()()()nnnnnnnnnnnxxxxxxxxxxxxp xyyxxxxxxxxxxxxxxxxxxyxxxxxx 2323也就是也就是011011()()()()()()()()jjnjjjjjj

16、jnxxxxxxxxyxxxxxxxx1202010102010121()()()()()()( )()()()()()()nnnnnx x x xx xx xx xx xp xyyxx xxxxxxxxxx011011()()()()()()nnnnnnxxxxxxyxxxxxx 0110011()()()()()()()()njjnjjjjjjjjnxxxxxxxxyxxxxxxxx 从而就得到了在区间从而就得到了在区间 a, ,b 上的关于函数上的关于函数 f(x) 的的n次多项式近似计算式：次多项式近似计算式： ( )( ), , nfxpxxa b 242401234012340(2

17、,0)(4,3)(6,5)(8,4)(10,1)24681003541(4)(6)(8)(10)( )(24)(26)(28)(2 10)1(4)(6)(8)(10)384xxxxxyyyyyxxxxlxxxxx例3：求过点的拉格朗日型插值多项式。解：用4次插值多项式对5个点插值：，； 2525123(2 ) (6 ) (8 ) (1 0 )()( 42 ) ( 46 ) ( 48 ) ( 41 0 )1(2 ) (6 ) (8 ) (1 0 )9 6(2 ) (4 ) (8 ) (1 0 )()( 62 ) ( 64 ) ( 68 ) ( 61 0 )1(2 ) (4 ) (8 ) (1 0

18、 )6 4(2 ) (4 ) (6 ) (1 0 )()( 82 ) ( 84 ) ( 86 ) ( 81 0 )1xxxxlxxxxxxxxxlxxxxxxxxxlx，4(2 ) (4 ) (6 ) (1 0 )9 6(2 ) (4 ) (6 ) (8 )()(1 02 ) (1 04 ) (1 06 ) (1 08 )1(2 ) (4 ) (6 ) (8 )3 8 4xxxxxxxxlxxxxx，； 262640 01 12 23 34 4( )( )( )( )( )( )3(2)(6)(8)(10)965(2)(4)(8)(10)644(2)(4)(6)(10)961(2)(4)(6)

19、(8)384P xy lxy l xy lxy l xy lxxxxxxxxxxxxxxxxx 所以 2727拉格朗日插值多项式的截断误差拉格朗日插值多项式的截断误差 ,( )( ),( )( )( )( )( )( )( )0nnnniniinia bP xf xR xR xf xP xxxR xf xP x我们在上用多项式来近似替代函数其截断误差记作，。当 在插值结点 上时，估计截断误差。 2828( )(1)01(1)011( )( ) , ( )( , )( ) , 1( )( )( )(1)!( , ),( )()()()nnnnnnnnnyf xnyfxa byfxa baxxxb

20、P xnxa bR xfxna bxxxxxxx 定理 ：设函数的 阶导数在上连续，在上存在；插值结点为，是次拉格朗日插值多项式；则对任意有：其中 2929对于线性插值，其余项为对于二次插值，其余项为 303012212lg121lg10lg20lg12P 121.0602,lg121.0792e1.0792 1.06020.0190( )lg1( ),10,20ln101|( )|0.043ln10 101|( )(12 10)(1220)| 8 0.0430.3442f xxfxxxff 分析例，例 中计算的截断误差在例中，用和计算，估计误差：当时 3131243242lg10 lg15l

21、g20lg12.P 121.0766,e1.0792 1.07660.00262( )8.686 10ln101|(12)| |( )(12 10)(12 15)(1220)|3!18.686 102 3 80.006956fxxRf 在例 中，用，和计算估计误差：故 3232作业1.已知函数表求抛物插值多项式及其余项。2已知插值基函数，证明：当m7时时,舍入误差亦会增大舍入误差亦会增大.)()!()()()()()(xwnfxLxfxRnnnn1 考考虑虑 6464,01n+1niiiiii+1i : a,bn+1a= x xxx =by = f(x )(i=1,2,.,n),x ,1 j

22、x = y ,(i=1,2,.,n)(2)x ,x,(i=0,1,2,.,n-1)xxx定义设在区间上给定个节点及节点上的函数值作一个插值函数（ ） 使其满足（） （ ）在每个小区间上，（ ）都是线性插值函数。称函数（ ）为a,b上关于数据（)(0,1,2,., )iyin的分段线性插值函数。分段分段线性线性插值插值 6565由Lagrange线性插值公式容易写出 的分段表达式也可以通过构造基函数的方法来求 。首先构造一组基函数 ,每个 都满足2） 在每个小区间 上都是线性函数。这组函数称为分段线性插值函数。其表达式也可以写为) x(111)=(0,1,2,., )i+1iiiiiii+1ii

23、x- xxxxyyxxxinx - xxx() x( )(0,1,2,., )il x in( )il x01) ()( ,0,1,2,., )1ijjil xi jnji()ilx1 ,(0,1,2,., )iix xin10100101( ),0,xxxxlxxx xxx x 6666111111011,( ) ,(1,2,.,1)0,)(,iiiiiiiiiiiiinxxxxxxxxxl xxx xinxxxx xxx11101,( )0,nnnnnnnxxxxxxxlxxx x 6767类似于lagarange插值多项式的构造，函数就是所求的分段插值函数。例已知函数y=f(x)=1/(

24、1+x2),在区间0,5上取等距节点xi=0+i(i=0,1,5)。求分段线性插值函数，并由此计算f(4.5)的近似值。解 节点处的函数值见下表0()()niiixy lxxi012345yi1.000000.500000.200000.100000.058820.03846 6868解解 分段插值基函数为分段插值基函数为5 , 1 (0 1 , 01)(0 xxxxl 1, 15 , 00 1,) 1(, 11)(iixiixixiixixxli)4 , 005 , 44)(xxxxln 6969所以，分段插值函数为所以，分段插值函数为)(10000. 0)(20000. 0)(50000.

25、 0)()(3210 xlxlxlxlx)(03846. 0)(05882. 054xlxl)5 . 4(03846. 0)5 . 4(05882. 0)5 . 4()5 . 4(54llf5 . 003846. 05 . 005882. 004864. 004706. 0)5 . 4(f比较，结果是比较精确的。比较，结果是比较精确的。 7070分段线性函数的优缺点：1）简便易行，而且当节点加密时，分段线性插值误差变小，收敛性可以得到保证。2）由于每个小区间的插值函数只依赖于本段的节点值，因而每个节点只影响相邻的一两个小区间，不会扩大误差，从而保证插值结果的稳定性。3）由于分段插值函数仅在a,

26、b上连续。一般的，在节点处插值函数不可微，因此不满足工程上光滑度的问题。 7171例 对sinx进行线性插值示例解在编辑窗口输入clear;x=0:10;y=sin(x);xi=0:0.25:10;yi=interp1(x,y,xi);plot(x,y,o,xi,yi) 7272 4.5 Hermite插值法插值法Lagrange插值虽然构造比较简单，但插值曲线只是在节点处与原函数吻合，若还要求在节点处两者相切，即导数值相等，使之与被插函数的”密切”程度更好，这就要用到带导数的插值.0101( ),nnf xax xxbfff设在节点处的函数值为值函数上的具有一阶导数的插的在区间为设,)()(

27、baxfxP处必须满足在节点即nxxxxP,)(10()()iiiP xf xf()()iiiP xfxf)(,)()1(一阶光滑度上具有一阶导数在若要求baxP-(1)0,1,in 7373处必须满足在节点即nxxxxP,)(10()()iiiP xf xf()()iiiP xfxf)(,)(,)2(阶光滑度阶导数上具有在若要求同样mmbaxP个待定的系数可以解出22n次的多项式可以是最高次数为因此12)(nxP次多项式作为插值函数两个节点就可以用3112()()iiiP xfxf()()()()()mmmiiiPxfxfni, 1 ,0-(2)22n 共个方程 7474定义1. 称满足(1

28、)或(2)式的插值问题为Hermite插值，称满足(1)或(2)式的插值多项式P(x)为Hermite插值多项式，记为 , 为多项式次数.,( )(),kkHermiteHxkRungek一般次插值多项式的次数 如果太高会影响收敛性和稳定性现象因此 不宜太大.一、两点三次Hermite插值先考虑只有两个节点的插值问题0101( ),f xxxff设在节点处的函数值为0101,xxff在节点处的的一阶导数值为作为插值函数多项式次两个节点最高可以用)(33xHHermite( )kHxk 7575应满足插值条件)(3xH300()Hxf311()Hxf300()Hxf311()Hxf示应用四个插值

29、基函数表)(3xH3 ,2 ,1 ,0),()(3ixhxHi的插值基函数为设)()()()()(332211003xhaxhaxhaxhaxH希望插值系数与Lagrange插值一样简单重新假设300110011( )( )( )( )( )Hxfxfxfxfx 7676其中1)(00 x0)(00 x1)(00 x300110011( )( )( )( )( )Hxfxfxfxfx 0)(10 x0)(01x1)(11x0)(10 x0)(01x0)(11x0)(00 x0)(10 x0)(01 x0)(11 x0)(10 x0)(01 x1)(11 x可知即可假设的二重零点是,)(01xx

30、)()()(210baxxxx1)(00 x0)(00 x由 7777可得310)(2xxa3100210)(2)(1xxxxxb)()()(210baxxxx21)(xx 310)(2xxx3100210)(2)(1xxxxx21021)()(xxxx102xxx10021xxx01021xxxx2101xxxx)()(21(201xlxlLagrange插值基函数 7878)(1x)()(21(210 xlxl类似可得)(0 x)()(200 xlxx)(1x)()(211xlxx10121xxxx2010 xxxx0 xx 2101xxxx2010 xxxx1xx )(0 x01021x

31、xxx2101xxxx)()(21(201xlxl即将以上结果代入 7979300110011( )( )( )( )( )Hxfxfxfxfx得两个节点的三次Hermite插值公式300110011( )( )( )( )( )Hxfxfxfxfx2101(12 ( )( )flxlx2000()( )fxxlx2111()( )f xxlx2010(12 ( )( )fl xlx110112xxfxx2010 xxxx00fxx2101xxxx2010 xxxx11fxx001012xxfxx2101xxxx 8080二、两点三次Hermite插值的余项两点三次Hermite插值的误差为)

32、()()(33xHxfxR0)()()(33iiixHxfxR0)()()(33iiixHxfxR1 , 0i因此可设的二重零点均为,)(,310 xRxx21203)()()(xxxxxKxR待定其中)(xK 818121203)()()()()(xtxtxKtHtft构造辅助函数0)()()()()(21203xxxxxKxHxfxiiiii1 , 0i0)()()()()(21203xxxxxKxHxfx均是二重零点个零点至少有因此5)(t连续使用4次Rolle定理，可得，,10 xx至少存在一点使得0)()4( 82820)(! 4)()()4()4(xKf即! 4)()()4(fxK

33、所以,两点三次Hermite插值的余项为2120)4(3)()(! 4)()(xxxxfxR10 xx以上分析都能成立吗？上述余项公式成立上存在时在当,)(10)4(xxxf 838311( )( ),0,1,( )( )nniiiiiiHfinHfxxxxyyHermite满足插值条件:Hermite值问题的解插值法定理1 存在且唯一.(22)221210( )( )( )() , , (2(2 !)nnnnjja bf xHfxnRxxx Hermite插值多项式的余项为Hermite插2:值法定理 8484例.1)2(,0)1(21)(3)2(,2)1(21)(ffxfffxf处的导数值

34、为，在节点处的函数值为，在节点已知.7 . 1 , 5 . 1)(,)(处的函数值在及的两点三次插值多项式求xxfxf解:2, 110 xx012,3ff010,1ff 300110011( )( )( )( )( )Hxfxfxfxfx110112xxfxx2010 xxxx00fxx2101xxxx2010 xxxx11fxx001012xxfxx2101xxxx 8585)2(213x21x21x2 x)1(212x22x)(3xH91713323xxx)5 . 1(f)5 . 1(3H625. 2)7 . 1(f)7 . 1(3H931. 2作为多项式插值作为多项式插值,三次已是较高的

35、次数，次数再高就有三次已是较高的次数，次数再高就有可能发生可能发生Runge现象现象,因此，对有因此，对有n+1节点的插值问题，节点的插值问题，我们可以使用分段两点三次我们可以使用分段两点三次Hermite插值插值 8686例 按下表求hermite插值多项式解 将表中数据代入hermite插值多项式 xi01yi01yi3922222320110H ()(12)()0(12)()11001011010(0)()3(1)()90110(32)3(1)9(1)10123xxxxxxxxxx xx xxxxxx 8787Heemite插值的matlab实现function f = Hermite(

36、x,y,y1,x0)syms t;f = 0.0;if(length(x) = length(y) if(length(y) = length(y1) n = length(x); else disp(y和y的导数的维数不相等); return; endelse disp(x和y的维数不相等！ ); return;end for i=1:n h = 1.0; a = 0.0; for j=1:n if( j = i) h = h*(t-x(j)2/(x(i)-x(j)2); a = a + 1/(x(i)-x(j); end end f = f + h*(x(i)-t)*(2*a*y(i)-y

37、1(i)+y(i);endif (nargin=4) f=subs(f,t,x0);end 8888例 hermite插值应用实例，根据下表计算hermite插值多项式，并计算当x=1.56时的y值。解：在matlab编辑窗口输入x=1.0:0.2:1.8;y=1 1.0954 1.1832 1.2649 1.3416;y1=0.5 0.4564 0.4226 0.3953 0.3727;f=Hermite(x,y,y1)运行结果的x1.01.21.41.61.8y11.4095 1.18321.24691.3416y 0.50000.4564 0.42260.39530.3727 8989f

38、 =390625/576*(t-6/5)2*(t-7/5)2*(t-8/5)2*(t-9/5)2*(-61/3+64/3*t)+390625/36*(t-1)2*(t-7/5)2*(t-8/5)2*(t-9/5)2*(-114418258618928321/10995116277760000+337232823972231/35184372088832*t)+390625/16*(t-1)2*(t-6/5)2*(t-8/5)2*(t-9/5)2*(6660373488918492043/11258999068426240000+7612884810106783/1801439850948198

39、4*t)+390625/36*(t-1)2*(t-6/5)2*(t-7/5)2*(t-9/5)2*(384779664999124623/21990232555520000-2855713758717179/281474976710656*t)+390625/576*(t-1)2*(t-6/5)2*(t-7/5)2*(t-8/5)2*(8968626627620006931/175921860444160000-7762319875242775/281474976710656*t) 9090下面计算x=1.56时的值f=Hermite(x,y,y1,1.56)运行结果f = 1.2490 例4

40、.16 p142页解t=0.1 0.5 1 1.5 2 2.5 3;y=0.95 0.84 0.86 1.06 1.5 0.72 1.9;y1=1 1.5 2 2.5 3 3.5 4;yy=hermite(t,y,y1,1.8)t1=0.1:0.01:3;yy1=hermite(t,y,y1,t1);plot(t,y,o,t,y1,*,t1,yy1) 9191运行结果yy = 1.3298 92924.6三次样条插值1.三次差值和三次样条插值的matlab命令对于三次差值和三次样条插值，matlab中都设有专用命令三次插值命令的调用格式如下：yk=pchip(x,y,xk)三次样条插值命令的调

41、用格式如下：yk=spline(x,y,xk)参数x,y, xk, yk的意义及要求与interp1中的完全相同，插值效果与interp1中的参数”method”分别选用pchip和spline等价。 9393例 在y=f(x)的函数关系中，当x=-6.2 -5 -3 -1.7 0 1 3 5 6.2 7.5 8.2时，y=1 0 -1 0 1 0 -1 0 0.5 -1 0。用“线性插值”、“最近插值”分别求出当xi=-6.2:0.2:8.2时所对应的函数值，并画出曲线图。解：x=-6.2 -5 -3 -1.7 0 1 3 5 6.2 7.5 8.2;y=1 0 -1 0 1 0 -1 0

42、0.5 -1 0 ;xi=-6.2:0.2:8.2;y1=interp1(x,y,xi,linear);y2=interp1(x,y,xi,nearest);plot(x,y,ro,xi,y1,xi,y2,-);legend(,); 9494 9595例 在y=f(x)的函数关系中，当x=0 2 3 4 5 6时，y=-4 1 -2 -4 -4 -5 。用分段三次插值和三次样条插值分别求出当xi=-0.3：0.3：6.8时所对应的函数值，并画出曲线。解x=0 2 3 4 5 6;y=-4 1 -2 -4 -4 -5 ;xi=-0.3:0.3:6.8;y1=interp1(x,y,xi,pchi

43、p) y1=pchip(x,y,xi)y2=interp1(x,y,xi,spline) y2=spline(x,y,xi)plot(x,y,ro,xi,y1,r,xi,y2,-,linewidth,2);legend(样本点,分段三次插值,三次样条插值);grid; 9696 9797例对离散的分布在y=exsinx函数曲线上的数据点进行样条插值计算。解 x=0 2 4 5 8 12 12.8 17.2 19.9 20;y=exp(x).*sin(x);xx=0:0.5:20;yy=spline(x,y,xx)plot(x,y,o,xx,yy) 9898得到的结果得到的结果yy共有共有41个

44、数据点个数据点 9999( )( ) , S xf xa b三三次次样样则则条条称称为为在在上上的的插插值值函函数数q 三次样条插值函数三次样条插值函数(2)( ),( ),( ) , S xS xSxa b 都都在在区区间间上上连连续续1(1)( ),iiS xxx 在在每每个个小小区区间间上上都都是是三三次次多多项项式式( ) , S xa b则则称称为为区区间间上上的的三三次次样样条条函函数数(3)( )S x而而三三次次样样条条函函数数满满足足(),0,1,iif xy in(),0,1,iiS xy in01, , naxxxba b 为为区区间间的的一一个个分分割割( ) , :S

45、 xa b如如果果函函数数在在区区间间上上满满足足条条件件 100100q 三次样条插值原理三次样条插值原理在在n n个小区间构造个小区间构造S(xS(x) )，共有，共有n n个三次多项式，需确定个三次多项式，需确定4n4n个参数个参数在所有节点上在所有节点上n+1个方程个方程在除端点外的节点上在除端点外的节点上3(n-1)个方程个方程这样就得到这样就得到4n - 2个方程，为保证待定参数的唯一性，还差两个个方程，为保证待定参数的唯一性，还差两个方程。为此，常用的方法是对边界节点除函数值外附加要求，方程。为此，常用的方法是对边界节点除函数值外附加要求，这就是所谓的边界条件。根据实际问题的不同

46、，三次样条插值这就是所谓的边界条件。根据实际问题的不同，三次样条插值常用到下列常用到下列三类边界条件三类边界条件。( ),0,1,iiS xy in1, 1,)()(limniyxSxSiixxi1, 1),()(limnixSxSixxi1, 1,)()(lim niMxSxSiixxi 101101 m边界条件：边界条件： 即给定端点处的一阶导数值。即给定端点处的一阶导数值。 M边界条件：边界条件： 即给定端点处的二阶导数值。即给定端点处的二阶导数值。 周期边界条件周期边界条件： 当当y = g(x)是以是以 b-a= x0 - xn 为周期的周期函数时，要求为周期的周期函数时，要求S(x

47、)也是周期函数，故端点处要满足也是周期函数，故端点处要满足S(a + 0) = S(b - 0), S(a + 0) = S(b - 0)，此条件称为周期条件此条件称为周期条件。00() , ()nnS xyS xy00() , ()nnSxySxy 102102三次样条插值的matlab实现实例例 在区间的两个端点处给出的条件，称为边界条件。给定两端点处的二阶导数值S(a)=y0,S(b)=yn,然后根据下面的数据点求出其三次样条插值多项式，并计算当x=3.5时y的值。解：在matlab命令窗口输入x=1:8;y=0.84 0.91 0.14 -0.76 -0.96 -0.28 0.66,0

48、.99f,f0=ThrSample1(x,y,0.54,-0.15,3.5)x12345678y0.840.910.14-0.76-0.96-0.280.660.99y0.5403-0.1455 103103程序运行结果如下f =(7/25*t-7/10)*(t-4)2+(-171/25+38/25*t)*(t-3)2+(-71823/72775*t+215469/72775)*(4-t)2-(-37788/14555+9447/14555*t)*(t-3)2f0 = -0.3522 104104程序运行时调用的三次样条函数如下：function f,f0 = ThrSample1 (x,y,

49、y_1, y_N,x0)syms t;f = 0.0;f0 = 0.0;if(length(x) = length(y) n = length(x);else disp(x和y的维数不相等！); return;end %维数检查for i=1:n if(x(i)=x0) index = i; break; endend %找到x0所在区间for i=1:n if(x(i)=x0) index = i; break; endend %找到x0所在区间 A = diag(2*ones(1,n); %求解m的系数矩阵u = zeros(n-2,1);lamda = zeros(n-1,1);c =

50、zeros(n,1);for i=2:n-1 u(i-1) = (x(i)-x(i-1)/(x(i+1)-x(i-1); lamda(i) = (x(i+1)-x(i)/(x(i+1)-x(i-1); c(i) = 3*lamda(i)*(y(i)-y(i-1)/(x(i)-x(i-1)+ . 3*u(i-1)*(y(i+1)-y(i)/(x(i+1)-x(i); A(i, i+1) = u(i-1); A(i, i-1) = lamda(i); %形成系数矩阵及向量cendc(1) = 2*y_1;c(n) = 2*y_N;m = followup(A,c); %用追赶法求解方程组h = x

51、(index+1) - x(index); %x0所在区间长度f = y(index)*(2*(t-x(index)+h)*(t-x(index+1)2/h/h/h + . y(index+1)*(2*(x(index+1)-t)+h)*(t-x(index)2/h/h/h + . m(index)*(t-x(index)*(x(index+1)-t)2/h/h - . m(index+1)*(x(index+1)-t)*(t-x(index)2/h/h; %x0所在区间的插值函数f0 = subs(f,t,x0); %x0处的插值 105105其中的followup函数如下：function

52、 x=followup(A,b)n=rank(A);for i=1:n if A(i,i)=0 disp (error:对角线有元素0！) return; endendd=ones(n,1);a=ones(n-1,1);c=ones(n-1);for i=1:n-1 a(i,1)=A(i+1,i); c(i,1)=A(i,i+1); d(i,1)=A(i,i);endd(n,1)=A(n,n);for i=2:n d(i,1)=d(i,1)-(a(i-1,1)/d(i-1,1)*c(i-1,1); b(i,1)=b(i,1)-(a(i-1,1)/d(i-1,1)*b(i-1,1);endx(n

53、,1)=b(n,1)/d(n,1);for i=n-1:-1:1 x(i,1)=(b(i,1)-c(i,1)*x(i+1,1)/d(i,1);end 106106例 在区间的两个端点处给出的条件，称为边界条件。给定两端点处的二阶导数值S(a)=y0,S(b)=yn,然后根据下面的数据点求出其三次样条插值多项式，并计算当x=3.5时y的值。解：x=1:8;y=0.84 0.91 0.14 -0.76 -0.96 -0.28 0.66,0.99f,f0=ThrSample2(x,y,0.84,-0.99,3.5)x12345678y0.840.910.14-0.76-0.96-0.280.660.

54、99y-0.8415-0.9894 107107运行结果如下f =(7/25*t-7/10)*(t-4)2+(-171/25+38/25*t)*(t-3)2+(-23849/23288*t+71547/23288)*(4-t)2-(-5730726096775419/2251799813685248+5730726096775419/9007199254740992*t)*(t-3)2f0 = -0.3585 108108程序运行调用的三次样条函数如下function f,f0 = ThrSample2 (x,y,y2_1, y2_N,x0)syms t;f = 0.0;f0 = 0.0;if

55、(length(x) = length(y) n = length(x);else disp(x和y的维数不相等！); return;end %维数检查for i=1:n if(x(i)=x0) index = i; break; endend %找到x0所在区间for i=1:n if(x(i)=x0) index = i; break; endend %找到x0所在区间 A = diag(2*ones(1,n); %求解m的系数矩阵A(1,2) = 1;A(n,n-1) = 1;u = zeros(n-2,1);lamda = zeros(n-1,1);c = zeros(n,1);for

56、 i=2:n-1 u(i-1) = (x(i)-x(i-1)/(x(i+1)-x(i-1); lamda(i) = (x(i+1)-x(i)/(x(i+1)-x(i-1); c(i) = 3*lamda(i)*(y(i)-y(i-1)/(x(i)-x(i-1)+ . 3*u(i-1)*(y(i+1)-y(i)/(x(i+1)-x(i); A(i, i+1) = u(i-1); A(i, i-1) = lamda(i); %形成系数矩阵及向量cendc(1) = 3*(y(2)-y(1)/(x(2)-x(1)-(x(2)-x(1)*y2_1/2;c(n) = 3*(y(n)-y(n-1)/(x(

57、n)-x(n-1)-(x(n)-x(n-1)*y2_N/2;m = followup(A,c); %用追赶法求解方程组h = x(index+1) - x(index); %x0所在区间长度f = y(index)*(2*(t-x(index)+h)*(t-x(index+1)2/h/h/h + . y(index+1)*(2*(x(index+1)-t)+h)*(t-x(index)2/h/h/h + . m(index)*(t-x(index)*(x(index+1)-t)2/h/h - . m(index+1)*(x(index+1)-t)*(t-x(index)2/h/h; %x0所在

58、区间的插值函数f0 = subs(f,t,x0); %x0处的插值 109109例4-25根据表中的数据，用程序求三次样条插值函数解 用程序进行计算x=0,1,2,3;y=0,0.5,2,1.5;dx0=0.2;dxn=-1; s=csfit(x,y,dx0,dxn)回车得到s = 0.4800 -0.1800 0.2000 0 -1.0400 1.2600 1.2800 0.5000 0.6800 -1.8600 0.6800 2.0000 x0123f(x)00.521.5f(x)0.2-1 1101103213223230.480.180.2 ,0,11.041.261.280.5,1,

59、20.681.860.682,2,3SxxxxSxxxxSxxxx 111111程序调用时的函数如下：function S=csfit(x,y,dx0,dxn)% x,y分别为n个节点的横坐标所组成的向量及纵坐标组成的向量% dx0,dn分别为S的导数在x0,xn处的值n=length(x)-1;h=diff(x);d=diff(y)./h;a=h(2:n-1);b=2*(h(1:n-1)+h(2:n);c=h(2:n);u=6*diff(d);b(1)=b(1)-h(1)/2;u(1)=u(1)-3*(d(1)-dx0);b(n-1)=b(n-1)-h(n)/2;u(n-1)=u(n-1)-

60、3*(dxn-d(n);for k=2:n-1 temp=a(k-1)/b(k-1); b(k)=b(k)-temp*c(k-1); u(k)=u(k)-temp*u(k-1);endm(n)=u(n-1)/b(n-1);for k=n-2:-1:1 m(k+1)=(u(k)-c(k)*m(k+2)/b(k);endm(1)=3*(d(1)-dx0)/h(1)-m(2)/2;m(n+1)=3*(dxn-d(n)/h(n)-m(n)/2;for k=0:n-1 S(k+1,1)=(m(k+2)-m(k+1)/(6*h(k+1); S(k+1,2)=m(k+1)/2; S(k+1,3)=d(k+1