數(shù)值分析課程設(shè)計(jì)---多項(xiàng)式插值的振蕩現(xiàn)象

1、<p><b>　　數(shù)值分析課程設(shè)計(jì)</b></p><p>　　課程設(shè)計(jì)的目的和意義：</p><p>　　《課程設(shè)計(jì)》是數(shù)值分析的同步課程，是《數(shù)值分析》的上機(jī)實(shí)習(xí)課。</p><p>　　《數(shù)值分析》課程中構(gòu)造了各種有效的算法和有效公式，同學(xué)們通過上機(jī)作課程設(shè)計(jì)，學(xué)習(xí)揣摩這些算法的思想和構(gòu)造，評(píng)價(jià)算法的優(yōu)劣性。</p>

2、;<p>　　通過上機(jī)，可以提高我們運(yùn)用數(shù)學(xué)軟件編程解決問題的能力，為今后從事科學(xué)計(jì)算和軟件開發(fā)打下良好的基礎(chǔ)。</p><p><b>　　課程設(shè)計(jì)的題目：</b></p><p>　　多項(xiàng)式插值的振蕩現(xiàn)象</p><p><b>　　設(shè)計(jì)目的：</b></p><p>　　通過對(duì)多

3、項(xiàng)式插值現(xiàn)象的觀察,了解多項(xiàng)式的次數(shù)與逼近效</p><p>　　的關(guān)系，提高同學(xué)們分析實(shí)驗(yàn)結(jié)果的能力。</p><p><b>　　問題提出：</b></p><p>　　考慮在一個(gè)固定區(qū)間上用插值逼近一個(gè)函數(shù)。顯然，Lagrange插值中使用的節(jié)點(diǎn)越多，插值多項(xiàng)式的次數(shù)就越高。我們自然關(guān)心插值多項(xiàng)式增加時(shí)，Ln(x)是否也更加靠近被逼近的函

4、數(shù)。龍格(Runge)給出的一個(gè)例子是極著名并富有啟發(fā)性的。設(shè)區(qū)間[-1,1]上的函數(shù)</p><p>　　區(qū)間[-考慮設(shè)計(jì)1,1]的一個(gè)等距劃分，節(jié)點(diǎn)為</p><p>　　則拉個(gè)朗日插值多項(xiàng)式為</p><p>　　其中的li(x),i=0,1,2,…,n是n次Lagrange插值基函數(shù)。</p><p><b>　　設(shè)計(jì)要求：

5、</b></p><p>　　1．選擇不斷增大的分點(diǎn)數(shù)n=2,3,…</p><p>　　*畫出原函數(shù)f(x)及插值多項(xiàng)式函數(shù)Ln(x)在[-1,1]上的圖像；</p><p>　　*給出每一次逼近的最大誤差；</p><p>　　*比較并分析實(shí)驗(yàn)結(jié)果。</p><p>　　2．選擇其它函數(shù)，例如定義在區(qū)間

6、[-5,5]上的函數(shù)。</p><p>　　重復(fù)上述實(shí)驗(yàn)看其結(jié)果如何。</p><p>　　區(qū)間[a,b]上切比雪夫(Chebychev)點(diǎn)的定義為</p><p>　　以x1,x2,…,xn+1為插值節(jié)點(diǎn)構(gòu)造上述各函數(shù)的Lagrange插值多項(xiàng)式，比較其結(jié)果。</p><p><b>　　設(shè)計(jì)過程：</b></p

7、><p>　　已知函數(shù)f(x)在n+1個(gè)點(diǎn)x0,x1,…,xn處的函數(shù)值為 y0,y1,…,yn 。求一n次多項(xiàng)式函數(shù)Pn(x)，使其滿足：</p><p>　　Pn(xi)=yi,i=0,1,…,n.</p><p>　　解決此問題的拉格朗日插值多項(xiàng)式公式如下</p><p>　　其中Li(x) 為n次多項(xiàng)式：</p><p

