1、精選優質文檔-傾情為你奉上分 數：評卷人：研究生（數值分析）實驗報告學 號 M姓 名 劉良專 業 光學工程指 導 教 師 劉小英院（系、所） 光學與電子信息學院2018年4月30日數值實驗2.1 多項式差值的振蕩現象問題提出：考慮在一個固定的區間上用插值逼近一個函數，顯然Lagrange插值中使用的節點越多，插值多項式的次數就越高，我們自然關心插值多項式的次數增加時，是否也更加靠近被逼近的函數。Runge給出的一個例子是極其著名并富有啟發性的，設區間-1,1上函數實驗內容：考慮區間-1,1的一個等距劃分，分點為則拉格朗日插值多項式為其中是n次Lagrange插值基函數。實驗要求：（1） 選擇不

2、斷增大的分點數目n=2,3,···，畫出原函數f(x)插值多項式函數(x)在-1,1上的圖像，比較并分析實驗結果。（2） 選擇其他的函數，例如定義在區間-5,5上的函數重復上述的實驗看其結果如何。1、 MATLAB實驗程序prompt='如選擇f(x),請輸入f,如選擇h(x),請輸入h,如選擇g(x),請輸入g:'word='請選擇試驗函數'answer=inputdlg(prompt,word,1,'f');choose=char(answer);prompt='請輸入插值多項式的次數：'word

