數(shù)值分析實驗報告

1、Xi= -1 + 2i/n , i=0 , 1, 2,，n, 實驗2.1多項式插值的振蕩現(xiàn)象 實驗?zāi)康模?Ln(X) 25x； li(x) 逼近的函數(shù)。 實驗內(nèi)容： Runge給出的一個例子是極者名并富有啟發(fā)性的。 在一個固定的區(qū)間上用插值逼近一個函數(shù)，顯然Lagrange插值中使用的節(jié)點越多，插 值多項式的次數(shù)就越高。我們自然關(guān)心插值多項式的次數(shù)增加時，Ln(x)是否也更加靠近被 設(shè)區(qū)間卜1，1上函數(shù)f(x)=1/(1+25x 2)。 考慮區(qū)間-1,1的一個等距劃分，分點為 則拉格朗日插值多項式為 其中，li(x), i=0 , 1, 2，n是n次Lagrange插值基函數(shù)。 實驗步驟與結(jié)果

2、分析： 實驗源程序 function Chap2 In terpolati on %數(shù)值實驗二：“實驗2.1 :多項式插值的震蕩現(xiàn)象” %輸入：函數(shù)式選擇，插值結(jié)點數(shù) %輸出：擬合函數(shù)及原函數(shù)的圖形 promps = 請選擇實驗函數(shù)，若選f(x),請輸入f,若選h(x),請輸入h,若選g(x)，請輸入g:; titles = charpt_2; result = in putdlg(promps,charpt 2,1,f); Nb_f = char(result); if(Nb_f = f return;e nd result = inputdlg(請輸入插值結(jié)點數(shù)N:,charpt_2,1,

3、10); Nd = str2 nu m(char(result); if(Nd 1)errordlg(結(jié)點輸入錯誤！);return;end switch Nb_f case f f=in li ne(1丿(1+25*x42); a = -1;b = 1; case h f=i nlin e(x./(1+x.A4); a = -5; b = 5; case g f=in li ne(ata n( x); a = -5; b= 5; end x0 = lin space(a, b, Nd+1); y0 = feval(f, x0); x = a:0.1:b; y = Lagra nge(xO,

4、y0, x); fplot(f, a b, co); hold on; plot(x, y, b-); xlabel(x); ylabel(y = f(x) o and y = Ln(x)-); % function y=Lagra nge(xO, 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 實

5、驗結(jié)果分析 增大分點n=2, 3，時，拉格朗日插值函數(shù)曲線如圖所示。 n=6 n=7 n=8 n=9 n=10 從圖中可以看出，隨著n的增大，拉格朗日插值函數(shù)在x=0附近較好地逼近了原來的 函數(shù) f(x) ，但是卻在兩端 x= -1 和 x=1 處出現(xiàn)了很大的振蕩現(xiàn)象。 并且，仔細分析圖形，可以看出，當n為奇數(shù)時，雖然有振蕩，但振蕩的幅度不算太大， n為偶數(shù)時，其振蕩幅度變得很大。通過思考分析，我認為，可能的原因是f(x)本身是偶函 數(shù)，如果n為奇數(shù)，那么Lagrange插值函數(shù)Ln(x)的最高次項xn-1是偶次幕，比較符合f(x) 本身是偶函數(shù)的性質(zhì)；如果 n為偶數(shù)，那么Lagrange插值

6、函數(shù)Ln(x)的最高次項xn-1是奇次 冪，與 f(x) 本身是偶函數(shù)的性質(zhì)相反，因此振蕩可能更劇烈。 將原來的f(x)換為其他函數(shù)如h(x)、g(x)，結(jié)果如圖所示。 其中 h(x), g(x) 均定義在 -5， 5區(qū)間上， h(x)=x/(1+x 4)， g(x)=arctan x 。 h(x), n=7 h(x), n=8 h(x), n=9 h(x), n=10 g(x), n=7 g(x), n=8 g(x), n=9 g(x), n=10 分析兩個函數(shù)的插值圖形，可以看出： 隨著n的增大，拉格朗日插值函數(shù)在x=0附近較好地逼近了原來的函數(shù)f(x)，但是卻在 兩端x= -5和x=5處

7、出現(xiàn)了很大的振蕩現(xiàn)象。 并且，仔細分析圖形，可以看出，當n為偶數(shù)時，雖然有振蕩，但振蕩的幅度不算太大， n為奇數(shù)時，其振蕩幅度變得很大。原因和上面f(x)的插值類似，h(x)、g(x)本身是奇函數(shù)， 如果n為偶數(shù)，那么Lagrange插值函數(shù)Ln(x)的最高次項xn-1是奇次幕，比較符合h(x)、g(x) 本身是奇函數(shù)的性質(zhì)；如果n為奇數(shù)，那么Lagrange插值函數(shù)Ln(x)的最高次項xn-1是偶次 幕，與h(x)、g(x)本身是奇函數(shù)的性質(zhì)相反，因此振蕩可能更劇烈。 實驗3.1多項式最小二乘擬合 實驗?zāi)康模?* lx- kx中的參數(shù) ak、平方誤差3 2,并作離散據(jù) 編制以函數(shù)xkk=O,

8、n；為基的多項式最小二乘擬合程序。 Xi -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 yi -4.447 -0.452 0.551 0.048 -0.447 0.549 4.552 實驗內(nèi)容： 對表中的數(shù)據(jù)作三次多項式最小二乘擬合。 取權(quán)函數(shù)Wi= 1,求擬合曲線 n xi, yi的擬合函數(shù)的圖形。 實驗源程序 fun ctio n Chap3CurveFitt ing alph和誤差r %數(shù)值實驗三：“實驗3.1” %輸出:原函數(shù)及求得的相應(yīng)插值多項式的函數(shù)的圖像以及參數(shù) x0 = -1:0.5:2; yO = -4.447 -0.452 0.551 0.048 -0.447

9、 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(x); ylabel(y0 * and polyfit.y-); hold on plot(x0, y0, *) grid on; disp(平方誤差:,num2str(r) disp(參數(shù) alph:, num2str(alph) 實驗結(jié)果 平方誤差：2.1762

10、e-005 參數(shù) alph:1.9991-2.9977 -3.9683e-0050.54912 實驗4.1 實驗?zāi)康模簭突蠓e公式計算定積分. 實驗題目：數(shù)值計算下列各式右端定積分的近似值 1)In 2 - In 3 = 2 rdx;(2)兀二 4一rrfr : 實驗要求: (1) 若用復化梯形公式、復化Simpson公式和復化 Gauss-Legendre I型公式做計算，要求 1 * 10 7 .、. 絕對誤差限為$，分別利用它們的余項對每種算法做出步長的事前估計 (2) 分別用復化梯形公式，復化Simpson公式和復化 Gauss-Legendre I型公式作計算. (3) 將計算結(jié)果與

11、精確解做比較，并比較各種算法的計算量 實驗程序： 1. 事前估計的 Matlab程序如下： (1).用復化梯形公式進行事前估計的Matlab程序 format long g x=2:0.01:3; f=-4*(3*x.A2+1)./(x.A2-1).A3; % 二階導函數(shù) %plot(x,f)%畫出二階導函數(shù)圖像 x=2.0;%計算導函數(shù)最大值 f=-4*(3*xA2+1)/(xA2-1)A3; h2=0.5*10A(-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.A

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

13、=ceil(1/h)+1 format long g x=1:0.01:2; %二階導函數(shù) %畫出二階導函數(shù)圖像 %計算導函數(shù)最大值 %步長 %選取的點數(shù) f=2.*exp(x)+x.*exp(x);% %plot(x,f) x=2; 二階導函數(shù) %畫出二階導函數(shù)圖像 %計算導函數(shù)最大值 f=2.*exp(x)+x.*exp(x); %步長 h2=0.5*10A(-7)*12/f; h=sqrt(abs(h2) n=1/h n=ceil(1/h)+1%選取的點數(shù) 估計結(jié)果步長 h 及結(jié)點數(shù) n 分別為 h = 0.000558156305651438 n =1793 h =0.00054772

14、2557505166 n =1827 h =0.000407071357304889 n =2458 h =0.000142479094906909 n =7020 (2)用復化 simpson 公式進行事前估計的 Matlab 程序 format long g x=2:0.01:3; f=-2*(-72*x.A2-24).*(x.A2-1)-192*x.A2.*(x.A2+1)./(x.A2-1).A5;% 四階導函數(shù) x=2.0; f=-2*(-72*xA2-24)*(xA2-1)-192*xA2*(xA2+1)/(xA2-1)A5; % 計算導函數(shù)最大值 h4=0.5*10A(-7)*1

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

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

17、值 h4=0.5*10A(-7)*180*16/f; h=sqrt(sqrt(abs(h4) n=1/h;%求分段區(qū)間個數(shù) n=2*ceil(1/h)+1 % 選取的點數(shù) 估計結(jié)果步長 h 及結(jié)點數(shù) n 分別為 h =0.0437490486013411 n =47 h =0.0588566191276542 n =35 h =0.0757645166218433 n =29 h =0.0424527247118546 n =49 2. 積分計算的 Matlab 程序： format long g promps=請選擇積分公式，若用復化梯形，請輸入 T,用復化simpson,輸入S, 用復化

18、Gauss_Legendre，輸入 GL : ; result=inputdlg(promps,charpt 4,1,T); Nb=char(result); if(Nb=T return; end result=inputdlg( 請輸入積分式題號 1-4： ,實驗 4.1,1,1); Nb_f=str2num(char(result); if(Nb_f4) errordlg( 沒有該積分式 ); return; end switch Nb_f case 1 fun=i nlin e(-2./(x.A2-1);a=2;b=3; case 2 fun=inline(4./(x.A2+1);a=

19、0;b=1; case 3 fun=inline(3.Ax);a=0;b=1; case 4 fun=inline(x.*exp(x);a=1;b=2; end if(Nb=T)% 用復化梯形公式 promps= 請輸入用復化梯形公式應(yīng)取的步長： ; result=inputdlg(promps, 實驗 4.2,1,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

20、+fun(xk); end t=h*(fun(a)+fun(b)+2*detsum)/2; time=toc; t end if(Nb=S)% 用復化 Simpson 公式 promps= 請輸入用復化 Simpson 公式應(yīng)取的步長： ; result=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; de

21、tsum_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=h*(fun(a)+fun(b)+2*detsum_1+4*detsum_2)/6; time=toc; t end if(Nb=GL)% 用復化 Gauss_Legendre I %先根據(jù)復化 Gauss_Legendre I 公式的余項估計步長 promps= 請輸入用復化 Gauss_Legendre I 公式應(yīng)取的步長： ; result=inputdlg(promps, 實驗 4.2,1,0.0

22、1); h=str2num(char(result); if(h=0) errordlg( 請輸入正確的步長 !); 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(精確解：In 2-1 n3=-0.4054651081) disp(絕對誤差：,n um2str(abs(t+0.4054651081); disp( 運行時

23、間：,num2str(time); case 2 disp(精確解：pi=3.14159265358979) disp( 絕對誤差：,num2str(abs(t-pi); disp( 運行時間：,num2str(time); case 3 disp(精確解：2/l n3=1.82047845325368) disp( 絕對誤差： ,num2str(abs(t-1.82047845325368); disp( 運行時間： ,num2str(time); case 4 disp(精確解：eA2=7.38905609893065) disp(絕對誤差：,nu m2str(abs(t-7.389056

24、09893065); disp( 運行時間： ,num2str(time); end 1. 當選用復化梯形公式時： (1)式運行結(jié)果為： t =-0.40546512204351 精確解： ln2-ln3=-0.4054651081 絕對誤差： 1.3944e-008 運行時間： 0.003 (2)式運行結(jié)果為： t =3.14159261385336 精確解： pi=3.14159265358979 絕對誤差： 3.9736e-008 運行時間： 0.005 ( 3)式運行結(jié)果為： t = 1.82047849690861 精確解： 2/ln3=1.82047845325368 絕對誤差：

25、4.3655e-008 運行時間： 0.016 (4)式運行結(jié)果為： t =7.38905611970610 精確解：eA2=7.38905609893065 絕對誤差： 2.0775e-008 運行時間： 0.007 2. 當選用復化 Simpson 公式進行計算時： (1)式運行結(jié)果為： t =-0.405465108127519 精確解： ln2-ln3=-0.4054651081 絕對誤差： 2.7519e-011 運行時間： 0.022 (2)式運行結(jié)果為： t =3.14159265358979 精確解： pi=3.14159265358979 絕對誤差： 0 運行時間： 0.02

26、1 ( 3)式運行結(jié)果為： t =1.82047845326288 精確解： 2/ln3=1.82047845325368 絕對誤差： 9.2018e-012 運行時間： 0.019 （4 ）式運行結(jié)果為： t =7.38905609902118 精確解：eA2=7.38905609893065 絕對誤差：9.0528e-011 運行時間：0.021 3. 當選用復化 Gauss-Legendre I型公式進行計算時： （1 ）式運行結(jié)果為： t =-0.405465108095262 精確解：In 2-l n3=-0.4054651081 絕對誤差：4.7385e-012 運行時間：0.02

27、3 （2 ）式運行結(jié)果為： t =3.14159265358979 精確解：pi=3.14159265358979 絕對誤差：1.3323e-015 運行時間：0.021 （3 ）式運行結(jié)果為： t =1.82047845324754 精確解：2/ln3=1.82047845325368 絕對誤差：6.1431e-012 運行時間：0.019 （4 ）式運行結(jié)果為： t =7.38905607315046 精確解：eA2=7.38905609893065 絕對誤差：1.441e-012 運行時間：0.021 結(jié)果分析： 當選用復化梯形公式時，對步長的事前估計所要求的步長很小，選取的節(jié)點很多，誤

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

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

30、絕對值不大的負值和兩個計算步長， 一個步長使參數(shù)ah在 經(jīng)典 R-K 法的穩(wěn)定域內(nèi)，另一個步長在經(jīng)典 R-K 法的穩(wěn)定域外。 分別用經(jīng)典 R-K 法計算并比較計算結(jié)果。取全域等距的 10 個點上的計算值，列表說明。 實驗程序： Matlab 程序如下： function charp5RK %數(shù)值試驗 5.1：常微分方程性態(tài)和 R-K 法穩(wěn)定性試驗 %輸入：參數(shù)a，步長h %輸出：精確解和數(shù)值解圖形對比 %clf; result=inputdlg( 請輸入 -50， 50間的參數(shù) a:, 實驗 5.1,1,-40); a=str2num(char(result); if (a50)errordlg(請輸入正確的參數(shù) a!); return;e nd result=inputdlg( 請輸入( 0 1)之間的步長 :, 實驗 5.1,1,0.01); h=str2num(char(r