數(shù)值分析課程設(shè)計(jì)--用拉格朗日插值法和牛頓插值法求近似值

1、<p>　　《數(shù)值分析課程設(shè)計(jì)》</p><p><b>　　報(bào) 告</b></p><p>　　【摘要】 本文簡介拉格朗日插值和牛頓插值。拉格朗日插值的算法及程序和拉格朗日在實(shí)際生活中的運(yùn)用。運(yùn)用了拉格朗日插值的公式，以及它在MATLAB中的算法程序，并用具體例子說明。拉格朗日插值在很多方面都可以運(yùn)用，具有很高的應(yīng)用價(jià)值。關(guān)于牛頓插值法，本文首先給出差商

2、的定義及性質(zhì)，由差商遞推得到Newton插值公式。在增加一個(gè)插值節(jié)點(diǎn)后，只需計(jì)算新增插值節(jié)點(diǎn)帶來的計(jì)算，而不必重新計(jì)算整個(gè)插值公式。然而并不是插值節(jié)點(diǎn)越多越好，插值多項(xiàng)式隨節(jié)點(diǎn)的增多而振動增多，反而不能更好的接近被插函數(shù)，這就是龍格現(xiàn)象。龍格現(xiàn)象從根本上否定了增多節(jié)點(diǎn)一提高插值多項(xiàng)式的次數(shù)來達(dá)到更好近似的可行性，從而產(chǎn)生了質(zhì)的飛躍。</p><p>　　【關(guān)鍵詞】 均差 ； 牛頓插值多項(xiàng)式 ； 龍格現(xiàn)象拉

3、格朗日；插值；公式；算法程序；應(yīng)用；科學(xué)。</p><p>　　題目：用拉格朗日插值法和牛頓插值法求近似值</p><p><b>　　二、理論</b></p><p>　　Lagrange插值法的理論：</p><p><b>　　1、基本概念</b></p><p>　　

4、已知函數(shù)y=f(x)在若干點(diǎn)的函數(shù)值=（i=0,1,,n）一個(gè)差值問題就是求一“簡單”的函數(shù)p(x)：p()=,i=0,1,,n, (1)</p><p>　　則p(x)為f(x)的插值函數(shù)，而f(x)為被插值函數(shù)會插值原函數(shù)，，，，...,為插值節(jié)點(diǎn)，式（1）為插值條件，如果對固定點(diǎn)求f()數(shù)值解，我們稱為一個(gè)插值節(jié)點(diǎn)，f()p()稱為點(diǎn)的插值，當(dāng)[min(，，，...,)，max(

5、，，，...,)]時(shí)，稱為內(nèi)插，否則稱為外插式外推，特別地，當(dāng)p(x)為不超過n次多項(xiàng)式時(shí)稱為n階Lagrange插值。</p><p>　　Lagrange插值公式</p><p><b>　　（1）線性插值</b></p><p>　　設(shè)已知 ， 及=f() ,=f(),為不超過一次多項(xiàng)式且滿足=,=,幾何上，為過（，），（，）的直線

6、，從而得到</p><p>　　=+（x-）. （2）</p><p>　　為了推廣到高階問題，我們將式（2）變成對稱式</p><p><b>　　=（x）+(x).</b></p><p><b>　　其中，</b></p><p>　?。▁）=,(x)=。均為

7、1次多項(xiàng)式且滿足</p><p>　?。▁）=1且(x)=0?；颍▁）=0且(x)=1。</p><p>　　兩關(guān)系式可統(tǒng)一寫成= 。 （3）</p><p>　　（2）n階Lagrange插值</p><p>　　設(shè)已知，，，...,及=f()(i=0,1,.....,n),為不超過n次多項(xiàng)式且滿足（i=0,1,...n）.</p&

8、gt;<p>　　易知=（x）+....+.</p><p>　　其中，均為n次多項(xiàng)式且滿足式（3）（i,j=0,1,...,n）,再由（ji）為n次多項(xiàng)式的n個(gè)根知=c.最后，由</p><p>　　c=,i=0,1,...,n.</p><p>　　總之，=，=式為n階Lagrange插值公式，其中，（i=0,1,...n）稱為n階Lagrange

