數(shù)值分析實(shí)習(xí)報(bào)告插值法.doc

1、x大學(xué)計(jì)算機(jī)與信息學(xué)院（數(shù)學(xué)類課程）課程實(shí)習(xí)報(bào)告課程名稱：數(shù)值分析課程實(shí)習(xí)實(shí)習(xí)題目：插值法姓 名：系：信息與計(jì)算科學(xué)專 業(yè)：信息與計(jì)算科學(xué)年 級：2009學(xué) 號：指導(dǎo)教師：職 稱：講師2012 年 05 月 1 日課程實(shí)習(xí)報(bào)告評分表序號項(xiàng)目評分權(quán)值評定成績1算法：正確，步驟合理0.32程序：畫流程圖，用語言或數(shù)學(xué)軟件編寫0.23質(zhì)量：分析處理科學(xué)；文字通順；計(jì)算及測試結(jié)果準(zhǔn)確；0.44工作量、工作態(tài)度：按期完成規(guī)定的任務(wù)，工作量飽滿；工作努力，遵守紀(jì)律0.1合 計(jì)1指導(dǎo)教師簽名： 年 月 日目錄1.實(shí)習(xí)的目的和任務(wù)12.實(shí)習(xí)要求13.實(shí)習(xí)地點(diǎn)14.主要儀器設(shè)備15.實(shí)習(xí)內(nèi)容15.1問題描述15

2、.2lagrange插值法15.3逐次線性插值法65.4牛頓插值法85.5等距節(jié)點(diǎn)插值法105.6 hermite插值法155.7分段低次插值法186問題討論與分析197結(jié)束語20參考文獻(xiàn)20數(shù)值分析課程實(shí)習(xí)1. 實(shí)習(xí)的目的和任務(wù)目的：通過課程實(shí)習(xí)能夠應(yīng)用matlab軟來計(jì)算函數(shù)的插值；了解函數(shù)插值方法。任務(wù)：通過具體的問題，利用matlab軟件來計(jì)算問題的結(jié)果，分析問題的結(jié)論。2. 實(shí)習(xí)要求能夠從案例的自然語言描述中，抽象出其中的數(shù)學(xué)模型；能夠熟練應(yīng)用所學(xué)的函數(shù)插值方法；能夠熟練使用matlab軟件；對函數(shù)的插值方法有所認(rèn)識，包括對不同插值法有所認(rèn)識。3. 實(shí)習(xí)地點(diǎn)學(xué)生宿舍4. 主要儀器設(shè)備

3、計(jì)算機(jī)microsoft windows 7matlab r2010b5. 實(shí)習(xí)內(nèi)容5.1 問題描述給出函數(shù)f(x)=lnx的數(shù)值表（見表5.1）用lagrange插值法，逐次線性插值法，牛頓插值法，等距節(jié)點(diǎn)插值法，hermite插值法，分段低次插值法分別求該函數(shù)0.4，0.8上的插值多項(xiàng)式并計(jì)算ln0.53的近似值。 表5.1x0.40.50.60.70.8lnx-0.916291-0.693147-0.510826-0.357765-0.2231445.2 lagrange插值法拉格朗日插值法是基于基函數(shù)的插值方法，插值多項(xiàng)式可表示為： 其中稱為i次基函數(shù)5.2.1線性插值和拋物插值線性插

