三次樣條擬合典型實(shí)例

1、1設(shè)計(jì)目的、要求 對(duì)龍格函數(shù)在區(qū)間-1,1上取的等距節(jié)點(diǎn)，分別作多項(xiàng)式插值、三次樣條插值和三次曲線擬合，畫出及各逼近函數(shù)的圖形，比較各結(jié)果。2設(shè)計(jì)原理（1） 多項(xiàng)式插值：利用拉格朗日多項(xiàng)式插值的方法，其主要原理是拉格朗日多項(xiàng)式，即：表示待插值函數(shù)的個(gè)節(jié)點(diǎn)，其中；（2） 三次樣條插值：三次樣條插值有三種方法，在本例中，我們選擇第一邊界條件下的樣條插值，即兩端一階導(dǎo)數(shù)已知的插值方法： （3）三次曲線擬合：本題中采用最小二乘法的三次多項(xiàng)式擬合。最小二乘擬合是利用已知的數(shù)據(jù)得出一條直線或者曲線，使之在坐標(biāo)系上與已知數(shù)據(jù)之間的距離的平方和最小。在本題中，n= 10，故有11個(gè)點(diǎn)，以這11個(gè)點(diǎn)的和值為已知

2、數(shù)據(jù)，進(jìn)行三次多項(xiàng)式擬合，設(shè)該多項(xiàng)式為，該擬合曲線只需的值最小即可。3采用軟件、設(shè)備 計(jì)算機(jī)、matlab軟件4設(shè)計(jì)內(nèi)容1、 多項(xiàng)式插值：在區(qū)間上取的等距節(jié)點(diǎn)，帶入拉格朗日插值多項(xiàng)式中，求出各個(gè)節(jié)點(diǎn)的插值，并利用matlab軟件建立m函數(shù)，畫出其圖形。在matlab中建立一個(gè)lagrange.m文件，里面代碼如下：%lagrange 函數(shù)function y=lagrange(x0,y0,x)n=length(x0);m=length(x);for i=1:m z=x(i); s=0.0; for k=1:n p=1.0; for j=1:n if j=k p=p*(z-x0(j)/(x0(k

3、)-x0(j); end end s=p*y0(k)+s; end y(i)=s;end建立一個(gè)polynomial.m文件，用于多項(xiàng)式插值的實(shí)現(xiàn)，代碼如下：%lagrange插值x=-1:0.2:1;y=1./(1+25*x.2);x0=-1:0.02:1;y0=lagrange(x,y,x0);y1=1./(1+25*x0.2);plot(x0,y0,'-r')%插值曲線hold on%原曲線plot(x0,y1,'-b')運(yùn)行duoxiangshi.m文件，得到如下圖形：2、 三次樣條插值：所謂三次樣條插值多項(xiàng)式是一種分段函數(shù)，它在節(jié)點(diǎn) EMBED Equ

4、ation.DSMT4 分成的每個(gè)小區(qū)間上是3次多項(xiàng)式，其在此區(qū)間上的表達(dá)式如下：因此，只要確定了的值，就確定了整個(gè)表達(dá)式，的計(jì)算方法如下：令：則滿足如下個(gè)方程：對(duì)于第一種邊界條件下有如果令那么解就可以為 求函數(shù)的二階導(dǎo)數(shù)：>> syms x>> f=sym(1/(1+25*x2) f =1/(1+25*x2) >> diff(f) ans = -(50*x)/(25*x2 + 1)2將函數(shù)的兩個(gè)端點(diǎn)，代入上面的式子中：f(-1)= 0.0740f(1)=-0.0740 求出從-1到1的n=10的等距節(jié)點(diǎn)，對(duì)應(yīng)的x，y值 對(duì)應(yīng)m文件代碼如下：for x=-1

5、:0.2:1 y=1/(1+25*x2)endy =得出x=-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1y=0.0385 0.0588 0.1 0.2 0.5 1 0.5 0.2 0.1 0.0588 0.0385 編寫m文件Three_Spline.mx=linspace(-1,1,11);y=1./(1+25*x.2);m,p=scyt1(x,y,0.0740,-0.0740);hold onx0=-1:0.01:1;y0=1./(1+25*x0.2);plot(x0,y0,'-b')得到如下圖像：.其中藍(lán)色曲線為原圖，紅色曲線為擬合

6、后的圖像。3、 三次曲線擬合：這里我們使用最小二乘法的3次擬合建立一個(gè)Three_fitting .m文件，代碼如下：%主要代碼x=-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1;y=0.0385 0.0588 0.1 0.2 0.5 1 0.5 0.2 0.1 0.0588 0.0385;a=polyfit(x,y,3);x1=-1:0.01:1;y1=a(4)+a(3)*x1+a(2)*x1.2+a(1)*x1.3;x0=-1:0.01:1;y0=1./(1+25*x0.2)%原曲線plot(x0,y0,'-r')hold on%三次擬

7、合曲線plot(x1,y1,'-b')上圖中，藍(lán)色部分為三次擬合曲線，紅色部分為原曲線6結(jié)果分析拉格朗日插值的優(yōu)點(diǎn)是對(duì)于某一區(qū)域，不限于被估計(jì)點(diǎn)周圍，公式簡單易實(shí)施。一般認(rèn)為n的次數(shù)越高，逼近的精度就越好，但在本題中，對(duì)龍格函數(shù)，中間部分插值效果比較好，而對(duì)于兩端，插值結(jié)果是非常不理想的，即龍格現(xiàn)象。樣條函數(shù)可以給出光滑的插值曲線，從本題中就能體現(xiàn)出來。從以上圖形可以看出，三次樣條插值的圖形是比較逼近于原圖的，收斂性相對(duì)而言是非常好的，但在本題中，僅將原區(qū)間分成10個(gè)等距區(qū)間，因此，逼近效果還不是特別理想，當(dāng)我們將n增大時(shí)，插值后的曲線越逼近于原曲線。總的來說，三次樣條插值的穩(wěn)

8、定性比較好，收斂性比較強(qiáng)。在這三種方法中，三次曲線擬合的效果是最差的，所得的圖形與原曲線差距甚遠(yuǎn)。最小二乘法中，并不要求擬合后的曲線經(jīng)過所有已知的點(diǎn)，只需要擬合多項(xiàng)式上的點(diǎn)在某種標(biāo)準(zhǔn)上與定點(diǎn)之間的差距最小即可，因此與原曲線的逼近程度是最差的。最小二乘法的多項(xiàng)式擬合只適用于多項(xiàng)式，而本題中的函數(shù)并不是一個(gè)多項(xiàng)式，因此，不建議使用最小二乘法擬合。參考文獻(xiàn):1 李慶揚(yáng) 王能超等.數(shù)值分析M.清華大學(xué)出版社2 吳振遠(yuǎn).科學(xué)計(jì)算實(shí)驗(yàn)指導(dǎo)書 基于MATLAB數(shù)值分析M. 中國地質(zhì)大學(xué)出版社3 宋葉志. MATLAB數(shù)值分析與應(yīng)用M. 機(jī)械工業(yè)出版社 , 2009.07附錄三次樣條插值主要代碼：functi

9、on m,p=scyt1(x,y,df0,dfn)n=length(x);r=ones(n-1,1);u=ones(n-1,1);d=ones(n,1);r(1)=1;d(1)=6*(y(2)-y(1)/(x(2)-x(1)-df0)/(x(2)-x(1);u(n-1)=1;d(n)=6*(dfn-(y(n)-y(n-1)/(x(n)-x(n-1)/(x(n)-x(n-1);for k=2:n-1 u(k-1)=(x(k)-x(k-1)/(x(k+1)-x(k-1); r(k)=(x(k+1)-x(k)/(x(k+1)-x(k-1); d(k)=6*(y(k+1)-y(k)/(x(k+1)-x

10、(k)-(y(k)-y(k-1)/(x(k)-x(k-1)/(x(k+1)-x(k-1);endA=eye(n,n)*2;for k=1:n-1 A(k,k+1)=r(k); A(k+1,k)=u(k);endm=Ad;ft=d(1);syms tfor k=1:n-1 %求s(x)即插值多項(xiàng)式 p(k,1)=m(k)/(6*(x(k+1)-x(k); p(k,2)=m(k+1)/(6*(x(k+1)-x(k); p(k,3)=(y(k)-m(k)*(x(k+1)-x(k)2/6)/(x(k+1)-x(k); p(k,4)=(y(k+1)-m(k+1)*(x(k+1)-x(k)2/6)/(x(k+1)-x(k); sx(k)=p(k,1)*(x(k+1)-t)3+p(k,2)*(t-x(k)3+p(k,3)*(x(k+1)-t)+p(k,4)*(t-x(k);endkmax=1000;xt=linspace(x(1),x(n),kmax);for i=1:n-1 %出點(diǎn)xt對(duì)應(yīng)的y值 for k=1:kmax if x(i)<=xt(k)&xt(k)<=x(i+1) fx(k)=subs(sx(i),xt(k); e