9、插值的基函數(shù)。</p><p>　　3，Lagrange插值余項(xiàng)</p><p>　　設(shè)，，，...,[a,b],f(x)在[a,b]上有連續(xù)的n+1階導(dǎo)數(shù)，為f(x)關(guān)于節(jié)點(diǎn)，，，...,的n階Lagrange插值多項(xiàng)式，則對任意x[a,b],</p><p>　　其中，位于，，，...,及x之間（依賴于x），(x)=</p><p>　　

10、牛頓插值多項(xiàng)式的理論依據(jù)：</p><p><b>　　1. 均差定義</b></p><p>　　利用插值基函數(shù)很容易得到拉格朗日插值多項(xiàng)式，公式結(jié)構(gòu)緊湊，在理論分析中甚為方便，但當(dāng)插值節(jié)點(diǎn)增減時(shí)全部插值基函數(shù)均要隨之變化，整個(gè)公式也將發(fā)生變化，這在實(shí)際計(jì)算中是很不方便的，可把插值多項(xiàng)式表示為如下便于計(jì)算的形式</p><p><b&g

11、t;　　，</b></p><p>　　其中為待定系數(shù)，可又插值條件</p><p><b>　　確定．</b></p><p><b>　　當(dāng)時(shí)，．</b></p><p><b>　　當(dāng)時(shí)，，推得 </b></p><p>　　，

12、 ··· ···</p><p>　　依此遞推可得到．為寫出系數(shù)的一般表達(dá)式，先引進(jìn)如下均差定義．</p><p>　　定義1.1　記為f的零階均差，零階均差的差商記為　　　　稱為函數(shù)關(guān)于點(diǎn)的一階均差.一般地，記(k-1)階均差的差商為　　　　　　　　(1.1)稱為f關(guān)于點(diǎn)的k階均差.</p>&l

13、t;p>　　2． 均差的基本性質(zhì)</p><p>　?。保A均差可表為函數(shù)值的線性組合，即</p><p><b>　?。　　?1.2)</b></p><p>　　這個(gè)性質(zhì)可用歸納法證明．這性質(zhì)也表明均差于節(jié)點(diǎn)的排列次序無關(guān)，稱為均差的對稱性．即</p><p><b>　?。?lt;/b>&l

14、t;/p><p>　?。　　　　　　　?1.3)</p><p>　　，　　　　　　　　　 (1.4)</p><p>　?。玻尚再|(zhì)１及（1.1）可得</p><p>　　３．若在上存在階導(dǎo)數(shù)，且節(jié)點(diǎn)則階均差與導(dǎo)數(shù)關(guān)系如下：</p><p>　　這公式可直接用羅爾定理證明。</p><p>　　

15、均差計(jì)算可列均差表如下（表1-1）。</p><p><b>　　表1-1</b></p><p><b>　　3． 牛頓插值公式</b></p><p>　　根據(jù)均差定義，把看成上一點(diǎn)，可得</p><p><b>　　，</b></p><p>&l

16、t;b>　　，</b></p><p><b>　　……</b></p><p><b>　?。?lt;/b></p><p>　　只要把后一式代入前一式，就得到</p><p><b>　　其中</b></p><p><b>　

17、　(1.5)</b></p><p>　　插值余項(xiàng)為 ，　　　　?。?.6）</p><p>　　多項(xiàng)式顯然滿足插值條件，且次數(shù)不超過，．</p><p>　　我們稱為牛頓(Newton)均差插值多項(xiàng)式.系數(shù)就是均差表1-1中加橫線的各階均差。</p><p>　　三、方法、算法與程序設(shè)計(jì)</p><

18、p>　　Lagrange插值算法和程序：</p><p>　　function yy=nalagr(x,y,xx) </p><p>　　%用途：Lagrange插值法數(shù)值求解；格式：yy=nalagr(x,y,xx)</p><p>　　%x是節(jié)點(diǎn)向量，y是節(jié)點(diǎn)上的函數(shù)值，xx是插值點(diǎn)（可以多個(gè)），yy返回插值</p><p>　　m

