1、【精品文檔】如有侵權(quán)，請(qǐng)聯(lián)系網(wǎng)站刪除，僅供學(xué)習(xí)與交流數(shù)值分析課程設(shè)計(jì).精品文檔.2013/2014第一學(xué)期數(shù)值分析課程設(shè)計(jì)設(shè)計(jì)題目： 插值算法總結(jié)與Matlab編程實(shí)現(xiàn) 信息與計(jì)算科學(xué) 專業(yè) 學(xué)號(hào) 112086202 姓名 胡波 學(xué)號(hào) 112086237 姓名 許世開 學(xué)號(hào) 112086235 姓名 郭海波 成績(jī) 指導(dǎo)老師 郭尊光 摘 要本論文首先分析用插值多項(xiàng)式表示復(fù)雜函數(shù)f（x）的必要性，然后運(yùn)用插值的方法根據(jù)給定的函數(shù)表構(gòu)造出一個(gè)既能反映函數(shù)f (x)的特征，又便于計(jì)算的簡(jiǎn)單函數(shù)p(x)，使得p(x)近似f (x)。文章通過對(duì)拉格朗日插值法存在的缺點(diǎn)的分析，引入了牛頓插值法、埃爾米特插值

2、法、分段低次插值法及三次樣條插值法，給出了每種插值算法的算法結(jié)構(gòu)、具體的程序MATLAB代碼及誤差估計(jì)的方法，并且對(duì)每種插值方法的優(yōu)缺點(diǎn)進(jìn)行了分析。關(guān)鍵詞：拉格朗日插值，牛頓插值，埃爾米特插值，分段低次插值，三次樣條插值問題提出 設(shè)f (x)的函數(shù)表為x =2,4,6,8,10;y =9.8,9.2,8.1,6.4,3.8，y=0.96,0.87,0.58,0.45,0.25,g（x）=，(-5x5)。請(qǐng)根據(jù)給出的函數(shù)及函數(shù)表解下面的問題。 (1)證明插值多項(xiàng)式的存在唯一性，以x ,y為節(jié)點(diǎn)計(jì)算其插值多項(xiàng)式。（2）用基函數(shù)構(gòu)造拉格朗日插值多項(xiàng)式，以x ,y為節(jié)點(diǎn)做拉格朗日插值曲線。（3）分析牛

3、頓插值多項(xiàng)式的算法及余項(xiàng)，利用牛頓插值多項(xiàng)式，求f (6.596)并估計(jì)誤差。 (4）埃爾米特插值多項(xiàng)式的構(gòu)造及余項(xiàng)，利用埃爾米特插值多項(xiàng)式，求f (6.596)。（5）在g（x）=，(-5x5)的條件下，分析分段線性插值的必要性（龍格現(xiàn)象）及誤差。（6）在g（x）=，(-5x5)的條件下，建立三次樣條函數(shù)，通過觀察圖象，得出結(jié)論“三次樣條插值s(x)較十次拉格朗日插值L(x)能更好地逼近被插函數(shù)f (x)”。1 插值多項(xiàng)式1.1插值多項(xiàng)式定義設(shè)函數(shù)y=f (x)在區(qū)間a,b上有定義，且已知在點(diǎn)上的值，若存在一個(gè)n次的多項(xiàng)式p(x)，使得 （1.1） 成立，則稱p(x)為f (x)的n次拉格朗

4、日插值多項(xiàng)式，點(diǎn)稱為插值節(jié)點(diǎn)，包含插值節(jié)點(diǎn)的區(qū)間a,b稱為插值區(qū)間。注： 1、上述定義中的(1.1)式通常稱為插值條件，它是n+1個(gè)條件。 2、插值多項(xiàng)式的幾何意義如下 求曲線y=p(x)，使其通過給定的n+1個(gè)點(diǎn)() ( i=0,1, 2,n )，并用它近似已知曲線y=f (x)。1.2插值多項(xiàng)式存在唯一性與計(jì)算首先我們證明插值多項(xiàng)式存在且唯一。由插值條件(1.1),有 + + = （1.2）式（1.2）是一個(gè)關(guān)于未知量的線性方程組，其系數(shù)矩陣的行列式是一個(gè)范德蒙行列式由于插值節(jié)點(diǎn)是互異的，所以V0。從而線性方程組(1.2)有唯一解，p(x)存在且唯一。以上的證明表明：由n+1個(gè)插值條件可唯

5、一地確定一個(gè)n個(gè)拉氏插值多項(xiàng)式。下面我們?cè)趩栴}1的條件下計(jì)算以x=2,4,6,8,10,y=9.8,9.2,8.1,6.4,3.8為插值節(jié)點(diǎn)的插值多項(xiàng)式。MATLAB代碼如下：format shortx=2,4,6,8,10'y=9.8,9.2,8.1,6.4,3.8'A=x.0,x.1,x.2,x.3B=Ay運(yùn)行結(jié)果為：B =10.1600,-0.1476,-0.0071,-0.0042' 這表明: 。2 用基函數(shù)構(gòu)造拉格朗日插值多項(xiàng)式2.1 拉格朗日插值多項(xiàng)式的構(gòu)造利用求解線性方程組(1.2)求拉格朗日插值多項(xiàng)式不便計(jì)算機(jī)實(shí)現(xiàn)，不易編程。在實(shí)際應(yīng)用中，一般采用基函數(shù)

