數值分析實驗報告

1、 實驗2.1 多項式插值的振蕩現象實驗目的：在一個固定的區(qū)間上用插值逼近一個函數，顯然Lagrange插值中使用的節(jié)點越多，插值多項式的次數就越高。我們自然關心插值多項式的次數增加時，Ln(x)是否也更加靠近被逼近的函數。Runge給出的一個例子是極著名并富有啟發(fā)性的。實驗容：設區(qū)間-1，1上函數 f(x)=1/(1+25x2)?？紤]區(qū)間-1,1的一個等距劃分，分點為 xi= -1 + 2i/n，i=0，1，2，n，則拉格朗日插值多項式為.其中，li(x)，i=0，1，2，n是n次Lagrange插值基函數。實驗步驟與結果分析：實驗源程序function Chap2Interpolation%

2、 數值實驗二：“實驗2.1：多項式插值的震蕩現象”% 輸入：函數式選擇，插值結點數% 輸出：擬合函數及原函數的圖形promps = '請選擇實驗函數，若選f(x),請輸入f,若選h(x),請輸入h,若選g(x),請輸入g:'titles = 'charpt_2'result = inputdlg(promps,'charpt 2',1,'f');Nb_f = char(result);if(Nb_f = 'f' & Nb_f = 'h' & Nb_f = 'g')er

