內(nèi)江師范學(xué)院數(shù)值分析實(shí)驗(yàn)三

1、內(nèi)江師范學(xué)院數(shù)值分實(shí)驗(yàn)報(bào)告冊(cè)編制 張莉 審定 牟廉明專業(yè)：班級(jí)：級(jí)班學(xué)號(hào)：姓名：數(shù)學(xué)與信息科學(xué)學(xué)院2013 年 9 月說明一、學(xué)生在做實(shí)驗(yàn)之前必須要準(zhǔn)備實(shí)驗(yàn)，主要包括預(yù)習(xí)與本次實(shí)驗(yàn)相關(guān)的 理論知識(shí)，熟練與本次實(shí)驗(yàn)相關(guān)的軟件操作，收集整理相關(guān)的實(shí)驗(yàn)參考資料， 要求學(xué)生在做實(shí)驗(yàn)時(shí)能帶上充足的參考資料；若準(zhǔn)備不充分，則學(xué)生不得參加 本次實(shí)驗(yàn)，不得書寫實(shí)驗(yàn)報(bào)告；二、要求學(xué)生要認(rèn)真做實(shí)驗(yàn)，主要是指不得遲到、早退和曠課，在做實(shí)驗(yàn) 過程中要嚴(yán)格遵守實(shí)驗(yàn)室規(guī)章制度，認(rèn)真完成實(shí)驗(yàn)內(nèi)容，極積主動(dòng)地向?qū)嶒?yàn)教 師提問等；三、各個(gè)實(shí)驗(yàn)按照學(xué)生水平分別設(shè)置了 A、B、C、D 四個(gè)等級(jí)，其中對(duì)應(yīng)的難度系數(shù)為 1、0.8、0

2、.7、0.6，也可根據(jù)實(shí)際完成情況制定相應(yīng)地的難度系數(shù)但總體保證難度排序?yàn)锳級(jí)難度最大，B級(jí)次之，C級(jí)較易，D級(jí)最簡(jiǎn)單。四、學(xué)生可以根據(jù)自己對(duì)各個(gè)實(shí)驗(yàn)涉及到的知識(shí)點(diǎn)掌握的程度自由選取A、B、 C、 D 等級(jí)的實(shí)驗(yàn)題目。五、學(xué)生要認(rèn)真工整地書寫實(shí)驗(yàn)報(bào)告，實(shí)驗(yàn)報(bào)告的內(nèi)容要緊扣實(shí)驗(yàn)的要求 和目的，不得抄襲他人的實(shí)驗(yàn)報(bào)告；四、根據(jù)實(shí)驗(yàn)準(zhǔn)備、實(shí)驗(yàn)態(tài)度、實(shí)驗(yàn)報(bào)告的書寫、實(shí)驗(yàn)報(bào)告的內(nèi)容進(jìn)行綜合評(píng)定，并給出實(shí)驗(yàn)成績(jī)?cè)u(píng)定分。實(shí)驗(yàn)名稱：實(shí)驗(yàn)三 數(shù)值積分指導(dǎo)教師： 吳開騰 張莉?qū)嶒?yàn)時(shí)數(shù)：2實(shí)驗(yàn)設(shè)備：安裝了 Matlab、C+、VF軟件的計(jì)算機(jī)實(shí)驗(yàn)日期：2013年10月_29_日 實(shí)驗(yàn)地點(diǎn)：第五教學(xué)樓北902實(shí)驗(yàn)?zāi)康?/p>

3、：掌握數(shù)值積分的基本思想和基本步驟。理解各類數(shù)值積分方法的優(yōu)缺點(diǎn)，并能自行編程求解。比較各類數(shù)值積分的代數(shù)精度，體會(huì)復(fù)化、變步長(zhǎng)及其加速的思想和實(shí)現(xiàn)步驟。 實(shí)驗(yàn)準(zhǔn)備：在開始本實(shí)驗(yàn)之前，請(qǐng)回顧教科書的相關(guān)內(nèi)容；需要一臺(tái)準(zhǔn)備安裝Windows XP Professional操作系統(tǒng)和裝有數(shù)學(xué)軟件的計(jì)算機(jī)。 實(shí)驗(yàn)內(nèi)容及要求A 題 計(jì)算積分I = j s exp(sinx x2 /100)dx。gB題 用Romberg法求函數(shù)積分I = j1sindx，精度為0.5e-6。0 xC 題 用復(fù)化梯形公式計(jì)算積分I Jdx01 + X2并分析剖分區(qū)間數(shù)對(duì)誤差的影響，取n二2,4,8,16,32,64,12

4、& 256,512積分精確值I =n = 3.141592653。D 題 計(jì)算 j30 2000udu。15 8.1u2 +1200說明：實(shí)驗(yàn)過程應(yīng)包括對(duì)問題的簡(jiǎn)要分析、求解方法、求解步驟、程序及其必要的 圖表等內(nèi)容。實(shí)驗(yàn)過程：A 題的實(shí)驗(yàn)過程1、實(shí)驗(yàn)中問題的重述計(jì)算積分 I 二exp(sinx x2 /100)dxg2、實(shí)驗(yàn)分析由于積分為I二卜exp(sinx x2/100)dx是從負(fù)無窮到正無窮積分，用常規(guī)的求積分的int函g數(shù)求解時(shí)，會(huì)顯示 Warning:Explicit integral could not be found。說明函數(shù) y=exp(sin(x)-x.A2/100)不是