19、=length(x);n=length(y);</p><p>　　if m~=n,error('向量x與y的長度必須一致');end</p><p><b>　　s=0;</b></p><p><b>　　for i=1:n</b></p><p>　　t=ones(1,leng

20、th(xx));</p><p>　　for j=1:n </p><p><b>　　if j~=i</b></p><p>　　t=t.*(xx-x(i))/(x(i)-x(j));</p><p><b>　　end</b></p><p><b>　　en

21、d</b></p><p>　　s=s+t*y(i);</p><p><b>　　end</b></p><p><b>　　yy=s;</b></p><p><b>　　牛頓插值：</b></p><p>　　1.newton插值多項(xiàng)式

22、的表達(dá)式如下：</p><p>　　其中每一項(xiàng)的系數(shù)ci的表達(dá)式如下：</p><p>　　即為f (x)在點(diǎn)處的i階差商，（，），由差商的性質(zhì)可知：</p><p>　　2.牛頓插值的程序?qū)崿F(xiàn)方法：</p><p><b>　　第一步：計(jì)算。</b></p><p>　　第二步：計(jì)算牛頓插值多項(xiàng)

23、式中，，得到n個(gè)多項(xiàng)式。</p><p>　　第三步：將第二步得到的n個(gè)多項(xiàng)式相加，得到牛頓插值多項(xiàng)式。</p><p>　　第四步：利用所得到的插值多項(xiàng)式，估算取其它值時(shí)的值。</p><p>　　第五步：作出所求多項(xiàng)式在插值結(jié)點(diǎn)周圍的函數(shù)圖像。</p><p>　　3.牛頓插值的程序設(shè)計(jì)：</p><p>　　fu

24、nction [p2,z]=newTon(x,y,t) </p><p>　　%輸入?yún)?shù)中x,y為元素個(gè)數(shù)相等的向量，t為待估計(jì)的點(diǎn)，可以為數(shù)字或向量。</p><p>　　%輸出參數(shù)中p2為所求得的牛頓插值多項(xiàng)式，z為利用多項(xiàng)式所得的t的函數(shù)值。</p><p>　　n=length(x);</p><p>　　chaS(1)=y(1);&

25、lt;/p><p><b>　　for i=2:n</b></p><p>　　x1=x;y1=y;</p><p>　　x1(i+1:n)=[];</p><p>　　y1(i+1:n)=[];</p><p>　　n1=length(x1);</p><p><b&g

26、t;　　s1=0;</b></p><p>　　for j=1:n1</p><p><b>　　t1=1;</b></p><p>　　for k=1:n1</p><p><b>　　if k==j</b></p><p><b>　　continu

27、e;</b></p><p><b>　　else</b></p><p>　　t1=t1*(x1(j)-x1(k));</p><p><b>　　end</b></p><p><b>　　end</b></p><p>　　s1=s1+

28、y1(j)/t1;</p><p><b>　　end</b></p><p>　　chaS(i)=s1;</p><p><b>　　end</b></p><p>　　b(1,:)=[zeros(1,n-1) chaS(1)];</p><p>　　cl=cell(1,n

29、-1);</p><p><b>　　for i=2:n</b></p><p><b>　　u1=1;</b></p><p>　　for j=1:i-1</p><p>　　u1=conv(u1,[1 -x(j)]);</p><p>　　cl{i-1}=u1;</

30、p><p><b>　　end</b></p><p>　　cl{i-1}=chaS(i)*cl{i-1};</p><p>　　b(i,:)=[zeros(1,n-i),cl{i-1}];</p><p><b>　　end</b></p><p>　　p2=b(1,:);&l

31、t;/p><p><b>　　for j=2:n</b></p><p>　　p2=p2+b(j,:);</p><p><b>　　end</b></p><p>　　if length(t)==1</p><p><b>　　rm=0;</b></

32、p><p><b>　　for i=1:n</b></p><p>　　rm=rm+p2(i)*t^(n-i);</p><p><b>　　end</b></p><p><b>　　z=rm;</b></p><p><b>　　else<

33、;/b></p><p>　　k1=length(t);</p><p>　　rm=zeros(1,k1);</p><p>　　for j=1:k1</p><p><b>　　for i=1:n</b></p><p>　　rm(j)=rm(j)+p2(i)*t(j)^(n-i);<