8、>　?。?）.在MATLAB6.5中輸入函數(shù)</p><p>　　當(dāng)取不同的分點(diǎn)數(shù)n時(shí)，所得圖象與原函數(shù)圖象對(duì)比如下：</p><p>　　n=2時(shí)；最大誤差為： Max[L(x)－f(x)]=0.6462</p><p>　　n=3時(shí)，max[L(X)-f(X)]=0.7070；</p><p>　　n=6時(shí)，max[L(X)-f(X

9、)]=0.6169；</p><p>　　n=8時(shí)，max[L(X)-f(X)]=1.0452；</p><p>　　n=10時(shí)；最大誤差為：max[L(X)-f(X)]=1.9156；</p><p>　　n=20時(shí)，最大誤差為：max[L(X)-f(X)]=58.5855。</p><p><b>　　圖象分析：</b&g

10、t;</p><p>　　從圖中可以看出當(dāng)插值節(jié)點(diǎn)很少時(shí)，插值的誤差很大，插值圖象與原圖象沒有很好重疊在一起，而當(dāng)隨著插值的節(jié)點(diǎn)增加，中間能很好的重疊，但是兩邊出現(xiàn)很大誤差，隨著n值的增多，總體上分散的越厲害，最大誤差也逐漸增加，在n＝3時(shí)，最大誤差為0.6462，但到了n＝20時(shí)，已經(jīng)變成了58.5855，這種隨著節(jié)點(diǎn)數(shù)增多依然不能很好的接近被插值函數(shù)的現(xiàn)象稱為“龍格現(xiàn)象”，亦稱為多項(xiàng)式插值的振蕩現(xiàn)象。因此通過

11、增加節(jié)點(diǎn)數(shù)從而提高插值多項(xiàng)式的次數(shù)來逼近被插函數(shù)是不可取的。</p><p>　?。?）.對(duì)于定義在定義在區(qū)間[-5,5]上的函數(shù)h(x)=x/(1+x^4);</p><p>　　輸入函數(shù)后取節(jié)點(diǎn)數(shù)n=4，8，10；得到圖象如下</p><p>　　n=4時(shí)，max[L(X)-h(X)]=0.4020；</p><p>　　n=8時(shí)，max

12、[L(X)-h(X)]=0.1708；</p><p>　　n=10時(shí)，max[L(X)-h(X)]=0.1092。</p><p><b>　　圖象分析：</b></p><p>　　從圖象可以看出隨著插值節(jié)點(diǎn)數(shù)的增加出現(xiàn)異常的擺動(dòng)，中間能較好的接近原函數(shù)，但兩邊卻出現(xiàn)很大的誤差。</p><p>　　(3).對(duì)定義在

13、（-5，5）上的函數(shù)</p><p>　　g(x)=arctanx</p><p>　　在matlab里輸入相應(yīng)的lagrange插值公式，當(dāng)n取不同的值時(shí)，得到相應(yīng)圖象如下：</p><p>　　n=4時(shí)，max[L(X)-g(X)]=0.4458；</p><p>　　n=8時(shí)，max[L(X)-g(X)]=0.3240；</p&g

14、t;<p>　　N=10時(shí)，max[L(X)-g(X)]=0.8066。</p><p><b>　　圖象分析：</b></p><p>　　從圖象可以看出隨著插值節(jié)點(diǎn)數(shù)的增加出現(xiàn)異常的擺動(dòng)，中間能較好的接近原函數(shù)，隨著插值節(jié)點(diǎn)增加有較好的插值重疊，但兩邊卻出現(xiàn)很大的誤差。</p><p>　　（3）. 以上實(shí)驗(yàn)是從對(duì)區(qū)間的等距劃