5、 初等原函數(shù)，可用數(shù)值積分進(jìn)行分析。為方便計(jì)算，將積分區(qū)間分為兩個(gè)區(qū)間，以0 為分界點(diǎn)，分 為(g,0 )和(0, g)，則有：I =exp(sin x x2 /100)dx 二 f exp(sin x x2 /100)dx + f exp(sin x x2 /100)dxgg所以在兩個(gè)區(qū)間分別計(jì)算，然后將結(jié)果相加。3、實(shí)驗(yàn)求解過程首先利用辛普生法，從結(jié)果得到：在積分區(qū)間(0, g)積分上限增加，積分值趨于12.3599 ；在積 分區(qū)間(g,0)積分下限減小，積分值趨于10.0804；然后利用牛頓一柯特斯法，觀察結(jié)果得到，(0,g)的積分區(qū)間，積分上限增加，積分值趨于 12.3599 ；在(g

6、,0 )的積分區(qū)間，積分下限減小，積分值趨于10.0804；4、實(shí)驗(yàn)求解結(jié)果從上述兩種方法的實(shí)驗(yàn)結(jié)果表明，辛普生法和牛頓一柯特斯法在(0,g)的積分區(qū)間，積分上限 增加，積分值均趨于12.3599 ；在(g,0 )的積分區(qū)間，積分下限減小，積分值均趨于10.0804，故在 兩個(gè)積分區(qū)間I = f g exp(sinx x2/100)dx都是可積的，所以gI = f g exp(sin x 一 x2 /100)dx = f exp(sin x 一 x2 /100) dx + f exp(sin x 一 x2 /100) dxgg0故 I=10.0804+12.3599=22.4403B 題的實(shí)驗(yàn)