34、;/p><p><b>　　end</b></p><p><b>　　z=rm;</b></p><p><b>　　end</b></p><p><b>　　end</b></p><p>　　plot(t,z,'y

35、9;,x,y,'*r')</p><p><b>　　四、算例、應(yīng)用實(shí)例</b></p><p>　　例1已知函數(shù)表sin=0.5000,sin=0.7071,sin=0.8660,分別由線性插值與拋物插值求sin的數(shù)值解，并由余項(xiàng)公式估計(jì)計(jì)算結(jié)果的精度。</p><p>　　解：（1）這里有三個(gè)節(jié)點(diǎn)，線性插值需要兩個(gè)節(jié)點(diǎn)，根據(jù)

36、余項(xiàng)公式，我們選取前兩個(gè)節(jié)點(diǎn)，易知：</p><p>　　sin()=0.5000+(-)</p><p>　　=0.5000+0.2071=0.6381</p><p><b>　　截?cái)嗾`差，</b></p><p><b>　　=，</b></p><p>　　得知結(jié)果至

37、少有1位有效數(shù)字。</p><p><b>　　易知sin</b></p><p>　　0.7071+=0.8660=0.6434</p><p><b>　　截?cái)嗾`差為：</b></p><p>　　得知結(jié)果至少有兩位數(shù)字。</p><p>　　比較本題精確解sin=0.6

38、42787609...,實(shí)際誤差限分別為0.0047和0.00062。</p><p>　　用上述程序的1的結(jié)果為</p><p>　　>> x=pi*[1/6 1/4];y=[0.5 0.7071];xx=2*pi/9;</p><p>　　>> yy1=nalagr(x,y,xx)</p><p><b>

39、;　　yy1 =</b></p><p><b>　　-0.5690</b></p><p>　　>> x=pi*[1/6 1/4 1/3];y=[0.5 0.7071 0.866];</p><p>　　>> yy2=nalagr(x,y,xx)</p><p><b>　

40、　yy2 =</b></p><p><b>　　0.8023</b></p><p>　　>> fplot('sin',[pi/6,pi/3]);hold on;</p><p>　　>> plot(x,y,'o',xx,0.6381,'g^',xx,0.64

41、34,'rv');hold off;</p><p><b>　　圖形為</b></p><p>　　Lagrange插值應(yīng)用實(shí)例：</p><p>　　在物理化學(xué)，資產(chǎn)價(jià)值鑒定工作和計(jì)算某一時(shí)刻的衛(wèi)星坐標(biāo)和鐘差等這些方面可以應(yīng)用Lagrange插值。采用拉格朗日插值法計(jì)算設(shè)備等功能重置成本，計(jì)算精度較高，方法快捷。但是這方法只

42、能針對可比性較強(qiáng)的標(biāo)準(zhǔn)設(shè)備，方法本身也只考慮了單一功能參數(shù)，它的應(yīng)用范圍因此受到了一定的限制。作為一種探索，我們可以將此算法以及其它算法集成與計(jì)算機(jī)評估分析系統(tǒng)中，作為傳統(tǒng)評估分析方法的輔助參考工具，以提高資產(chǎn)價(jià)值鑒定工作的科學(xué)性和準(zhǔn)確性。</p><p>　　2. Newton基本插值公式</p><p>　　由均差可以導(dǎo)出n次Newton插值多項(xiàng)式</p><p&g

43、t;<b>　　其余項(xiàng)公式為</b></p><p>　　2.Newton基本插值公式的程序</p><p><b>　　2.1 程序Ⅰ</b></p><p>　　按上述公式編寫Matlab程序（文件名New_Int.m）.</p><p>　　function yi = New_Int(x,y,

44、xi)</p><p>　　% Newton基本插值公式，其中，</p><p>　　% x向量，全部的插值節(jié)點(diǎn)，按行輸入；</p><p>　　% y向量，插值節(jié)點(diǎn)處的函數(shù)值，按行輸入；</p><p>　　% xi向量，自變量x；</p><p>　　% yi為xi處的函數(shù)估計(jì)值.</p><p

