关于三次样条插值函数的学习报告(研究生)资料

学习报告——三次样条函数插值问题的讨论班级：数学二班学号：152111033姓名：刘楠楠样条函数：由一些按照某种光滑条件分段拼接起来的多项式组成的函数；最常用的样条函数为三次样条函数，即由三次多项式组成，满足处处有二阶连续导数。

一、三次样条函数的定义：对插值区间[,]a b 进行划分，设节点011n n a x x x x b -=<<<<=，若函数2()[,]s x c a b ∈在每个小区间1[,]i i x x +上是三次多项式，则称其为三次样条函数。

如果同时满足()()i i s x f x = (0,1,2)i n =，则称()s x 为()f x 在[,]a b 上的三次样条函数。

二、三次样条函数的确定：由定义可设：1012121(),[,](),[,]()(),[,]n n n s x x x x s x x x x s x s x x x x -∈⎧⎪∈⎪=⎨⎪⎪∈⎩其中()k s x 为1[,]k k x x -上的三次多项式，且满足11(),()k k k k k k s x y s x y --== (1,2,,k n =由2()[,]s x C a b ∈可得：''''''()(),()(),k k k k s x s x s x s x -+-+== 有''1()(),k k k k s x s x -++=''''1()(),(1,2,,1)k k k k s x s x k n -++==-，已知每个()k s x 均为三次多项式，有四个待定系数，所以共有4n 个待定系数，需要4n 个方程才能求解。

前面已经得到22(1)42n n n +-=-个方程，因此要唯一确定三次插值函数，还要附加2个条件，一般上，实际问题通常对样条函数在端点处的状态有要求，即所谓的边界条件。

1、第一类边界条件：给定函数在端点处的一阶导数，即 ''''00(),()n n s x f s x f ==2、第二类边界条件：给定函数在端点处的二阶导数，即''''''''00(),()n n s x f s x f ==特别地，当''''0()()0n s x s x ==时，称为自然边界条件，此时的样条函数称为自然样条函数。

3、第三类边界条件：设()f x 是周期函数，并设0n x x -是一个周期，于是()s x 满足''''''00()(),()()n n s x s x s x s x == 三、三次样条函数的分析计算：设出()s x 在每个节点处的二阶导数值，即''(),(0,1,2,,)j j s x M j n == 考虑区间1[,]j j x x -，在此区间上，()()j s x s x =是三次多项式，故''()j s x 为线性函数，且 ''''''''111()(),()()j j j j j j j j s x s x M s x s x M ---==== 利用线性插值公式，可得''()j s x 的表达式： 1''11(),j j j j j j j j jjx x x x s x M M h x x h h -----=+=-积分两次后可得()j s x 的表达式： 331112()()()66j j j j j jjx x x x s x M M c x c h h ----=+++将插值条件11(),()j j j j j j s x y s x y --==代入可确定积分常数1c 和2c ，整理上式得332211111()()()6666j j j j j j jj j j j j j jjj j x x x x M h x x M h x x s x M M y y h h h h -----⎛⎫⎛⎫----=++-+- ⎪ ⎪ ⎪ ⎪⎝⎭⎝⎭只需确定01,,,n M M M 即可给出()s x 的表达式。

()s x 在各个节点处的一阶导数存在''()()j j s x s x -=即有''1()()j j j j s x s x -++=对()j s x 求导得：22111'1()()()226j j j j j j jj j j jjjx x x x y y M M s x M M h h h h --------=-++-则1111111()()2626j j j j j j j jj j j j j j jj h h y y h h y y M M M M M M h h -+++-++----+=---+整理得关于1,j j M M -和1j M +的方程：112j j j j j j M M M d μλ-+++=其中111111111,,166[,,]j j j j j j j j j j j j j j j j j j j j j j h h h h h h y y y y d f x x x h h h h μλμλ++++--+++⎧==+=⎪++⎪⎨⎛⎫--⎪=-= ⎪⎪ ⎪+⎝⎭⎩1,2,,1j n =-共1n -个方程，附加边界条件后，补充两个方程，便可确定1n +个未知量01,,,n M M M 。

第一类边界条件：''''00(),()n ns x y s x y ==直接代入()j s x 的一阶导数表达式即得 ''1100101011()()26,26n n n n n n nny y y y y y h h M M d M M d h h --⎛⎫⎛⎫---- ⎪ ⎪⎝⎭⎝⎭+=≡+=≡与前面的1n -个方程联立得到1n +阶线性方程组：001111222211112122212n n n n n n M d M d M d M d M d μλμλμλ----⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦类似的有第二类和第三类边界条件得到的方程组：''1111022222222''11112222n n n n n n n n n M d y M d M d M d y λμμλμλμλ--------⎡⎤-⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥-⎣⎦⎣⎦⎣⎦111222211112222n n n n nnn n M d M d M d M d λμλμλλμ----⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦又知上述三个方程组的系数矩阵严格对角占优，故矩阵可逆，则方程组的解存在且唯一。

四、三次样条函数的一致收敛性和稳定性：1、收敛性：设2()[,],()f x C a b S x ∈是以(0,1,,)k x k n =为节点，满足任意边界条件得三次样条插值函数，设10101,max,min i i i i i i n i n h x x h h h δ+≤≤-≤≤-=-==，则当hc δ≤<∞时，()S x 和'()S x 在[,]a b 上一致收敛到()f x 和'()f x ，即''00lim()(),lim ()(),[,]h h S x f x S x f x x a b →→==∈从而我们可以知道：当插值节点越多时，三次样条插值函数就越接近被插函数，也就不会产生龙格现象中不收敛的情况。

2、稳定性：在三次样条插值法中，插值节点处函数值的波动，只对这个节点两边的分段有影响，而对离该点较远的分段的影响会逐渐变小，因此样条插值法具有较好的稳定性。

五、三次样条插值函数的算例与实际应用：例1、求满足下面函数表所给出的插值条件的三次样条插值函数。

解：已知031,0M M ==，算得030λμ==，111112λ==+，211112λ==+， 11112μλ=-=，22112μλ=-=，02d =，1230d d d ===，所以由 第二类边界条件可得：0123202122120122120020M M M M ⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦计算M 及样条系数得：01231,415,115,0M M M M ==-==x 0 1 2 3 y0 0 0 0 ''y10,10,20,31,11,22690,12,1990,790,215C C C C C =-==-==-1,32,12,22,3118,290,130,190C C C C ==-==-所以三次样条插值函数为32323219126,[0,1]90290127()(1)(1)(1),[1,2]181590112(2)(2)(2),[2,3]903090x x x x s x x x x x x x x x ⎧-+-∈⎪⎪⎪=---+-∈⎨⎪⎪--+---∈⎪⎩例2、已知一组数据点，求解满足这些数据点的三次样条插值函数解：由于数据点含有多位小数，不便计算，所以本题采用Matlab 编程，源代码见附录，结果为2.8620,2.7541,2.4130,2.2421,2.1450g =2.0000 1.00000.6429 2.00000.35710.4000 2.00000.60000.5714 2.00000.42861.0000 2.0000Q ⎡⎤⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎣⎦ 0.9697,0.9227,0.7992,0.7424,0.7013m =插值函数为s1=264523075327265143/56294995342131200*X^2-23087199381998953/112589990684262400*X+40023205577025431/112589990684262400-17634871688484217/2814749767106560*X^3s2=-160081743506404901/60798594969501696*X^2+404209705972252i x 0.250.30.390.450.53iy0.50000.54770.62450.67080.7280727/202661983231672320*X+429142243010323951/3166593487994880000+257361898089296225/136796838681378816*X^3s3=1437374409830143/13510798882111488*X^2+10427488839800859/11258999068426240*X+12358231431982259943/45035996273704960000-3495912536773825/7599824371187712*X^3s4=26732501105874704913/720575940379279360000+38626753769033575/18014398509481984*X^3-245666153971053021/72057594037927936*X^2+3614707928905781673/1441151880758558720*X0.250.30.350.40.450.50.550.60.50.550.60.650.70.75xy图像eval例3、已知汽车门曲线型值点的数据如下： i x 0 12 3 4 5 6 7 8 9 10i y 2.51 3.30 4.04 4.70 5.22 5.54 5.78 5.40 5.575.70 5.80求车门的三次样条插值函数，并打印出(4.5)s ，端点条件为''0100.8,0.2y y ==解：源代码见附录，插值函数为3232323232(1)0.00851400.00148600.80000 2.5100(2)0.00445790.0136540.81217 2.5059(3)0.00365440.0184750.82181 2.4995(4)0.0409240.316960.18448 3.5058(5)0.10735 1.46246s t t t s t t t s t t t s t t t s t t =--++=--++=--++=-+-+=-+3232323232.9328 5.9839(6)0.26848 4.175221.25540.996(7)0.426598.336153.813109.14(8)0.26786 6.247448.272129.06(9)0.054872 1.498313.69436.184(10)0.058376 1.592914.545t s t t t s t t t s t t t s t t t s t t -=-+-+=-+-=-+-+=-+-=-+38.738t -插值点处值为 (4.5) 5.3833s =0123456789102.533.544.555.56附录例2、源代码：x=[0.25 0.3 0.39 0.45 0.53]y=[0.5000 0.5477 0.6245 0.6708 0.7280]n=length(x);for i=1:n-1h(i)=x(i+1)-x(i);endfor i=1:n-2k(i)=h(i+1)/(h(i)+h(i+1));u(i)=h(i)/(h(i)+h(i+1));endfor i=1:n-2gl(i)=3*(u(i)*(y(i+2)-y(i+1))/h(i+1)+k(i)*(y(i+1)-y(i))/h(i)); endg0=3*(y(2)-y(1))/h(1);g00=3*(y(n)-y(n-1))/h(n-1);g=[g0 gl g00];g=transpose(g)k1=[k 1];u1=[1 u];Q=2*eye(5)+diag(u1,1)+diag(k1,-1)m=transpose(Q\g)syms X;for i=1:n-1p1(i)=(1+2*(X-x(i))/h(i))*((X-x(i+1))/h(i))^2*y(i);p2(i)=(1-2*(X-x(i+1))/h(i))*((X-x(i))/h(i))^2*y(i+1);p3(i)=(X-x(i))*((X-x(i+1))/h(i))^2*m(i);p4(i)=(X-x(i+1))*((X-x(i))/h(i))^2*m(i+1);p(i)=p1(i)+p2(i)+p3(i)+p4(i);p(i)=simple(p(i));ends1=p(1)s2=p(2)s3=p(3)s4=p(4)for k=1:4for z=x(k):0.001:x(k+1)q=eval(subs(p(k),'X','z'));plot(z,q,'b')hold onendendgrid onlegend('eval')title('图像')xlabel('x')ylabel('y')例3、源代码：function []=spline3(X,Y,dY,x0,m)N=size(X,2);s0=dY(1);sN=dY(2);interval=0.025;disp('x0为插值点')x0h=zeros(1,N-1);for i=1:N-1h(1,i)=X(i+1)-X(i);endd(1,1)=6*((Y(1,2)-Y(1,1))/h(1,1)-s0)/h(1,1);d(N,1)=6*(sN-(Y(1,N)-Y(1,N-1))/h(1,N-1))/h(1,N-1);for i=2:N-1d(i,1)=6*((Y(1,i+1)-Y(1,i))/h(1,i)-(Y(1,i)-Y(1,i-1))/h(1,i-1))/(h(1,i)+h(1,i-1)); endmu=zeros(1,N-1);md=zeros(1,N-1);md(1,N-1)=1;mu(1,1)=1;for i=1:N-2u=h(1,i+1)/(h(1,i)+h(1,i+1));mu(1,i+1)=u;md(1,i)=1-u;endp(1,1)=2;q(1,1)=mu(1,1)/2;for i=2:N-1p(1,i)=2-md(1,i-1)*q(1,i-1);q(1,i)=mu(1,i)/p(1,i);endp(1,N)=2-md(1,N-1)*q(1,N-1);y=zeros(1,N);y(1,1)=d(1)/2;for i=2:N三次样条函数插值问题的讨论- 10 - - 10 - y(1,i)=(d(i)-md(1,i-1)*y(1,i-1))/p(1,i);endx=zeros(1,N);x(1,N)=y(1,N);for i=N-1:-1:1x(1,i)=y(1,i)-q(1,i)*x(1,i+1);endfprintf('M 为三对角方程的解\n');M=x;fprintf('\n');syms t;digits(m);for i=1:N-1pp(i)=M(i)*(X(i+1)-t)^3/(6*h(i))+M(i+1)*(t-X(i))^3/(6*h(i))+(Y(i)-M(i)*h(i)^2/6)*(X(i+1)-t)/h(i)+(Y(i+1)-M(i+1)*h(i)^2/6)*(t-X(i))/h(i);pp(i)=simplify(pp(i));coeff=sym2poly(pp(i));if length(coeff)~=4tt=coeff(1:3);coeff(1:4)=0;coeff(2:4)=tt;endif x0>X(i)&x0<X(i+1)L=i;y0=coeff(1)*x0^3+coeff(2)*x0^2+coeff(3)*x0+coeff(4);endval=X(i):interval:X(i+1);for k=1:length(val)fval(k)=coeff(1)*val(k)^3+coeff(2)*val(k)^2+coeff(3)*val(k)+coeff(4);endif mod(i,2)==1plot(val,fval,'r+')elseplot(val,fval,'b.')endhold onclear val fvalans=sym(coeff,'d');ans=poly2sym(ans,'t');fprintf('在区间[%f,%f]内\n',X(i),X(i+1));fprintf('三次样条函数S(%d)=',i);pretty(ans);endfprintf('x0所在区间为[%f,%f]\n',X(L),X(L+1));fprintf('函数在插值点x0=%f 的值为\n',x0);y0。