15、分而作出的拉格朗日多項(xiàng)式，也因此產(chǎn)生了龍格振蕩現(xiàn)象，現(xiàn)在我們通過用切比雪夫點(diǎn)來對(duì)上述各函數(shù)的區(qū)間進(jìn)行劃分，以此來建立拉格朗日多項(xiàng)式，看這樣建立的拉格朗日多項(xiàng)式是否還會(huì)出現(xiàn)龍格振蕩現(xiàn)象。 還是考慮函數(shù) ，為了與等距節(jié)點(diǎn)進(jìn)行比較，我們?nèi)匀贿x取節(jié)點(diǎn)數(shù)不斷增多的拉格朗日多項(xiàng)式。得到不同圖象：</p><p><b>　　n=3時(shí)；</b></p><p

16、><b>　　n=6時(shí)；</b></p><p><b>　　n=11時(shí)；</b></p><p><b>　　n=21時(shí)。</b></p><p><b>　　圖象分析：</b></p><p>　　可見利用切比雪夫點(diǎn)來構(gòu)造的拉格朗日多項(xiàng)式比等距節(jié)

17、點(diǎn)下的拉格朗日多項(xiàng)式更為有效，隨著節(jié)點(diǎn)數(shù)的增加，逼近程度就越好，它能很好的消除了“龍格現(xiàn)象”。</p><p><b>　　心得體會(huì)：</b></p><p>　　深入了解matlab運(yùn)行環(huán)境和操作環(huán)境，初步學(xué)會(huì)調(diào)試程序，運(yùn)用繪圖命令制作函數(shù)圖象。</p><p>　　了解lagrange插值法龍格的異?，F(xiàn)象，以及數(shù)值分析的解決方案。</

18、p><p>　　懂得如何運(yùn)用已有的知識(shí)更進(jìn)一步了解未知的問題。</p><p>　　獨(dú)立解決和思考問題的能力有了一定的提高。</p><p><b>　　參考文獻(xiàn)：</b></p><p>　　《數(shù)值分析》華中科技大學(xué)出版社。李紅 著</p><p>　　《matlab從入門到精通》人民郵電出版社 求

19、是科技 編著</p><p>　　《數(shù)值分析》北京理工大學(xué)出版社 史萬明 孫新 等 編著</p><p><b>　　附錄</b></p><p><b>　　程序代碼1：</b></p><p>　　function t_charpt</p><p>　　result=in

20、putdlg({'請(qǐng)選擇實(shí)驗(yàn),若選2.1,請(qǐng)輸入1,否則輸入2:'},'charpt_2',1,{'1'});</p><p>　　Nb=str2num(char(result));</p><p>　　if(Nb~=1)&(Nb~=2) errordlg('實(shí)驗(yàn)選擇錯(cuò)誤!');return;end</p>