6、法來構(gòu)造。下面，我們介紹這一方法，這里，我們將n次拉格朗日插值多項(xiàng)式用來記。設(shè),其中是待定的n次多項(xiàng)式。首先構(gòu)造個(gè)插值節(jié)點(diǎn)上的插值基函數(shù)，對(duì)任一點(diǎn)所對(duì)應(yīng)的插值基函數(shù)，由于在所有取零值，因此有因子。又因是一個(gè)次數(shù)不超過的多項(xiàng)式，所以只可能相差一個(gè)常數(shù)因子，固可表示成：利用得：于是因此滿足 的插值多項(xiàng)式可表示為：從而次拉格朗日插值多項(xiàng)式為：2.2拉格朗日插值多項(xiàng)式在MATLAB中的實(shí)現(xiàn)及使用在MATLAB中構(gòu)造拉格朗日插值函數(shù)代碼如下：function yy=lagelangrichazhi(x,y,xx)m=length(x);n=length(y);if m=n;end s=0; for i=

7、1:n t=ones(1,length(xx); for j=1:n if j=i,t=t.*(xx-x(j)/(x(i)-x(j); end end s=s+t*y(i); end yy=s;利用該函數(shù)做題一中x，y的插值曲線代碼如下：clearx=2,4,6,8,10;y=9.8,9.2,8.1,6.4,3.8;xi=2:0.1:10;yi=lagelangrichazhi(x,y,xi);plot(x,y,'o',xi,yi,'k');title('拉格朗日插值')得到拉格朗日插值曲線如圖2.1所示： 圖2.1 拉格朗日插值曲線3牛頓多項(xiàng)式

8、lagrange插值多項(xiàng)式作為一種計(jì)算方案，公式簡(jiǎn)潔，做理論分析也方便。其缺點(diǎn)是，當(dāng)節(jié)點(diǎn)改變時(shí)，公式需要重建，計(jì)算量大；如果還要根據(jù)精度要求，選取合適的節(jié)點(diǎn)和插值多項(xiàng)式的次數(shù)，則只好逐次計(jì)算出 , 等，并做誤差試算，才可以做到，這當(dāng)然是不理想的。為此，人們從改進(jìn)插值多項(xiàng)式的形式入手，給出另一種風(fēng)格的插值公式，這就是Newton(牛頓)插值公式。3.1均差概念及其性質(zhì)設(shè)函數(shù)在互異節(jié)點(diǎn)上的值為 , 等，定義：(1) 在上的1階均差為 (2) 在上的2階均差為 （3）遞推地，在上的階均差為同時(shí)規(guī)定在上的零階均差為 。均差具有的性質(zhì)如下：性質(zhì)1 階均差可以表示成個(gè)函數(shù)值的線性組合，即或記為 性質(zhì)2 均

9、差對(duì)其節(jié)點(diǎn)是對(duì)稱的，即節(jié)點(diǎn)按任意順序排列，其值不變。如這個(gè)性質(zhì)稱為對(duì)稱性。性質(zhì)3 如果是的次多項(xiàng)式，則其1階均差是的次多項(xiàng)式，且由此遞推可得階均差為零階均差，階均差為零。差商表的構(gòu)造如表1.3所示：表3.1 差商表 零階均差1階均差2階均差 3階均差 利用MATLAB計(jì)算題1中x，y的各階均差代碼如下：x=2,4,6,8,10'y=9.8,9.2,8.1,6.4,3.8'n=max(size(x); %計(jì)算節(jié)點(diǎn)個(gè)數(shù)z=zeros(n,n+1);%建一個(gè)n 行n+1 列的零矩陣z(:,1)=x; %將x 賦給z 的第1 列z(:,2)=y;%將y 賦給z 的第2 列for j=1

10、:n-1for i=j+1:nz(i,j+2)=(y(i)-y(i-1)/(x(i)-x(i-j); %計(jì)算j 階差商endfor i=j:ny(i)=z(i,j+2); %將j 階差商值賦給向量y 的對(duì)應(yīng)分量endendz計(jì)算結(jié)果為：3.2牛頓插值多項(xiàng)式的算法及余項(xiàng)由均差的定義可知對(duì)上述第二式兩端乘以，第三式兩端乘以, 最后一式兩端乘以，然后由后一式的左端代入前一式的右端(第二項(xiàng))，即可得 顯然，是次多項(xiàng)式，又有，故有，從而滿足條件，這就說明是關(guān)于的次插值多項(xiàng)式，通常稱它為Newton插值公式，為插值余項(xiàng)。在MATLAB中使用牛頓插值法求得題1中以x，y為插值節(jié)點(diǎn)的條件下f(6.596)的值

11、并估計(jì)誤差的代碼如下：format longx=2,4,6,8,10;y=9.8,9.2,8.1,6.4,3.8;n=max(size(x);%計(jì)算差商for j=1:1:n-1for i=n:-1:j+1y(i)=(y(i)-y(i-1)/(x(i)-x(i-j);endend%計(jì)算牛頓插值多項(xiàng)式的值z(mì)=6.596;%需求其函數(shù)值的點(diǎn)p=y(n);for i=n-1:-1:1p=y(i)+p*(z-x(i);endfprintf('f(%8.5f)=%8.5fn',z,p);%誤差估計(jì)r=y(n);for i=1:nr=r*(z-x(i);endr=abs(r);fprint

12、f('截?cái)嗾`差R(%8.5f)=%15.14fn',z,r);計(jì)算結(jié)果如下：4 埃爾米特插值多項(xiàng)式 拉氏插值多項(xiàng)式的還存在著一個(gè)缺陷：在插值節(jié)點(diǎn)處，被插函數(shù)曲線與插值多項(xiàng)式曲線的切線方向可能會(huì)不一致。不少實(shí)際的插值問題不但要求在節(jié)點(diǎn)上函數(shù)值相等，而且還要求對(duì)應(yīng)的導(dǎo)數(shù)值也相等，甚至要求高階導(dǎo)數(shù)也相等，滿足這種要求的插值多項(xiàng)式稱為埃爾米特（Hermite）插值多項(xiàng)式。4.1埃爾米特插值多項(xiàng)式定義設(shè)在節(jié)點(diǎn)上，若存在一個(gè)2n+1次多項(xiàng)式H(x)，使得 成立，則稱H(x)為f (x)的2n+1次埃爾米特插值多項(xiàng)式。4.2埃爾米特插值多項(xiàng)式的構(gòu)造及余項(xiàng)設(shè)計(jì)2n+2個(gè)插值基函數(shù)，每一個(gè)基函數(shù)

13、都是2n+1次多項(xiàng)式，并且滿足以下條件于是滿足條件(2.5)的插值多項(xiàng)式H(x)可以表示成為下面的問題就是求滿足條件(2.6)的基函數(shù)。a a 利用拉氏插值基函數(shù)以及條件(2.6)，可知的形式應(yīng)為又整理得解出由于于是類似地，可推導(dǎo)出。仿照拉格朗日余項(xiàng)的證明方法，若f (x)在(a,b)內(nèi)的2n+2階導(dǎo)數(shù)存在，則其插值余項(xiàng)為在MATLAB中構(gòu)造埃爾米特函數(shù)代碼如下： function f=Hermite(x,y,y_1,x0) syms t; f=0.0; if(length(x)=length(y) if(length(y)=length(y_1) n=length(x); else disp

14、('y和y的導(dǎo)數(shù)的維數(shù)不相等！'); return; end else disp('y和y的導(dǎo)數(shù)的維數(shù)不相等！'); 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_1(i)+y(i); if(i=n) if(nargin=4) f=subs(f,'t',x0); else f=vpa(f,6); end end end在MATLA

15、B中計(jì)算以x，y為插值節(jié)點(diǎn)，且對(duì)應(yīng)的一階導(dǎo)數(shù)為y的埃爾米特多項(xiàng)式，并計(jì)算在x=6.596點(diǎn)處函數(shù)的值代碼如下：clearx=2,4,6,8,10;y=9.8,9.2,8.1,6.4,3.8;y1=0.96,0.87,0.58,0.45,0.25;f=Hermite(x,y,y1)f1=Hermite(x,y,y1,6.596)計(jì)算結(jié)果為：f=0.000108507*(8.53667*t - 24.9467)*(t - 6.0)2*(t - 10.0)2*(t - 2.0)2*(t - 8.0)2 - 0.000108507*(4.88333*t - 45.4667)*(t - 6.0)2*(t

16、 - 10.0)2*(t - 2.0)2*(t - 4.0)2 + 0.00000678168*(21.3767*t - 32.9533)*(t - 6.0)2*(t - 10.0)2*(t - 4.0)2*(t - 8.0)2 - 0.00000678168*(7.66667*t - 80.4667)*(t - 6.0)2*(t - 2.0)2*(t - 4.0)2*(t - 8.0)2 + 0.000244141*(t - 10.0)2*(t - 2.0)2*(t - 4.0)2*(t - 8.0)2*(0.58*t + 4.62)f1 = 8.19065 分段線性插值和分段三次hermi

17、te插值5.1龍格現(xiàn)象 lagrange插值方法使根據(jù)區(qū)間a,b上給出的節(jié)點(diǎn)構(gòu)造插值多項(xiàng)式近似f(x)。一般認(rèn)為的次數(shù)n越高，逼近f(x)的效果越好，龍格現(xiàn)象表明 時(shí)，不一定收斂到f(x)。我們以函數(shù)為例研究龍格現(xiàn)象，MATLAB代碼如下：x=-5:1:5.1;y=1./(1+x.2);plot(x,y,'r');n=max(size(x);hold onz=-5:0.005:5; m=max(size(z); for k=1:mL(k)=0;for i=1:nt=1;for j=1:n if j=it=t*(z(k)-x(j)/(x(i)-x(j);endendL(k)=L(

18、k)+t*y(i);endendplot(z,L,'b')hold offtitle('龍格現(xiàn)象 ')圖5.1 龍格現(xiàn)象這說明高次插值多項(xiàng)式近似的效果并不好，因而通常并不采用高次插值，而采用分段低次插值。5.2分段線性插值分段線性插值就是通過插值點(diǎn)用折線段連接起來逼近。設(shè)已知節(jié)點(diǎn)a=x0 <x1 <<xn =b上的函數(shù)值yk =f(xk)(k=0,1,n)，已知函數(shù)y= f(x) 在這些插值結(jié)點(diǎn)的函數(shù)值為yk =f(xk)(k=0,1,n)求一個(gè)分段函數(shù)Ih(x), 使其滿足：  (1) Ih(xk )=yk ，(k=0,1

19、,n);  (2) 在每個(gè)區(qū)間xk ,xk+1 上,Ih (x)是個(gè)一次函數(shù)。易知,Ih(x)是個(gè)折線函數(shù), 在每個(gè)區(qū)間xk ,xk+1 上,(k=0,1,n) 于是, Ih (x)在a,b上是連續(xù)的，但其一階導(dǎo)數(shù)是不連續(xù)的。分段線性函數(shù)的基函數(shù)就是一階lagrange插值多項(xiàng)式。其誤差為（ 其中）在matalb中一為插值函數(shù)interp1默認(rèn)的就是線性插值。在MATLAB中做函數(shù)g（x）=，(-5x5)插值節(jié)點(diǎn)數(shù)不同的插值曲線，可以看到分段低次插值不會(huì)產(chǎn)生龍格現(xiàn)象，代碼如下：clear%輸入插值區(qū)間的等分?jǐn)?shù)disp('給出插值區(qū)間的等分?jǐn)?shù)n(5 的倍數(shù))'

20、)n=input('n = ');%作被插函數(shù)圖象u=-5:(10/200):5;v=2*u./(1+u.2);plot(u,v,'r');hold on%固化圖形屏幕%給出插值條件x=-5:(10/n):5;y=2*x./(1+x.2);%取需用分段線性插值函數(shù)計(jì)算其函數(shù)值的點(diǎn)zz=-5:10/(2*n):5;m=max(size(z);%計(jì)算分段線性插值函數(shù)在z處的值，并作圖for k=1:mi=1;while i<=nif z(k)>=x(i) && z(k)<=x(i+1)Ih(k)=y(i)*(z(k)-x(i+1)/

21、(x(i)-x(i+1)+y(i+1)*(z(k)-x(i)/(x(i+1)-x(i);breakelsei=i+1;endendendplot(z,Ih,'b')hold off%釋放固化的圖形屏幕輸入等分?jǐn)?shù)為5時(shí)，拉格朗日插值曲線如圖5.2所示，圖5.2 等分?jǐn)?shù)為5時(shí)的拉格朗日插值曲線輸入等分?jǐn)?shù)為800時(shí)，拉格朗日插值曲線如圖5.3所示， 圖5.3 等分?jǐn)?shù)為800時(shí)的拉格朗日插值曲線6 三次樣條插值多項(xiàng)式通過平面上在n+1個(gè)不同的已知點(diǎn)()，可以作一條光滑的n次代數(shù)曲線，但當(dāng)點(diǎn)很多時(shí)，除了插值多項(xiàng)式的次數(shù)高，計(jì)算復(fù)雜外，從逼近的效果來看，隨著插值多項(xiàng)式的次數(shù)增高，逼近的效果

22、并不是隨之越來越理想，相反還會(huì)產(chǎn)生龍格現(xiàn)象。采用分段線性插值函數(shù)去逼近被插函數(shù)，能保持較好的逼近效果，但其曲線在節(jié)點(diǎn)處僅僅只能保證連續(xù)，卻并不是光滑的。在生產(chǎn)和科學(xué)實(shí)驗(yàn)中，對(duì)所構(gòu)造的插值曲線，既要求簡(jiǎn)單，又要求在連接處比較光滑。即所構(gòu)造的分段插值函數(shù)在分段上要求多項(xiàng)式次數(shù)低，而在分點(diǎn)（節(jié)點(diǎn)）上還具有連續(xù)的低階導(dǎo)數(shù)，滿足這樣條件的插值函數(shù)，我們稱為樣條插值。6.1三次樣條插值多項(xiàng)式算法組織所謂三次樣條插值多項(xiàng)式是一種分段函數(shù)，它在節(jié)點(diǎn)分成的每個(gè)小區(qū)間上是3次多項(xiàng)式，其在此區(qū)間上的表達(dá)式如下：因此，只要確定了的值，就確定了整個(gè)表達(dá)式，的計(jì)算方法如下：令：則滿足如下n-1個(gè)方程：對(duì)于第一種邊界條件下

23、有如果令那么解就可以為 6. 2三次樣條插值多項(xiàng)式算例分析做以x，y為插值節(jié)點(diǎn)的三次樣條插值曲線代碼如下：x = 2,4,6,8,10;y = 9.8,9.2,8.1,6.4,3.8;xx=2:0.2:10;yy = spline(x,y,xx);plot(x,y,'o',xx,yy)title('三次樣條插值曲線')grid on所得結(jié)果如圖6.1所示：圖6.1 三次樣條插值曲線6.3三次樣條插值函數(shù)與拉格朗日插值多項(xiàng)式的效果比較為顯示三次樣條插值函數(shù)與拉格朗日插值多項(xiàng)式的效果比較，我們用MATLAB分別做g（x）=，(-5x5)的函數(shù)曲線、拉格朗日插值曲線、

24、三次樣條插值曲線，代碼如下：clearformat short%作y=2*x/(1+x2)的圖象u=-5:(10/200):5;v=2*u./(1+u.2);plot(u,v,'r')%紅線表示函數(shù)圖象hold on%給出插值條件x=-5:1:5;y=2*x./(1+x.2);y_bar(1)=152/(26)2; %在-5,5的一階導(dǎo)數(shù)值y_bar(2)=152/(26)2;%計(jì)算每個(gè)小區(qū)間長(zhǎng)度for i=1:10h(i)=x(i+1)-x(i);endfor i=2:10lmd(i)=h(i)/(h(i-1)+h(i);mu(i)=1-lmd(i);b(i)=6*(y(i+

25、1)-y(i)/h(i)-(y(i)-y(i-1)/h(i-1)/(h(i)+h(i-1);endlmd(1)=1;mu(11)=1;b(1)=(6/h(1)*(y(2)-y(1)/h(1)-y_bar(1);b(11)=(6/h(10)*(y_bar(2)-(y(11)-y(10)/h(10);%用追趕法解三對(duì)角方程,求矩mbata(1)=lmd(1)/2;for i=2:10bata(i)=lmd(i)/(2-mu(i)*bata(i-1);endu(1)=b(1)/2;for i=2:11u(i)=(b(i)-mu(i)*u(i-1)/(2-mu(i)*bata(i-1);endm(11

26、)=u(11);for i=10:-1:1m(i)=u(i)-bata(i)*m(i+1);end%給出需計(jì)算函數(shù)值的點(diǎn)z=-5:0.5:5;n=max(size(z);%計(jì)算樣條函數(shù)的值for j=1:ni=1;while i<=10if z(j)>=x(i) & z(j)<=x(i+1)s(j)=m(i)*(x(i+1)-z(j)3/(6*h(i)+m(i+1)*(z(j)-x(i)3/(6*h(i)+(y(i)-m(i)*h(i).2/6)*(x(i+1)-z(j)/h(i)+(y(i+1)-m(i+1)*h(i)2/6)*(z(j)-x(i)/h(i);bre

27、akelsei=i+1;endendend%計(jì)算拉氏插值多項(xiàng)式在z處的函數(shù)值for k=1:nL(k)=0;for i=1:11t(i)=1;for j=1:11if j=it(i)=t(i)*(z(k)-x(j)/(x(i)-x(j);endendL(k)=L(k)+t(i)*y(i);endend%計(jì)算被插函數(shù)在z處的值q=2*z./(1+z.2);%輸入值與圖象,作比較fprintf(' 節(jié)點(diǎn) 函數(shù)值 樣條值 拉氏插值n')disp(z',q',s',L')plot(z,L,'b')%藍(lán)線表示拉氏插值圖象plot(z,s,&

28、#39;y')%黃線表示三次插值樣條圖象hold off所得結(jié)果如圖6.2所示圖6.2 三次插值樣條與十次拉格朗日插值效果比較輸出結(jié)果如下： 節(jié)點(diǎn) 函數(shù)值 樣條值 拉氏插值 -5.0000 -0.3846 -0.3846 -0.3846 -4.5000 -0.4235 -0.3757 0.2572 -4.0000 -0.4706 -0.4706 -0.4706 -3.5000 -0.5283 -0.5455 -0.7007 -3.0000 -0.6000 -0.6000 -0.6000 -2.5000 -0.6897 -0.6682 -0.5970 -2.0000 -0.8000 -0.8000 -0.8000 -1.5000 -0.9231 -0.9905 -1.0195 -1.0000 -1.0000 -1.0000 -1.0000 -0.5000 -0.8000 -0.6198 -0.6264 0 0 0 0 0.5000