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

2、法、分段低次插值法及三次樣條插值法，給出了每種插值算法的算法結構、具體的程序MATLAB代碼及誤差估計的方法，并且對每種插值方法的優缺點進行了分析。關鍵詞：拉格朗日插值，牛頓插值，埃爾米特插值，分段低次插值，三次樣條插值問題提出 設f (x)的函數表為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)。請根據給出的函數及函數表解下面的問題。 (1)證明插值多項式的存在唯一性，以x ,y為節點計算其插值多項式。（2）用基函數構造拉格朗日插值多項式，以x ,y為節點做拉格朗日插值曲線。（3）分析牛

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

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

5、一地確定一個n個拉氏插值多項式。下面我們在問題1的條件下計算以x=2,4,6,8,10,y=9.8,9.2,8.1,6.4,3.8為插值節點的插值多項式。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運行結果為：B =10.1600,-0.1476,-0.0071,-0.0042' 這表明: 。2 用基函數構造拉格朗日插值多項式2.1 拉格朗日插值多項式的構造利用求解線性方程組(1.2)求拉格朗日插值多項式不便計算機實現，不易編程。在實際應用中，一般采用基函數

6、法來構造。下面，我們介紹這一方法，這里，我們將n次拉格朗日插值多項式用來記。設,其中是待定的n次多項式。首先構造個插值節點上的插值基函數，對任一點所對應的插值基函數，由于在所有取零值，因此有因子。又因是一個次數不超過的多項式，所以只可能相差一個常數因子，固可表示成：利用得：于是因此滿足 的插值多項式可表示為：從而次拉格朗日插值多項式為：2.2拉格朗日插值多項式在MATLAB中的實現及使用在MATLAB中構造拉格朗日插值函數代碼如下：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;利用該函數做題一中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牛頓多項式

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

9、差對其節點是對稱的，即節點按任意順序排列，其值不變。如這個性質稱為對稱性。性質3 如果是的次多項式，則其1階均差是的次多項式，且由此遞推可得階均差為零階均差，階均差為零。差商表的構造如表1.3所示：表3.1 差商表 零階均差1階均差2階均差 3階均差 利用MATLAB計算題1中x，y的各階均差代碼如下：x=2,4,6,8,10'y=9.8,9.2,8.1,6.4,3.8'n=max(size(x); %計算節點個數z=zeros(n,n+1);%建一個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 階差商endfor i=j:ny(i)=z(i,j+2); %將j 階差商值賦給向量y 的對應分量endendz計算結果為：3.2牛頓插值多項式的算法及余項由均差的定義可知對上述第二式兩端乘以，第三式兩端乘以, 最后一式兩端乘以，然后由后一式的左端代入前一式的右端(第二項)，即可得 顯然，是次多項式，又有，故有，從而滿足條件，這就說明是關于的次插值多項式，通常稱它為Newton插值公式，為插值余項。在MATLAB中使用牛頓插值法求得題1中以x，y為插值節點的條件下f(6.596)的值

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

12、f('截斷誤差R(%8.5f)=%15.14fn',z,r);計算結果如下：4 埃爾米特插值多項式 拉氏插值多項式的還存在著一個缺陷：在插值節點處，被插函數曲線與插值多項式曲線的切線方向可能會不一致。不少實際的插值問題不但要求在節點上函數值相等，而且還要求對應的導數值也相等，甚至要求高階導數也相等，滿足這種要求的插值多項式稱為埃爾米特（Hermite）插值多項式。4.1埃爾米特插值多項式定義設在節點上，若存在一個2n+1次多項式H(x)，使得 成立，則稱H(x)為f (x)的2n+1次埃爾米特插值多項式。4.2埃爾米特插值多項式的構造及余項設計2n+2個插值基函數，每一個基函數

13、都是2n+1次多項式，并且滿足以下條件于是滿足條件(2.5)的插值多項式H(x)可以表示成為下面的問題就是求滿足條件(2.6)的基函數。a a 利用拉氏插值基函數以及條件(2.6)，可知的形式應為又整理得解出由于于是類似地，可推導出。仿照拉格朗日余項的證明方法，若f (x)在(a,b)內的2n+2階導數存在，則其插值余項為在MATLAB中構造埃爾米特函數代碼如下： 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的導數的維數不相等！'); return; end else disp('y和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_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中計算以x，y為插值節點，且對應的一階導數為y的埃爾米特多項式，并計算在x=6.596點處函數的值代碼如下：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)計算結果為：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龍格現象 lagrange插值方法使根據區間a,b上給出的節點構造插值多項式近似f(x)。一般認為的次數n越高，逼近f(x)的效果越好，龍格現象表明 時，不一定收斂到f(x)。我們以函數為例研究龍格現象，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('龍格現象 ')圖5.1 龍格現象這說明高次插值多項式近似的效果并不好，因而通常并不采用高次插值，而采用分段低次插值。5.2分段線性插值分段線性插值就是通過插值點用折線段連接起來逼近。設已知節點a=x0 <x1 <<xn =b上的函數值yk =f(xk)(k=0,1,n)，已知函數y= f(x) 在這些插值結點的函數值為yk =f(xk)(k=0,1,n)求一個分段函數Ih(x), 使其滿足：  (1) Ih(xk )=yk ，(k=0,1

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

20、)n=input('n = ');%作被插函數圖象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);%取需用分段線性插值函數計算其函數值的點zz=-5:10/(2*n):5;m=max(size(z);%計算分段線性插值函數在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%釋放固化的圖形屏幕輸入等分數為5時，拉格朗日插值曲線如圖5.2所示，圖5.2 等分數為5時的拉格朗日插值曲線輸入等分數為800時，拉格朗日插值曲線如圖5.3所示， 圖5.3 等分數為800時的拉格朗日插值曲線6 三次樣條插值多項式通過平面上在n+1個不同的已知點()，可以作一條光滑的n次代數曲線，但當點很多時，除了插值多項式的次數高，計算復雜外，從逼近的效果來看，隨著插值多項式的次數增高，逼近的效果

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

23、有如果令那么解就可以為 6. 2三次樣條插值多項式算例分析做以x，y為插值節點的三次樣條插值曲線代碼如下：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所得結果如圖6.1所示：圖6.1 三次樣條插值曲線6.3三次樣條插值函數與拉格朗日插值多項式的效果比較為顯示三次樣條插值函數與拉格朗日插值多項式的效果比較，我們用MATLAB分別做g（x）=，(-5x5)的函數曲線、拉格朗日插值曲線、

24、三次樣條插值曲線，代碼如下：clearformat short%作y=2*x/(1+x2)的圖象u=-5:(10/200):5;v=2*u./(1+u.2);plot(u,v,'r')%紅線表示函數圖象hold on%給出插值條件x=-5:1:5;y=2*x./(1+x.2);y_bar(1)=152/(26)2; %在-5,5的一階導數值y_bar(2)=152/(26)2;%計算每個小區間長度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);%用追趕法解三對角方程,求矩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%給出需計算函數值的點z=-5:0.5:5;n=max(size(z);%計算樣條函數的值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%計算拉氏插值多項式在z處的函數值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%計算被插函數在z處的值q=2*z./(1+z.2);%輸入值與圖象,作比較fprintf(' 節點 函數值 樣條值 拉氏插值n')disp(z',q',s',L')plot(z,L,'b')%藍線表示拉氏插值圖象plot(z,s,&

28、#39;y')%黃線表示三次插值樣條圖象hold off所得結果如圖6.2所示圖6.2 三次插值樣條與十次拉格朗日插值效果比較輸出結果如下： 節點 函數值 樣條值 拉氏插值 -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