21、<p>　　promps={'請(qǐng)選擇實(shí)驗(yàn)函數(shù)，若選f(x),請(qǐng)輸入f,若選h(x),請(qǐng)輸入h,若選g(x),請(qǐng)輸入g:'};</p><p>　　% 選擇所要運(yùn)算的函數(shù)</p><p>　　titles='charpt_2';</p><p>　　result=inputdlg(promps,'charpt 2&

22、#39;,1,{'f'});</p><p>　　Nb_f=char(result);</p><p>　　if(Nb_f~='f'&Nb_f~='h'&Nb_f~='g') errordlg('實(shí)驗(yàn)選擇錯(cuò)誤！');</p><p>　　return;end</p&

23、gt;<p>　　result=inputdlg({'請(qǐng)輸入插值點(diǎn)數(shù)N:'},'charpt_2',1,{'10'});</p><p>　　Nd=str2num(char(result));</p><p>　　if(Nd<1)errordlg('結(jié)點(diǎn)輸入錯(cuò)誤！');return;end</p>

24、;<p>　　switch Nb_f</p><p><b>　　case 'f'</b></p><p>　　f=inline('1./(1+25*x.^2)');a=-1;b=1;</p><p><b>　　case 'h'</b></p>

25、<p>　　f=inline('x./(1+x.^4)');a=-5;b=5;</p><p><b>　　case 'g'</b></p><p>　　f=inline('atan(x)');a=-5;b=5;</p><p><b>　　end</b></

26、p><p><b>　　if(Nb==1)</b></p><p>　　x0=linspace(a,b,Nd+1);y0=feval(f,x0);</p><p>　　x=a:0.1:b;y=lagrange(x0,y0,x);</p><p>　　fplot(f,[a,b],'co');</p>

27、<p><b>　　hold on;</b></p><p>　　plot(x,y,'b--');</p><p>　　xlabel('x') ;ylabel('y=f(x)o and y=ln(x)--');</p><p><b>　　hold on;</b>

28、;</p><p>　　f1=1./(1+25*x.^2);</p><p>　　err=max(abs(y-f1));</p><p>　　result=inputdlg({'請(qǐng)輸入插值點(diǎn)數(shù)N:'},'charpt_2',1,{'10'});</p><p>　　xc=-cos(pi*[0:4

29、]/4);</p><p>　　x=(a+b)*0.5+(b-a)*xc*0.5;</p><p>　　elseif(Nb==2)</p><p>　　x0=linspace(a,b,Nd+1);y0=feval(f,x0);</p><p>　　x=a:0.1:b;</p><p>　　cs=spline(x0,y0

30、);y=ppval(cs,x);</p><p>　　plot(x0,y0,'o');hold on;plot(x,y,'k-');</p><p>　　xlabel('x');ylabel('y=f(x) o and y=spline(x)-');</p><p><b>　　end<

31、/b></p><p>　　function y=lagrange(x0,y0,x)</p><p>　　n=length(x0); m=length(x);</p><p><b>　　for i=1:m</b></p><p><b>　　z=x(i);</b></p>&l

32、t;p><b>　　s=0.0;</b></p><p><b>　　for k=1:n</b></p><p><b>　　p=1.0;</b></p><p><b>　　for j=1:n</b></p><p><b>　　if j

33、~=k</b></p><p>　　p=p*(z-x0(j))/(x0(k)-x0(j));</p><p><b>　　end</b></p><p><b>　　end</b></p><p>　　s=p*y0(k)+s;</p><p><b>　

34、　end</b></p><p><b>　　y(i)=s;</b></p><p><b>　　end</b></p><p>　　t1=clock %獲得系統(tǒng)時(shí)間</p><p><b>　　程序代碼2：</b></p><p>　　x

35、=[-1:0.2:1];</p><p>　　y=1./(1+25.*x.^2);</p><p>　　x0=[-1:0.01:1];</p><p>　　y0=lagrange(x,y,x0);</p><p>　　y1=1./(1+25.*x0.^2); </p><p>　　plot(x0,y0,'--

36、r');</p><p><b>　　hold on;</b></p><p>　　plot(x0,y1,'-b');</p><p>　　x2=abs(y0-y1); </p><p><b>　　max(x2) ;</b></p><p><

37、b>　　程序代碼3：</b></p><p><b>　　n=3; </b></p><p><b>　　for i=1:n</b></p><p>　　x(i)=cos(((2.*i-1).*pi)./(2.*(n+1)));</p><p>　　y(i)=1./(1+25.

38、*x(i).*x(i)); </p><p><b>　　end</b></p><p>　　x0=-1:0.01:1;</p><p>　　y0=lagrange(x,y,x0);</p><p>　　y1=1./(1+25.*x0.^2);</p><p>　　plot(x0,y0,'

39、--r')</p><p><b>　　hold on</b></p><p>　　plot(x0,y1,'-b')</p><p>　　計(jì)算機(jī)配置：AMD athlon（tm） 1.73GHz，512MB內(nèi)存。</p><p>　　CPU時(shí)間：（4.880000-0.020000，4.88000