45、>　　n = length(x);m=length(y);</p><p>　　% s輸入的插值點(diǎn)與他的函數(shù)值應(yīng)有相同的個(gè)數(shù).</p><p><b>　　if n ~= m</b></p><p>　　error('The lengths of X and Y must be equal!');</p>&

46、lt;p><b>　　return;</b></p><p><b>　　end</b></p><p><b>　　%計(jì)算均差表Y.</b></p><p>　　Y=zeros(n);Y(:,1) = y';</p><p>　　for k = 1:n-1<

47、;/p><p>　　for i = 1:n-k</p><p>　　% 輸入的插值節(jié)點(diǎn)必須互異.</p><p>　　if abs(x(i+k)-x(i))< eps</p><p>　　error('the DATA is error!');</p><p><b>　　return;&l

48、t;/b></p><p><b>　　end</b></p><p>　　Y(i,k+1)=(Y(i+1,k)-Y(i,k))/(x(i+k)-x(i));</p><p><b>　　end</b></p><p><b>　　end</b></p>&

49、lt;p>　　% 計(jì)算Newton插值公式N(xi)</p><p><b>　　yi=0;</b></p><p>　　for i = 1:n</p><p><b>　　z=1;</b></p><p>　　for k = 1:i-1</p><p>　　z = z

50、*(xi-x(k));</p><p><b>　　end</b></p><p>　　yi = yi+Y(1,i)* z;</p><p><b>　　End</b></p><p><b>　　Matlab演示:</b></p><p>　　例2：已

51、知=1，=2,=3,用Newton插值公式求的近似值。</p><p>　　在matlab軟件中，程序執(zhí)行如下：</p><p>　　>>x = [1 4 9]; y = [1 2 3]; xi = 5;</p><p>　　>> yi = New_Int(x,y,xi)</p><p><b>　　yi =

52、</b></p><p><b>　　2.2667</b></p><p><b>　　2.2 程序Ⅱ</b></p><p>　　建立一個(gè)M文件，命名為newpoly.m.Matlab命令程序如下：</p><p>　　function[c,d]=newpoly(x,y)</p&

53、gt;<p>　　% c為所求的Newton插值多項(xiàng)式的系數(shù)構(gòu)成的向量</p><p>　　n = length(x);</p><p>　　d = =zeros(n,n);</p><p>　　d(:,1)= y';</p><p>　　for j = 2:n</p><p>　　for k =

54、 j:n</p><p>　　d(k,j)= (d(k,j-1)-d(k-1,j-1))/(x(k)-x(k-j+1));</p><p><b>　　end</b></p><p><b>　　end</b></p><p>　　c = d(n,n);</p><p>　　

55、for k =(n-1):-1:1</p><p>　　c = conv(c,poly(x(k)));</p><p>　　m = length(c);</p><p>　　c(m)=c(m)+d(k,k);</p><p><b>　　end</b></p><p><b>　　Mat

56、lab演示：</b></p><p>　　例3：已知函數(shù)的觀測數(shù)據(jù)見下表，試用上述程序求牛頓插值多項(xiàng)式。</p><p><b>　　在命令窗口輸入</b></p><p>　　>> newpoly([1 2 3 4 5 6],[-3 0 15 48 105 192])</p><p><b

57、>　　ans =</b></p><p>　　0 0 1 0 -4 0</p><p>　　所以牛頓插值多項(xiàng)式為 </p><p>　　牛頓插值法的應(yīng)用實(shí)例：與拉格朗日插值法相比較更為方便，它不僅克服了“增加一個(gè)節(jié)點(diǎn)時(shí)整個(gè)計(jì)算工作必須重新開始”的缺點(diǎn)，而且可以節(jié)省乘除法運(yùn)算次數(shù)。同時(shí)，在牛頓插值多項(xiàng)式中用到的差分與差商等概念，與數(shù)值計(jì)

58、算的其他方面有著密切的關(guān)系。牛頓插值法來源于生產(chǎn)實(shí)踐，其應(yīng)用也日益廣泛。特別是在計(jì)算機(jī)的使用和航空、造船、精密機(jī)械加工，冶金工程等領(lǐng)域都獲得了極為廣泛的應(yīng)用。</p><p>　　例4（分別用拉格朗日插值法和牛頓插值法求近似值） 已知</p><p>　　分別用Lagrange插值法和牛頓插值法來求log0.78的近似值。</p><p>　　解：用拉格朗日方法在m