4、值選取x=0.4 x=0.5 在matlab命令窗口中輸入：x=0.4 0.5;y=-0.916291 -0.693147;l0=poly(x(2)/(x(1)-x(2);l1=poly(x(1)/(x(2)-x(1);z=y(1)*l0+y(2)*l1;a=poly2sym(z) 回車結(jié)果為 a=(27893*x)/12500-1.8088670b=polyval(z,0.53) 回車結(jié)果為 b=-0.6262038在matlab命令窗口中輸入：x1=0.3:0.001:0.6;y1=(27893*x1)/12500-1.8088670;x2=0.4 0.5 0.53;y2=-0.91629

5、1 -0.693147 -0.6262038;plot(x1,y1,b-,x2,y2,o);title(圖 5-2-1 線性插值)回車后得到如圖所示結(jié)果拋物插值選取x=0.4 x=0.5 x=0.6在matlab命令窗口中輸入：x=0.4 0.5 0.6;y=-0.916291 -0.693147 -0.510826;l0=poly(x(2) x(3)/(x(1)-x(2)*(x(1)-x(3);l1=poly(x(1) x(3)/(x(2)-x(1)*(x(2)-x(3);l2=poly(x(1) x(2)/(x(3)-x(1)*(x(3)-x(2);z=y(1)*l0+y(2)*l1+y(

6、3)*l2;a=poly2sym(z)回車結(jié)果為 a=-(40823*x2)/20000+4.068475*x-2.217097b=polyval(z,0.53) 回車結(jié)果為 b=-0.6341643在matlab命令窗口中輸入：x=0.35:0.001:0.65;y=-(40823*x.2)/20000+4.068475*x-2.217097;x1=0.4 0.5 0.53 0.6;y1=-0.916291 -0.693147 -0.6262038 -0.510826;plot(x,y,g+,x1,y1,o);title(圖 5-2-2 拋物插值)回車后得到如圖所示結(jié)果5.2.2拉格朗日插值

7、法的matlab實(shí)現(xiàn)建立m文件language.m:function f=language(x,y,x0)%求已知數(shù)據(jù)點(diǎn)的拉格朗日插值多項(xiàng)式%已知數(shù)據(jù)點(diǎn)的x坐標(biāo)向量: x%已知數(shù)據(jù)點(diǎn)的y坐標(biāo)向量: y%插值點(diǎn)的x坐標(biāo): x0%求得的拉格朗日插值多項(xiàng)式或在x0處的插值: f syms t;if(length(x)=length(y) n=length(x);else disp(x和y的維度不相等！); return;%檢錯endf=0.0;for(i=1:n) l=y(i); for(j=1:i-1) l=l*(t-x(j)/(x(i)-x(j); end for(j=i+1:n) l=l*(t

8、-x(j)/(x(i)-x(j); end; f=f+l; simplify(f);%化簡 if(i=n) if(nargin=3) f=subs(f,t,x0);%計(jì)算插值點(diǎn)的函數(shù)值 else f=collect(f); f=vpa(f,8);%將插值多項(xiàng)式的系數(shù)化成8位精度的小數(shù) end endend在matlab工作窗口輸入如下程序：x=0.4 0.5 0.6 0.7 0.8;y=-0.916291 -0.693147 -0.510826 -0.357765 -0.223144;f=lagrange(x,y)運(yùn)行的結(jié)果為：f =- 0.30958333*t4 + 2.60825*t3 -

9、 5.4860542*t2 + 5.6920925*t - 2.474362再次輸入f=lagrange(x,y,0.53)結(jié)果為 f= -0.6347048在matlab命令窗口中輸入：t=0.35:0.001:0.85;y=- 0.30958333*t.4+ 2.60825*t.3 - 5.4860542*t.2 + 5.6920925*t - 2.474362;x1=0.4 0.5 0.53 0.6 0.7 0.8;y1=-0.916291 -0.693147 -0.6347048 -0.510826 -0.357765 -0.223144;plot(t,y,r*,x1,y1,o);ti

10、tle(圖 5-2-3 lagrange插值)回車后得到如圖所示結(jié)果5.3逐次線性插值法5.3.1逐次線性插值法基本原理定義為，一階定義均差為 一般地，函數(shù)f的k階均差定義為：5.3.2逐次線性插值法的matlab實(shí)現(xiàn)建立m文件aitken.mfunction f = aitken(x,y,x0)%求已知數(shù)據(jù)點(diǎn)的aitken逐次線性插值多項(xiàng)式%已知數(shù)據(jù)點(diǎn)的x 坐標(biāo)向量:x%已知數(shù)據(jù)點(diǎn)的y 坐標(biāo)向量:y%為插值點(diǎn)的x坐標(biāo):x0%求得的aitken插值多項(xiàng)式或x0處的插值:fsyms t;if(length(x) = length(y) n = length(x);else disp(x和y的維數(shù)

11、不相等！); return; %檢錯endy1(1:n)=t; for(i=1:n-1) for(j=i+1:n) y1(j) = y(j)*(t-x(i)/(x(j)-x(i)+y(i)*(t-x(j)/(x(i)-x(j); end y=y1; simplify(y1); %化簡endif(nargin = 3) f = subs(y1(n),t,x0); %計(jì)算插值點(diǎn)的函數(shù)值else simplify(y1(n); f = collect(y1(n); %將插值多項(xiàng)式展開 f = vpa(f,8); %將插值多項(xiàng)式的系數(shù)化成8位精度的小數(shù)end在matlab工作窗口輸入：x=0.4 0.

12、5 0.6 0.7 0.8;y=-0.916291 -0.693147 -0.510826 -0.357765 -0.223144;f= aitken(x,y)運(yùn)行的結(jié)果為：f =- 0.30958333*t4 + 2.60825*t3 - 5.4860542*t2 + 5.6920925*t - 2.474362再次輸入f= aitken(x,y,0.53)結(jié)果為f= -0.6347048在matlab命令窗口中輸入：t=0.35:0.001:0.85;y=- 0.30958333*t.4 + 2.60825*t.3 - 5.4860542*t.2 + 5.6920925*t - 2.474

13、362x1=0.4 0.5 0.53 0.6 0.7 0.8;y1=-0.916291 -0.693147 -0.6347048 -0.510826 -0.357765 -0.223144;plot(t,y,cd,x1,y1,o);title(圖 5-3 aitken插值)回車后得到如圖所示結(jié)果5.4牛頓插值法5.4.1牛頓插值法的基本原理零階均差定義為，一階定義均差為 一般地，函數(shù)f的k階均差定義為：利用均差的牛頓插值法多項(xiàng)式為5.4.2牛頓插值法的matlab實(shí)現(xiàn)建立m文件newton.m：function f = newton(x,y,x0)%求已知數(shù)據(jù)點(diǎn)的均差形式牛頓插值多項(xiàng)式%已知數(shù)

14、據(jù)點(diǎn)的x 坐標(biāo)向量:x%已知數(shù)據(jù)點(diǎn)的y 坐標(biāo)向量:y%為插值點(diǎn)的x坐標(biāo):x0%求得的均差形式牛頓插值多項(xiàng)式或x0處的插值:fsyms t;if(length(x) = length(y) n = length(x); c(1:n) = 0.0;else disp(x和y的維數(shù)不相等！); return; %檢錯endf = y(1);y1 = 0;l = 1; for(i=1:n-1) for(j=i+1:n) y1(j) = (y(j)-y(i)/(x(j)-x(i); end c(i) = y1(i+1); l = l*(t-x(i); f = f + c(i)*l; simplify(f

15、); %化簡 y = y1; if(i=n-1) if(nargin = 3) f = subs(f,t,x0); %計(jì)算插值點(diǎn)的函數(shù)值 else f = collect(f); %將插值多項(xiàng)式展開 f = vpa(f,8); %將插值多項(xiàng)式的系數(shù)化成8位精度的小數(shù) end endend在matlab工作窗口輸入如下程序：x=0.4 0.5 0.6 0.7 0.8;y=-0.916291 -0.693147 -0.510826 -0.357765 -0.223144;f= newton (x,y)運(yùn)行的結(jié)果為f = - 0.30958333*t4 + 2.60825*t3 - 5.486054

16、2*t2 + 5.6920925*t - 2.474362再次輸入 f= newton (x,y,0.53)結(jié)果為 f= -0.6348783在matlab命令窗口中輸入：t=0.35:0.001:0.85;y=- 0.30958333*t.4 + 2.60825*t.3 - 5.4860542*t.2 + 5.6920925*t - 2.474362x1=0.4 0.5 0.53 0.6 0.7 0.8;y1=-0.916291 -0.693147 -0.6348783 -0.510826 -0.357765 -0.223144;plot(t,y,m- ,x1,y1,o);title(圖 5

17、-4 newton插值)回車后得到如圖所示結(jié)果5.5等距節(jié)點(diǎn)插值法5.5.1向前牛頓插值法的基本原理向前牛頓插值多項(xiàng)式可表示如下：其中h叫做步長，且t的取值范圍為 。5.5.2向前牛頓插值法的matlab實(shí)現(xiàn)建立m文件newtonforward.m：function f = newtonforward(x,y,x0)%求已知數(shù)據(jù)點(diǎn)的向前差分牛頓插值多項(xiàng)式%已知數(shù)據(jù)點(diǎn)的x 坐標(biāo)向量:x%已知數(shù)據(jù)點(diǎn)的y 坐標(biāo)向量:y%為插值點(diǎn)的x坐標(biāo):x0%求得的向前差分牛頓插值多項(xiàng)式或x0處的插值:fsyms t;if(length(x) = length(y) n = length(x); c(1:n) =

18、0.0;else disp(x和y的維數(shù)不相等！); return; %檢錯endf = y(1);y1 = 0;xx =linspace(x(1),x(n),(x(2)-x(1);if(xx = x) disp(節(jié)點(diǎn)之間不是等距的！); return;endfor(i=1:n-1) for(j=1:n-i) y1(j) = y(j+1)-y(j); end c(i) = y1(1); l = t; for(k=1:i-1) l = l*(t-k); end; f = f + c(i)*l/factorial(i); simplify(f); %化簡 y = y1; if(i=n-1) if(

19、nargin = 3) f = subs(f,t,(x0-x(1)/(x(2)-x(1); else f = collect(f); %將插值多項(xiàng)式展開 f = vpa(f,8); %將插值多項(xiàng)式的系數(shù)化成8位精度的小數(shù) end endend在matlab工作窗口輸入如下程序：x=0.4 0.5 0.6 0.7 0.8;y=-0.916291 -0.693147 -0.510826 -0.357765 -0.223144;f= newtonforward (x,y)運(yùn)行的結(jié)果為f = - 0.000030958333*t4 + 0.0021129167*t3 - 0.026533542*t2

20、+ 0.24759558*t - 0.916291再次輸入 f= newtonforward (x,y,0.53)結(jié)果為 f= -0.6348783在matlab命令窗口中輸入：x=0.4 0.5 0.6 0.7 0.8;y=-0.916291 -0.693147 -0.510826 -0.357765 -0.223144;f= newtonforward(x,y)t=0.4:0.01:0.8;yy1=newtonforward(x,y,t);plot(x,y,ro,t,yy1, : )title(圖 5-5-1 向前牛頓插值法)回車后得到如圖所示結(jié)果5.5.3向后牛頓插值法的基本原理向后牛頓

21、插值多項(xiàng)式可表示如下：其中h叫做步長，且t的取值范圍為。5.5.4向后牛頓插值法的matlab實(shí)現(xiàn)建立m文件newtonback.m：function f = newtonback(x,y,x0)%求已知數(shù)據(jù)點(diǎn)的向后差分牛頓插值多項(xiàng)式%已知數(shù)據(jù)點(diǎn)的x 坐標(biāo)向量:x%已知數(shù)據(jù)點(diǎn)的y 坐標(biāo)向量:y%為插值點(diǎn)的x坐標(biāo):x0%求得的向前差分牛頓插值多項(xiàng)式或x0處的插值:fsyms t;if(length(x) = length(y) n = length(x); c(1:n) = 0.0;else disp(x和y的維數(shù)不相等！); return;endf = y(n);y1 = 0;xx =lins

22、pace(x(1),x(n),(x(2)-x(1);if(xx = x) disp(節(jié)點(diǎn)之間不是等距的！); return; %檢錯endfor(i=1:n-1) for(j=i+1:n) y1(j) = y(j)-y(j-1); end c(i) = y1(n); l = t; for(k=1:i-1) l = l*(t+k); end; f = f + c(i)*l/factorial(i); simplify(f); %化簡 y = y1;if(i=n-1) if(nargin = 3) f = subs(f,t,(x0-x(n)/(x(2)-x(1); else f = collect

23、(f); %將插值多項(xiàng)式展開 f = vpa(f,8); %將插值多項(xiàng)式的系數(shù)化成8位精度的小數(shù) end endend在matlab工作窗口輸入：x=0.4 0.5 0.6 0.7 0.8;y=-0.916291 -0.693147 -0.510826 -0.357765 -0.223144;f= newtonback (x,y)運(yùn)行的結(jié)果為f = - 0.000030958333*t4 + 0.0016175833*t3 - 0.0041505417*t2 + 0.12882192*t - 0.223144再次輸入 f= newtonback (x,y,0.53)結(jié)果為 f= -0.6347

24、048在matlab命令窗口中輸入：x=0.4 0.5 0.6 0.7 0.8;y=-0.916291 -0.693147 -0.510826 -0.357765 -0.223144;f= newtonback (x,y)t=0.4:0.01:0.8;yy1=newtonback (x,y,t);plot(x,y,mo,t,yy1, :x )title(圖 5-5-2 向后牛頓插值法)回車后得到如圖所示結(jié)果5.6 hermite插值法5.6.1 hermite插值法的基本原理hermite插值滿足在節(jié)點(diǎn)上等于給定函數(shù)值，而且在節(jié)點(diǎn)上的導(dǎo)數(shù)值也等于給定導(dǎo)數(shù)值。對于高階導(dǎo)數(shù)情況，hermite插值

25、多項(xiàng)式比較復(fù)雜，在實(shí)際中，常遇到的是函數(shù)值與一階導(dǎo)數(shù)給定的情況。在此情況下，n個(gè)節(jié)點(diǎn)x1，x2，xn的hermite插值多項(xiàng)式的表達(dá)式如下：其中，。5.6.2 hermite插值法的matlab實(shí)現(xiàn)建立m文件hermite.m：function f = hermite(x,y,y_1,x0)%已知數(shù)據(jù)點(diǎn)的x 坐標(biāo)向量:x%已知數(shù)據(jù)點(diǎn)的y 坐標(biāo)向量:y%已知數(shù)據(jù)點(diǎn)的導(dǎo)數(shù)向量:y_1%求得的hermite插值多項(xiàng)式或x0處的插值:fsyms t;f = 0.0;if(length(x) = length(y) if(length(y) = length(y_1) n = length(x); el

26、se disp(y和y的導(dǎo)數(shù)的維數(shù)不相等！); return; endelse disp(x和y的維數(shù)不相等！); return; endfor 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,8); %將插值多項(xiàng)式的系數(shù)化成8位精度的小數(shù) e

27、nd endend在matlab工作窗口輸入：x=0.4 0.5 0.6 0.7 0.8;y=-0.916291 -0.693147 -0.510826 -0.357765 -0.223144;y_1=5/2 2 5/3 10/7 5/4;f= hermite (x,y, y_1)運(yùn)行的結(jié)果為f = 6250000.0*(1.6666667*t - 1.510826)*(t - 0.4)2*(t - 0.8)2*(t - 0.7)2*(t - 0.5)2 + 173611.11*(10.547667*t - 8.6612773)*(t - 0.4)2*(t - 0.7)2*(t - 0.6)2

28、*(t - 0.5)2 + 2777777.8*(7.3913214*t - 5.53169)*(t - 0.4)2*(t - 0.8)2*(t - 0.6)2*(t - 0.5)2 - 173611.11*(35.678792*t - 13.355226)*(t - 0.8)2*(t - 0.7)2*(t - 0.6)2*(t - 0.5)2 - 2777777.8*(9.55245*t - 4.083078)*(t - 0.4)2*(t - 0.8)2*(t - 0.7)2*(t - 0.6)2再次輸入 f= hermite (x,y, y_1,0.53)結(jié)果為 f= -0.6349413

29、在matlab命令窗口中輸入：x=0.4 0.5 0.6 0.7 0.8;y=-0.916291 -0.693147 -0.510826 -0.357765 -0.223144;y_1=5/2 2 5/3 10/7 5/4;f= hermite (x,y,y_1)t=0.4:0.01:0.8;yy1=hermite(x,y,y_1,t);plot(x,y,o,t,yy1)title(圖 5-6 hermite插值)回車后得到如圖所示結(jié)果5.7分段低次插值法5.7.1分段線性插值在matlab命令窗口中輸入：x=0.4 0.5 0.6 0.7 0.8;y=-0.916291 -0.693147 -0.510826 -0.357765 -0.223144;xi=0.53;yi=interp1(x,y,xi,linear)結(jié)果為 yi =-0.63845070在matlab命令窗口中輸入：x=0.4 0.5 0.6 0.7 0.8;y=-0.916291 -0.693147 -0.510826 -0.357765 -0.223144;xi=0.4:0.1:0.8;yi=interp1(x,y,xi,linear);plot(x,y,ro,x,yi);title(圖 5-7 分