7、過程1、實(shí)驗(yàn)中問題的重述sin x用Romberg法求函數(shù)積分I = f 1dx，精度為05e - 60 x2、實(shí)驗(yàn)求解程序M文件function quad,R = Romberg(f,a,b,eps) %f表示被積函數(shù)%a, b表示積分區(qū)間的端點(diǎn) %eps表示精度%quad是用Romberg加速算法求得的積分值%只為Romberg表%err表示誤差的估計(jì) h=b_a;R(l,l)=h*(feval(f,a)+feval(f,b)/2;M=l;J=0;err=l; while erreps J=J+1; h=h/2; S=0; for p=1:M x=a+h*(2*p-1); S=S+feva

8、l(f,x); endR(J+1,1)=R(J,1)/2+h*S;M=2*M;for k=1:JR(J+1,k+1)=R(J+1,k)+(R(J+1,k)-R(J,k)/(4Ak-1); enderr=abs(R(J+1,J)-R(J+1,J+1);end quad=R(J+1,J+1)主程序 f=(x)(sin(x)./x);quad,R=Romberg(f,eps,1,0.5*10A(-6)3、實(shí)驗(yàn)求解結(jié)果用Romberg加速算法求得的積分值quad = 0.9461Romberg 表R 二0.9207000-0.93980.9461000.94450.94610.946100.94570

9、.94610.94610.9461C 題的實(shí)驗(yàn)過程1、實(shí)驗(yàn)中問題的重述14用復(fù)化梯形公式計(jì)算積分I J 1dx并分析剖分區(qū)間數(shù)對(duì)誤差的影響，取0 1 + X2n 二 2,4,8,16,32,64,128,256,512，積分精確值I 二兀二 3.141592653。2、實(shí)驗(yàn)分析在許多實(shí)際問題中，常常需要計(jì)算定積分I bf (x)dx的值。根據(jù)微積分學(xué)基本定理，若被積函a數(shù)f(x)在區(qū)間a,b上連續(xù)，只要能找到f(x)的一個(gè)原函數(shù)F(x)，便可利用牛頓-萊布尼茲公式 I bf (x) = F(b) - F(a)求得積分值。a 但是在實(shí)際使用中，往往遇到如下困難，而不能使用牛頓-萊布尼茲公式。找不

10、到用初等函數(shù)表示的原函數(shù)雖然找到了原函數(shù)，但因表達(dá)式過于復(fù)雜而不便計(jì)算f(x)是由測(cè)量或計(jì)算得到的表格函數(shù)由于以上種種困難，有必要研究積分的數(shù)值計(jì)算問題。利用插值多項(xiàng)式P (x)沁f (x)則積分1 bf (x)dx轉(zhuǎn)化為IbP (x)dx，顯然易算。JbP (x)dx稱為 n a a n a n 插值型求積公式。最簡(jiǎn)單的插值型求積公式是梯形公式和Simpson公式,。當(dāng)求積結(jié)點(diǎn)提供較多，可 以分段使用少結(jié)點(diǎn)的梯形公式和Simpson公式，并稱為復(fù)化梯形公式、復(fù)化Simpson公式。如步長(zhǎng)未知，可以通過誤差限的控制用區(qū)間逐次分半的策略自動(dòng)選取步長(zhǎng)的方法稱自適應(yīng)算法。 梯形遞推公式給出了區(qū)間分半

11、前后的遞推關(guān)系。由梯形遞推公式求得梯形序列，相鄰序列值作線性 組合得Simpson序列,Simpson序列作線性組合得柯特斯序列，柯特斯序列作線性組合的龍貝格序列。 若IR2-RIV，則輸出R2;否則依此類推。如此加工數(shù)據(jù)的過程叫龍貝格算法。3、實(shí)驗(yàn)結(jié)果利用龍貝格算法，得到當(dāng)分別取n利用龍貝格算法，得到當(dāng)分別取n二2,4,8,16,32,64,128,256,512時(shí)，產(chǎn)生的誤差見下表: 表 1 誤差結(jié)果表誤差-0.013239352830249-0.003315574025695-0.000828929586313 TOC o 1-5 h z n163264誤差-0.000207232961

12、173-0.000051808249116-0.000012952062417n128256512誤差-0.000003238015606-0.000000809503902-0.000000202375974D 題的實(shí)驗(yàn)過程1、實(shí)驗(yàn)分析302000u為了計(jì)算1du，設(shè)將積分區(qū)間劃分為n等份，則一共有n +1個(gè)等分點(diǎn)15 8.1u2 +1200 b-ax. = a + ih, h =，i = 0,1,n則可i從0到n -1累加求和，即可導(dǎo)出如下遞推算式:T2 廣 T2 廣 2 Tn + 2 空 f (:+1)hi=0+2b-a式中h二為二分前的步長(zhǎng)，而x= a + (i +ni+丄2、實(shí)驗(yàn)求解

13、過程已知被積函數(shù)f (X)，求積區(qū)間為a,b。首先，準(zhǔn)備初值:h u b - a, n u 1,T u 2 f (a) + f (b)1 h n -11然后，根據(jù)上式推導(dǎo)公式：有T二怎T +- f(a + (i + )h)2 2 1 22i=1b-a最后，利用復(fù)化梯形公式T =- n2 nf (a) + 4藝f (x ) + f (b)求積分。ii=1f f (a) + 4藝f (x ) + f (b)求出積分為:ii=13、實(shí)驗(yàn)結(jié)果b-a利用復(fù)化梯形公式T = 一n2 nJ 30152000u8.1u 2 +1200du =1.275067526872575e+002附錄A題的求解程序%辛普

14、生法clcclearfun=inline(exp(sin(x)-xJ2/100),x);Isim=quad(fun,0,100)Isim=quad(fun,0,1000)Isim=quad(fun,0,1000000)Isim=quad(fun,-100,0)Isim=quad(fun,-1000,0)Isim=quad(fun,-10000,0)%牛頓一柯特斯法clcclearfun=inline(exp(sin(x)-x.A2/100),x);IL=quadl(fun,0,100)IL=quadl(fun,0,1000)IL=quadl(fun,0,10000)IL=quadl(fun,-

15、10000,0)IL=quadl(fun,-1000,0)IL=quadl(fun,-100,0)B題的求解程序 function quad,R = Romberg(f,a,b,eps)%f表示被積函數(shù)%a, b表示積分區(qū)間的端點(diǎn) %eps表示精度%quad是用Romberg加速算法求得的積分值%只為Romberg表%err表示誤差的估計(jì)h=b_a;R(l,l)=h*(feval(f,a)+feval(f,b)/2;M=l;J=0;err=l;while errepsJ=J+1;h=h/2;S=0; for p=1:M x=a+h*(2*p-1); S=S+feval(f,x);endR(J+

16、1,1)=R(J,1)/2+h*S;M=2*M;for k=1:JR(J+1,k+1)=R(J+1,k)+(R(J+1,k)-R(J,k)/(4Ak-1); enderr=abs(R(J+1,J)-R(J+1,J+1);end quad=R(J+1,J+1)主程序：f=(x)(sin(x)./x);quad,R=Romberg(f,eps,1,0.5*10A(-6)C題的求解程序 function T = Complexification(a,b,Y,n)%復(fù)化梯形公式求解積分% a:下限% b：上限% n：點(diǎn)的個(gè)數(shù)c = 2*sum(Y(2:n-1);T = (b-a)/(2*n-2)*(Y(1)+c+Y(n);clcclear%主程序a = 0; b = 1;n = 2 4 8 16 32 64 128 256 512;T =;for i = 1:length(n)x = a:(b-a)/n(i):b;y = 4./(l+x.*x);tt = Complexification(a,b,y,length(x); T(i) = tt;end fprintf(分別產(chǎn)生的誤差為：) format longErr = (T - pi)./piD題的求解程序 function T = Complex