3、rordlg('實驗函數選擇錯誤！');return;endresult = inputdlg('請輸入插值結點數N:','charpt_2',1,'10');Nd = str2num(char(result);if(Nd <1)errordlg('結點輸入錯誤！');return;endswitch Nb_f case 'f' f=inline('1./(1+25*x.2)'); a = -1;b = 1; case 'h' f=inline('x./

4、(1+x.4)'); a = -5; b = 5; case 'g' f=inline('atan(x)'); a = -5; b= 5;end x0 = linspace(a, b, Nd+1); y0 = feval(f, x0); x = a:0.1:b; y = Lagrange(x0, y0, x); fplot(f, a b, 'co'); hold on; plot(x, y, 'b-'); xlabel('x'); ylabel('y = f(x) o and y = Ln(x)-&

5、#39;);%-function y=Lagrange(x0, y0, x);n= length(x0); m=length(x);for i=1:m z=x(i); s=0.0; for k=1:n p=1.0; for j=1:n if(j = k) p = p*(z - x0(j)/(x0(k) - x0(j); end end s = s + p*y0(k); end y(i) = s;end實驗結果分析(1)增大分點n=2，3，時，拉格朗日插值函數曲線如圖所示。 n=6 n=7 n=8 n=9 n=10從圖中可以看出，隨著n的增大，拉格朗日插值函數在x=0附近較好地逼近了原來的函數f

6、(x)，但是卻在兩端x= -1和x=1處出現了很大的振蕩現象。并且，仔細分析圖形，可以看出，當n為奇數時，雖然有振蕩，但振蕩的幅度不算太大，n為偶數時，其振蕩幅度變得很大。通過思考分析，我認為，可能的原因是f(x)本身是偶函數，如果n為奇數，那么Lagrange插值函數Ln(x)的最高次項xn-1是偶次冪，比較符合f(x)本身是偶函數的性質；如果n為偶數，那么Lagrange插值函數Ln(x)的最高次項xn-1是奇次冪，與f(x)本身是偶函數的性質相反，因此振蕩可能更劇烈。(2)將原來的f(x)換為其他函數如h(x)、g(x)，結果如圖所示。其中h(x), g(x)均定義在-5，5區(qū)間上，h(

7、x)=x/(1+x4)，g(x)=arctan x。h(x), n=7 h(x), n=8h(x), n=9 h(x), n=10g(x), n=7g(x), n=8g(x), n=9g(x), n=10分析兩個函數的插值圖形，可以看出：隨著n的增大，拉格朗日插值函數在x=0附近較好地逼近了原來的函數f(x)，但是卻在兩端x= -5和x=5處出現了很大的振蕩現象。并且，仔細分析圖形，可以看出，當n為偶數時，雖然有振蕩，但振蕩的幅度不算太大，n為奇數時，其振蕩幅度變得很大。原因和上面f(x)的插值類似，h(x)、g(x)本身是奇函數，如果n為偶數，那么Lagrange插值函數Ln(x)的最高次項

8、xn-1是奇次冪，比較符合h(x)、g(x)本身是奇函數的性質；如果n為奇數，那么Lagrange插值函數Ln(x)的最高次項xn-1是偶次冪，與h(x)、g(x)本身是奇函數的性質相反，因此振蕩可能更劇烈。實驗3.1 多項式最小二乘擬合實驗目的： 編制以函數xkk=0,n；為基的多項式最小二乘擬合程序。實驗容： 對表中的數據作三次多項式最小二乘擬合。xi-1.0-0.50.00.51.01.52.0yi-4.447-0.4520.5510.048-0.4470.5494.552取權函數wi1，求擬合曲線中的參數k、平方誤差2，并作離散據xi, yi的擬合函數的圖形。實驗源程序function

9、 Chap3CurveFitting% 數值實驗三:“實驗3.1”% 輸出:原函數及求得的相應插值多項式的函數的圖像以及參數alph和誤差rx0 = -1:0.5:2;y0 = -4.447 -0.452 0.551 0.048 -0.447 0.549 4.552;n = 3; % n為擬合階次alph = polyfit(x0, y0, n);y = polyval(alph, x0);r = (y0 -y)*(y0 -y)' % 平方誤差x = -1:0.01:2;y = polyval(alph, x);plot(x, y, 'k-');xlabel('

10、x'); ylabel('y0 * and polyfit.y-');hold onplot(x0, y0, '*')grid on;disp('平方誤差:', num2str(r)disp('參數alph:', num2str(alph)實驗結果平方誤差:2.1762e-005參數alph:1.9991 -2.9977 -3.9683e-005 0.54912實驗4.1 實驗目的：復化求積公式計算定積分. 實驗題目：數值計算下列各式右端定積分的近似值. 實驗要求： （1）若用復化梯形公式、復化Simpson公式和復化Ga

11、uss-Legendre I 型公式做計算，要求絕對誤差限為，分別利用它們的余項對每種算法做出步長的事前估計. （2）分別用復化梯形公式，復化Simpson 公式和復化Gauss-Legendre I 型公式作計算. （3）將計算結果與精確解做比較，并比較各種算法的計算量. 實驗程序：1.事前估計的Matlab程序如下： （1）用復化梯形公式進行事前估計的Matlab程序 format long g x=2:0.01:3; f=-4*(3*x.2+1)./(x.2-1).3; %二階導函數 %plot(x,f) %畫出二階導函數圖像 x=2.0; %計算導函數最大值 f=-4*(3*x2+1)

12、/(x2-1)3; h2=0.5*10(-7)*12/f; h=sqrt(abs(h2) %步長 n=1/h; n=ceil(1/h)+1 %選取的點數 format long g x=0:0.01:1; f=8.*(3*x.2-1)./(x.2+1).3; %二階導函數 %plot(x,f) %畫出二階導函數圖像 x=1; %計算導函數最大值 f=8.*(3*x.2-1)./(x.2+1).3; h2=0.5*10(-7)*12/f; h=sqrt(abs(h2) %步長 n=1/h n=ceil(1/h)+1 %選取的點數 format long g x=0:0.01:1; f=log(3

13、).*log(3).*3.x; %二階導函數 %plot(x,f); %畫出二階導函數圖像 x=1; %計算導函數最大值 f=log(3)*log(3)*3x; h2=0.5*10(-7)*12/f; h=sqrt(abs(h2) %步長 n=1/h n=ceil(1/h)+1 %選取的點數 format long g x=1:0.01:2; f=2.*exp(x)+x.*exp(x);%二階導函數 %plot(x,f) %畫出二階導函數圖像 x=2; %計算導函數最大值 f=2.*exp(x)+x.*exp(x); h2=0.5*10(-7)*12/f; h=sqrt(abs(h2) %步長

14、 n=1/h n=ceil(1/h)+1 %選取的點數估計結果步長h及結點數n分別為 h = 0.8 n =1793 h =0.6 n =1827 h =0.9 n =2458 h =0.9 n =7020 （2）用復化simpson公式進行事前估計的Matlab程序 format long g x=2:0.01:3; f=-2*(-72*x.2-24).*(x.2-1)-192*x.2.*(x.2+1)./(x.2-1).5;%四階導函數 x=2.0; f=-2*(-72*x2-24)*(x2-1)-192*x2*(x2+1)/(x2-1)5; %計算導函數最大值 h4=0.5*10(-7)

15、*180*16/f; h=sqrt(sqrt(abs(h4) %步長 n=1/h; %求分段區(qū)間個數 n=2*ceil(1/h)+1 %選取的點數 format long g x=0:0.01:1; f=4*(-72*x.2+24).*(x.2+1)-192*x.2.*(-x.2+1)./(x.2+1).5;%四階導函數 x=1; f=4*(-72*x2+24)*(x2+1)-192*x2*(-x2+1)/(x2+1)5; %計算導函數最大值 h4=0.5*10(-7)*180*16/f; h=sqrt(sqrt(abs(h4)%步長 n=1/h; %求分段區(qū)間個數 n=2*ceil(1/h)

16、+1 %選取的點數 format long g x=0:0.01:1; f=log(3)4*3.x;%四階導函數 x=1; f=log(3)4*3.x;%計算導函數最大值 h4=0.5*10(-7)*180*16/f; h=sqrt(sqrt(abs(h4)%步長 n=1/h; %求分段區(qū)間個數 n=2*ceil(1/h)+1 %選取的點數 format long g x=1:0.01:2; f=4*exp(x)+x.*exp(x);%四階導函數 plot(x,f) %畫出原函數 x=2; f=4*exp(x)+x.*exp(x); %計算導函數最大值 h4=0.5*10(-7)*180*16

17、/f; h=sqrt(sqrt(abs(h4) n=1/h; %求分段區(qū)間個數 n=2*ceil(1/h)+1 %選取的點數 估計結果步長h及結點數n分別為 h =0.13411 n =47 h =0.76542 n =35 h =0.18433 n =29 h =0.18546 n =49 2. 積分計算的Matlab程序： format long g promps='請選擇積分公式，若用復化梯形，請輸入T，用復化simpson，輸入S，用復化Gauss_Legendre，輸入GL：' result=inputdlg(promps,'charpt 4',1,&

18、#39;T'); Nb=char(result); if(Nb='T'&Nb='S'&Nb='GL') errordlg('積分公式選擇錯誤'); return; end result=inputdlg('請輸入積分式題號1-4：','實驗4.1',1,'1'); Nb_f=str2num(char(result); if(Nb_f<1|Nb_f>4) errordlg('沒有該積分式'); return; end switch N

19、b_f case 1 fun=inline('-2./(x.2-1)');a=2;b=3; case 2 fun=inline('4./(x.2+1)');a=0;b=1; case 3 fun=inline('3.x');a=0;b=1; case 4 fun=inline('x.*exp(x)');a=1;b=2; end if(Nb='T')%用復化梯形公式 promps='請輸入用復化梯形公式應取的步長：' result=inputdlg(promps,'實驗4.2',1,&

20、#39;0.01'); h=str2num(char(result); if(h<=0) errordlg('請輸入正確的步長！'); return; end tic; N=floor(b-a)/h); detsum=0; for i=1:N-1 xk=a+i*h; detsum=detsum+fun(xk); end t=h*(fun(a)+fun(b)+2*detsum)/2; time=toc; t end if(Nb='S')%用復化Simpson公式 promps='請輸入用復化Simpson公式應取的步長：' resul

21、t=inputdlg(promps,'實驗4.2',1,'0.01'); h=str2num(char(result); if(h<=0) errordlg('請輸入正確的步長！'); return; end tic; N=floor(b-a)/h); detsum_1=0; detsum_2=0; for i=1:N-1 xk_1=a+i*h; detsum_1=detsum_1+fun(xk_1); end for i=1:N xk_2=a+h*(2*i-1)/2; detsum_2=detsum_2+fun(xk_2); end t=

22、h*(fun(a)+fun(b)+2*detsum_1+4*detsum_2)/6; time=toc; t end if(Nb='GL')%用復化Gauss_Legendre I %先根據復化Gauss_Legendre I公式的余項估計步長 promps='請輸入用復化Gauss_Legendre I 公式應取的步長：' result=inputdlg(promps,'實驗4.2',1,'0.01'); h=str2num(char(result); if(h<=0) errordlg('請輸入正確的步長!

23、9;); return; end tic; N=floor(b-a)/h);t=0; for k=0:N-1 xk=a+k*h+h/2; t=t+fun(xk-h/(2*sqrt(3)+fun(xk+h/(2*sqrt(3); end t=t*h/2; time=toc; t end switch Nb_f case 1 disp('精確解：ln2-ln3=-0.4054651081') disp('絕對誤差：',num2str(abs(t+0.4054651081); disp('運行時間：',num2str(time); case 2 dis

24、p('精確解：pi=3.979') disp('絕對誤差：',num2str(abs(t-pi); disp('運行時間：',num2str(time); case 3 disp('精確解：2/ln3=1.368') disp('絕對誤差：',num2str(abs(t-1.368); disp('運行時間：',num2str(time); case 4 disp('精確解：e2=7.065') disp('絕對誤差：',num2str(abs(t-7.065); d

25、isp('運行時間：',num2str(time); end1. 當選用復化梯形公式時： （1）式運行結果為： t =-0.351 精確解：ln2-ln3=-0.4054651081 絕對誤差：1.3944e-008 運行時間：0.003 （2）式運行結果為： t =3.336 精確解：pi=3.979 絕對誤差：3.9736e-008 運行時間：0.005 （3）式運行結果為： t = 1.861 精確解：2/ln3=1.368 絕對誤差：4.3655e-008 運行時間：0.016 （4）式運行結果為： t =7.610 精確解：e2=7.065 絕對誤差：2.0775e-

26、008 運行時間：0.007 2. 當選用復化Simpson公式進行計算時： （1）式運行結果為： t =-0.7519 精確解：ln2-ln3=-0.4054651081 絕對誤差：2.7519e-011 運行時間：0.022 （2）式運行結果為： t =3.979 精確解：pi=3.979 絕對誤差：0 運行時間：0.021 （3）式運行結果為： t =1.288 精確解：2/ln3=1.368 絕對誤差：9.2018e-012 運行時間：0.019 （4）式運行結果為： t =7.118 精確解：e2=7.065 絕對誤差：9.0528e-011 運行時間：0.021 3. 當選用復化G

27、auss-Legendre I型公式進行計算時： （1）式運行結果為： t =-0.5262 精確解：ln2-ln3=-0.4054651081 絕對誤差：4.7385e-012 運行時間：0.023 （2）式運行結果為： t =3.979 精確解：pi=3.979 絕對誤差：1.3323e-015 運行時間：0.021 （3）式運行結果為： t =1.754 精確解：2/ln3=1.368 絕對誤差：6.1431e-012 運行時間：0.019 （4）式運行結果為： t =7.046 精確解：e2=7.065 絕對誤差：1.441e-012 運行時間：0.021 結果分析： 當選用復化梯形公

28、式時，對步長的事前估計所要求的步長很小，選取的節(jié)點很多，誤差絕對限要達到時，對不同的函數n 的取值需達到1000-10000之間，計算量是很大。 用復化simpson公式對步長的事前估計所要求的步長相對大些，選取的節(jié)點較少，誤差絕對限要達到時，對不同的函數n的取值只需在10-100之間，計算量相對小了很多，可滿足用較少的節(jié)點達到較高的精度，比復化梯形公式的計算量小了很多。用復化simpson公式計算所得的結果比用復化梯形公式計算所得的結果精度高很多，而且計算量小。 當選用Gauss-Lagrange I型公式進行計算時，選用較少的節(jié)點就可以達到很高的精度。實驗5.1 常微分方程性態(tài)和R-K法穩(wěn)

29、定性試驗實驗目的：考察下面微分方程右端項中函數y前面的參數對方程性態(tài)的影響（它可使方程為好條件的或壞條件的）和研究計算步長對R-K法計算穩(wěn)定性的影響。實驗容及要求： 實驗題目：常微分方程初值問題其中，。其精確解為。實驗要求：本實驗題都用4階經典R-K法計算。（1）對參數a分別取4個不同的數值：一個大的正值，一個小的正值，一個絕對值小的負值和一個絕對值大的負值。取步長h=0.01，分別用經典的R-K法計算，將四組計算結果畫在同一圖上，進行比較并說明相應初值問題的性態(tài)。（2）取參數a為一個絕對值不大的負值和兩個計算步長，一個步長使參數ah在經典R-K法的穩(wěn)定域，另一個步長在經典R-K法的穩(wěn)定域外。

30、分別用經典R-K法計算并比較計算結果。取全域等距的10個點上的計算值，列表說明。實驗程序：Matlab程序如下：function charp5RK%數值試驗5.1：常微分方程性態(tài)和R-K法穩(wěn)定性試驗%輸入：參數a，步長h%輸出：精確解和數值解圖形對比%clf;result=inputdlg('請輸入-50，50間的參數a:','實驗5.1',1,'-40');a=str2num(char(result); if (a<-50|a>50) errordlg('請輸入正確的參數a!'); return;endresult=inputdlg('請輸入（0 1）之間的步長:','實驗5.1'