3、='插值次數'answer1=inputdlg(prompt,word,1,'10');numb=str2num(char(answer1);if(numb<1) errordlg('請輸入大于0的自然數！') ruturn;endswitch choose case'f' f=(x) 1./(1+25*x.2);a=-1;b=1; case'h' f=(x) x./(1+x.4);a=-5;b=5; case'g' f=(x) atan(x);a=-5;b=5;endx0=linspace(

4、a,b,numb+1);y0=f(x0);x=a:0.1:b;y=lagrange(x0,y0,x);clf;fplot(f,a b,'r-');hold on;plot(x,y,'b-');xlabel('x');ylabel('精確函數（紅色），插值函數（藍色）');str1=strcat('插值次數為',answer1);str2=strcat('試驗函數為',answer)title(str1,str2);%lagrange插值函數function y=lagrange(x0,y0,x)n=

5、length(x0);m=length(x); for k=1:m z=x(k); s=0.0; for i= 1:n p=1.0; for j=1:n if (j=i) p=p*(z-x0(j)/(x0(i)-x0(j); end end s=s+p*y0(i); end y(k)=s; endend2、 實驗結果對f(x)，g(x)，h(x)，依次輸入分點數目為n=2,3,5,8,10,13,15,20,25,其原函數與其插值多項式函數Ln的圖像如下圖所示：當實驗函數選擇為f(x)時，取不同分點數目的原函數和插值多項式函數的圖像：當實驗函數選擇為h(x)時，取不同分點數目的原函數和插值多項

6、式函數的圖像： 當實驗函數選擇為g(x)時，取不同分點數目的原函數和插值多項式函數的圖像： 3、 結果分析（1） 對于f(x)的實驗結果分析可知： 首先觀察到，無論插值次數選擇為奇數還是偶數，其插值函數都為偶函數，這是因為原函數是偶函數，在計算插值函數L(x)的過程中，其奇數次的項全為0。隨著分點數目的增加，插值函數的次數和項數也在增加，其復雜性也在增加，同時計算次數也在增加，由于計算機位數限制而引入的計算誤差也在增加。其次，觀察插值次數由高到低的插值函數擬合情況，當插值次數較低（n=2,3,5）時，拉格朗日插值函數不能很好的與原函數曲線吻合，隨著次數的增加插值函數越來越逼近與原函數，次數越高

7、插值函數越吻合原函數，但是這只局限于插值區間的中間部分。后面可以觀察到插值次數增加的時候，在靠近插值區間端點會出現極大的波動。當n=13時，其偏離程度為0.5左右，n=15時達到0.7，n=20時擬合情況最好，n=25時在端點的偏離程度達到4.5。因此，在拉格朗日基函數的高次插值中，在插值區間的邊界部分插值函數會出現很大波動，明顯偏離原函數，所以插值次數不宜過高，容易出現龍格現象，但是也不能選擇的過低，不然拉格朗日插值函數在區間中部的擬合情況不好。（2） 對于h(x)、g(x)的實驗結果分析可知: h(x)和g(x)均為奇函數，插值函數也為奇函數。與f(x)的實驗分析結果類似，在低次插值時，插

8、值函數不能很好的擬合，在選擇過高的插值次數時，同樣出現龍格振蕩現象。在等距節點的拉格朗日插值中，插值次數不易過高，在高次插值中龍格現象無法避免，而且超過一定次數后，龍格現象會越來越嚴重。為解決龍格現象可以采用分段低次插值或采用切比雪夫零點為插值點。數值實驗3 函數逼近與曲線擬合（1）編制以函數為基的多項式最小二乘擬合程序，并用與對表3.11中的數據作3次多項式最小二乘擬合。 -1.0-0.50.00.51.01.52.0 -4.447-0.4520.5510.048-0.4470.5494.552取權數，求擬合曲線中的參數、平方誤差，并做離散數據的擬合函數的圖形。（3） 編制正交化多項式最小二

9、乘擬合程序，并用于求解上題中的3次多項式最小二乘擬合問題，作擬合曲線的圖形，計算平方誤差，并與上題結果進行比較。1、 MATLAB實驗程序answer=inputdlg('請輸入1選擇實驗3.1，否則輸入2','實驗選擇窗口',1,'1');nb=str2num(char(answer);x0=-1.0:0.5:2.0;y0=-4.447 -0.452 0.551 0.048 -0.447 0.549 4.552;n=3;%擬合次數w= 1 1 1 1 1 1 1 ;if (nb=1) alph=wlsf(x0,y0,w,n);%按照降次排列的系

10、數矩陣 y=polyval(alph,x0);%擬合多項式在x0處的值 r=(y0-y)*(y0-y)'%平方誤差 x=-1:0.01:2; y=polyval(alph,x); plot(x,y,'r-'); xlabel('x');ylabel('擬合曲線（紅色）離散數據（*）'); hold on plot(x0,y0,'*'); title('離散數據的多項式擬合'); grid on; else a,b,c,alph,r=oplsf(x0,y0,w,n);enddisp('平方誤差：

11、9;,sprintf('%g',r);disp('系數為(按降次)：',sprintf('%gt',alph)function p=wlsf(x0,y0,w,n)%帶權重的最小二乘擬合子函數format shortA=zeros(n+1,n+1);for i=0:n for j=0:n A(i+1,j+1)=sum(x0.(i+j)*w'); end b(i+1)=sum(x0.i.*y0)*w');enda=Ab'p=fliplr(a');end % g1:初始值為0矩陣；g2:初始值設置為一個長度1*m的單位一

12、矩陣；g3:根據遞推關系%v1:初始值：g(k-1)的內積；v2:g(k)的內積；v3:作為一個中間變量。%d(k):y與的內積；c(k)=存儲的是正交多項式的系數。function a,b,c,alph,r=oplsf(x,y,w,n)%正交化多項式最小二乘擬合子函數 m=length(x);g1=0; g2=ones(1,m);%g2表示默認（1,1,1,1,1···，m+1） v2=g2.2*w' d(1)=y*w'c(1)=d(1)/v2; for k=1:n xs=x.*g2.2*w'a(k)=xs/v2;%表示gk（s2）的內積

13、與x的乘積，v2表示g2的內積 if (k=1) b(k)=0; else b(k)=v2/v1; end g3=(x-a(k).*g2-b(k).*g1; v3=g3.2*w' d(k+1)=y.*g3*w' c(k+1)=d(k+1)/v3; g1=g2;g2=g3;v1=v2;v2=v3;%改變遞推值 end r=y.*y*w'-c*d'%平方誤差 alph=zeros(1,n+1);t=zeros(n+1,n+2); t(:,2)=ones(n+1,1);t(2,3)=-a(1); if (n>=2) for k=3:n+1 for i=3:k+1

14、 t(k,i)=t(k-1,i)-a(k-1)*t(k-1,i-1)-b(k-1)*t(k-2,i-2); end end end for i=1:n+1 for k=i:n+1 alph(n+2-i)=alph(n+2-i)+c(k)*t(k,k+2-i); end end xmin=min(x);xmax=max(x);dx=(xmax-xmin)/(25*m); t=(xmin-dx):dx:(xmax+dx); s=alph(1); for k=2:n+1 s=s.*t+alph(k); end plot(x,y,'*',t,s,'-'); title(

15、'正交多項式擬合'); xlabel('x');ylabel('y'); grid on ;end2、 實驗結果對表中的數據作3次最小二乘擬合，其離散數據的擬合函數的圖形如圖3.1所示，輸出的平方誤差為：2.17619e-05，系數為(按降次)：1.99911-2.99767-3.96825e-050.。對表中的數據通過正交多項式3次最小二乘擬合，其離散數據的擬合函數的圖形如圖3.2所示，輸出平方誤差：2.17619e-05。系數為(按降次)：1.99911-2.99767-3.96825e-050.當擬合次數n取不同的值時，其平方誤差得到如下表

16、格：n23456平方誤差13.4882.17619e-051.6671e-051.63734e-051.2374e-253、 結果分析最小二乘擬合是在離散情況下的最佳平方逼近。從圖中我們可以看到，在擬合次數為3次的情況下，各個離散點均能落在擬合曲線上，擬合效果很好，擬合的平方誤差為10-5量級。利用多項式進行最小二乘擬合和利用正交多項式進行最小二乘擬合的結果對比發現，兩種算法算出的擬合函數系數一致，平方誤差一致，所得的結果完全一致，實際上這也是必然的。在算法原理上，用一般多項式進行最小二乘擬合時，其系數矩陣是一個正規方程組，而利用正交多項式的解法時，無非是為了避免多項式擬合時求解病態法方程組，

17、通過尋找一組正交的基函數，使得法方程的系數矩陣變成一個對角陣，使得法方程組得到簡化，利于求解。本質上還是同一種算法，所得的結果必然是一致的。由于題中給出了7個離散點，通過依次改變擬合次數n，觀察計算得到平方誤差的表，在擬合次數越高時（小于離散點數目），其平方誤差越小，擬合效果也越好。數值實驗4.1 數值積分實驗目的：復化積分公式計算定積分。實驗題目：數值計算下列各式右端定積分的近似值。 實驗要求：（1） 若用復化梯形公式、復化Simpson公式和復化Gauss-legendre I型公式做計算，要求絕對誤差限為，分別利用他們的余項對每種算法做出步長的事前估計。（2） 分別用復化梯形公式、復化S

18、impson公式和復化Gauss-legendre I型公式做計算。（3） 將計算結果與精確解作比較，并比較各種算法的計算量。1、 MATLAB實驗程序answer=inputdlg('請選擇積分公式，復化梯形公式輸入T，復化simpson公式輸入S，高斯公式輸入G','積分公式選擇窗口',1,'T');answer1=inputdlg('請選擇積分題目14','題目選擇窗口',1,'1');nb=char(answer);nb1=str2num(char(answer1);r=0.5e-7;swi

19、tch nb1 case 1 f=(x) -2./(x.2-1);a=2;b=3; case 2 f=(x) 4./(1+x.2);a=0;b=1; case 3 f=(x) 3.x;a=0;b=1; case 4 f=(x) x.*exp(x);a=1;b=2;endswitch nb case 'T' tic; t=(f(a)+f(b)*(b-a)/2; %分成一等分的復化梯形公式計算值。 k=1;t0=0; while(abs(t-t0)>=r*3) t0=t;k=k+1; h=(b-a)/k; t=h/2*(f(a)+f(b)+2*sum(f(a+h:h:b-h)

20、;%梯形公式的一般計算式 end time=toc; case 'S' tic; t=(f(a)+f(b)+4*f(a+b)/2)*(b-a)/6; k=1;t0=0; while (abs(t-t0)>=r*3)%事前估計誤差大約為絕對誤差的3倍 t0=t;k=k+1;h=(b-a)/k; t=h/6*(f(a)+4*sum(f(a+h/2:h:b-h/2)+2*sum(f(a+h:h:b-h)+f(b);%simpson公式的一般計算式 end time=toc; case 'G' tic;h=b-a; t=h/2*(f(a+h/2-h/(2*sqrt

21、(3)+f(a+h/2+h/(2*sqrt(3); k=1;t0=0; while(abs(t-t0)>=r*3) t0=t;k=k+1;h=(b-a)/k; t=h/2*sum(f(a+h/2-h/(2*sqrt(3):h:b-h/2-h/(2*sqrt(3)+f(a+h/2+h/(2*sqrt(3):h:b-h/2+h/(2*sqrt(3);%高斯公式的一般計算式 end time=toc;end switch nb1 case 1 disp('精確解，In2-In3=-0.'); disp('絕對誤差：',num2str(abs(t+0.); dis

22、p('運行時間：',num2str(time); disp('所選擇的k等分',num2str(k); case 2 disp('精確解，pi=3.979'); disp('絕對誤差：',num2str(abs(t-3.979); disp('運行時間：',num2str(time); disp('所選擇的k等分',num2str(k); case 3 disp('精確解，2/In3=1.368'); disp('絕對誤差：',num2str(abs(t-1.368)

23、; disp('運行時間：',num2str(time); disp('所選擇的k等分',num2str(k); case 4 disp('精確解，e2=7.065'); disp('絕對誤差：',num2str(abs(t-7.065); disp('運行時間：',num2str(time); disp('所選擇的k等分',num2str(k); end2、 實驗結果(1)題目一所示的函數：采用復化梯形公式得到的結果如下：精確解，In2-In3=-0.絕對誤差：6.7578e-06運行時間：0.所選擇的k等分93采用復化simpson公式得到的結果如下：精確解，In2-In3=-0.絕對誤差：1.9236e-07運行時間：0.所選擇的k等分10采用復化Gauss-legendre I型公式得到的結果如下：精確解，In2-In3=-0.絕對誤差：1.9505e-07運行時間：0.所選擇的k等分9(2) 題目二所示的函數：采用復化梯形公式得到的結果如下：精確解，pi=3.979絕對誤差：9.7119e-06運行時間：0.所選擇的k等分131采用復化simpson公式得到的結果如下：精確解，pi=3.979絕對誤差：3.9651e-08運行時間：0.所選擇的k等分5