59、atlab中編寫程序，得到結(jié)果如下：</p><p>　　>> clear all;</p><p>　　>> x=[0.4,0.5,0.6,0.7,0.8,0.9];</p><p>　　>> y=[log(0.4),log(0.5),log(0.6),log(0.7),log(0.8),log(0.9)];</p>

60、<p>　　>> x0=[0.78];</p><p>　　>> disp('差值點(diǎn)')</p><p><b>　　差值點(diǎn)</b></p><p>　　>> du=[0.78]</p><p><b>　　du =</b><

61、/p><p><b>　　0.7800</b></p><p>　　>> disp('差值結(jié)果')</p><p><b>　　差值結(jié)果</b></p><p>　　>> yt=Lagrange(x,y,x0)</p><p><

62、b>　　yt =</b></p><p><b>　　-0.2485</b></p><p>　　>> disp('log函數(shù)值')</p><p><b>　　log函數(shù)值</b></p><p>　　>> yreal=[log(0.78)

63、]</p><p><b>　　yreal =</b></p><p><b>　　-0.2485</b></p><p>　　>> disp('差值與函數(shù)值誤差')</p><p><b>　　差值與函數(shù)值誤差</b></p><

64、;p>　　>> dy=yt-yreal</p><p><b>　　dy =</b></p><p>　　8.1720e-006</p><p>　　用牛頓插值法在matlab中編寫程序，結(jié)果如下：</p><p>　　>> clear all;</p><p>　

65、　>> x=[[0.4,0.5,0.6,0.7,0.8,0.9];];</p><p>　　>> y=[-0.916291-0.693147-0.510826-0.357762-0.223144-0.105361];</p><p>　　>> x0=0.78;</p><p>　　>> disp('用

66、牛頓插值法在x=1.4處插值：')</p><p>　　用牛頓插值法在x=1.4處插值：</p><p>　　>> yx=Newtoninter(x,y,x0)</p><p><b>　　yx =</b></p><p><b>　　-0.7555</b></p>

67、<p><b>　　圖形如下：</b></p><p><b>　　五、附錄</b></p><p>　　本文簡介拉格朗日插值和牛頓插值。拉格朗日插值的算法及程序和拉格朗日在實(shí)際生活中的運(yùn)用。運(yùn)用了拉格朗日插值的公式，以及它在MATLAB中的算法程序，并用具體例子說明。拉格朗日插值在很多方面都可以運(yùn)用，具有很高的應(yīng)用價(jià)值。關(guān)于牛頓插值

68、法，本文首先給出差商的定義及性質(zhì)，由差商遞推得到Newton插值公式。在增加一個(gè)插值節(jié)點(diǎn)后，只需計(jì)算新增插值節(jié)點(diǎn)帶來的計(jì)算，而不必重新計(jì)算整個(gè)插值公式。然而并不是插值節(jié)點(diǎn)越多越好，插值多項(xiàng)式隨節(jié)點(diǎn)的增多而振動增多，反而不能更好的接近被插函數(shù)，這就是龍格現(xiàn)象。龍格現(xiàn)象從根本上否定了增多節(jié)點(diǎn)一提高插值多項(xiàng)式的次數(shù)來達(dá)到更好近似的可行性，從而產(chǎn)生了質(zhì)的飛躍。</p><p><b>　　六、參考文獻(xiàn)</

69、b></p><p>　　1，約瑟夫·拉格朗日</p><p>　　2，作者，張玲。文章名，拉格朗日插值在資產(chǎn)評估中的應(yīng)用。</p><p>　　3，作者，宮厚誠，李全海。文章名，基于IGS精密星歷的衛(wèi)星坐標(biāo)和鐘差插值。</p><p>　　4，作者，吳法倫，趙占芬。文章名，利用計(jì)算機(jī)繪制物理化學(xué)實(shí)驗(yàn)中的曲線——拉格朗日插值&