計(jì)算方法第二次上機(jī)作業(yè)

1、計(jì)算方法第二次上機(jī)作業(yè)Gegebao摘要：程序基于MATLAB，包括問(wèn)題陳述、算法與程序、結(jié)果與分析、討論四個(gè)部分。一 、問(wèn)題陳述 數(shù)學(xué)上已經(jīng)證明了：0141+x2dx=成立，所以可以通過(guò)積分來(lái)計(jì)算的近似值。（1） 分別使用矩形、梯形和Simpson復(fù)合求積公式計(jì)算的近似值。選擇不同的h，對(duì)于每種求積公式，試將誤差刻畫(huà)成h的函數(shù)，并比較各方面的精度。是否存在某個(gè)h值，當(dāng)?shù)陀谶@個(gè)值之后，再繼續(xù)減少h的值，計(jì)算不再有所改進(jìn)？為什么？（2） 實(shí)現(xiàn)Romberg求積方法，并重復(fù)上面的計(jì)算。（3） 使用自適應(yīng)求積方法重復(fù)上面的計(jì)算。二、 算法與程序1矩形求積法考慮到對(duì)于具體的某一個(gè)h，不一定能整數(shù)個(gè)地覆

2、蓋積分空間，將會(huì)嚴(yán)重地影響算法精度。于是我們用n代替h，即將積分區(qū)間劃分為n個(gè)區(qū)間，h=1/n。函數(shù)將輸出積分結(jié)果Rec和這個(gè)結(jié)果對(duì)于精確的的誤差的對(duì)數(shù)值r。function Rec r = rec( n ) h=1/n; %求出對(duì)應(yīng)的hRec=0;for i=1:n x0=h*(i-1/2); %求出每一區(qū)間的中心點(diǎn)的橫坐標(biāo) Rec=Rec+h*4/(1+x02);endr=log10(abs(Rec-pi); %取其誤差的對(duì)數(shù)end2.梯形求積法 和矩形法一樣，采用n來(lái)標(biāo)度取點(diǎn)密度。同樣輸出積分結(jié)果Lad和這個(gè)結(jié)果對(duì)于精確的的誤差的對(duì)數(shù)值r。function Lad r = lad( n)

3、 h=1/n;Lad=0;for i=1:n x0=h*(i-1); x1=x0+h; 取得每個(gè)區(qū)間的端點(diǎn)的橫坐標(biāo) Lad=Lad+h/2*(4/(1+x02)+4/(1+x12);endr=log10(abs(Lad-pi);end3.Simpson求積法 輸入取點(diǎn)數(shù)目n，h=1/n, 輸出積分結(jié)果Simp和這個(gè)結(jié)果對(duì)于精確的的誤差的對(duì)數(shù)值r。function Simp ,r = simp( n)h=1/n; Simp=0;for i=1:n x0=h*(i-1); x1=x0+h/2; x2=x1+h/2; Simp=Simp+h/6*(4/(1+x02)+4*4/(1+x12)+4/(1

4、+x22);endr=log10(abs(Simp-pi);end4.Romberg求積法 輸入最初取點(diǎn)數(shù)目n（h=1/n）和需求遞推下去的步數(shù)k, 輸出積分結(jié)果Rom和這個(gè)結(jié)果對(duì)于精確的的誤差的對(duì)數(shù)值r。function Rom,r = rom(n,k)for i=1:k; a(i,1) c=lad(2(i-1)*n); %調(diào)用lad函數(shù)，求出遞推矩陣第一列的數(shù)值endfor i=2:k for j=1:k-i+1 a(j,i)=(a(j+1,i-1)-4(1-i)*a(j,i-1)/(1-4(1-i); endendRom=a(1 ,k);r=log10(abs(Rom-pi);end5.

5、分別基于前幾種算法的自適應(yīng)算法（1）基于矩形法 輸入初始的取點(diǎn)數(shù)n（h=1/n）和所需要的精度的對(duì)數(shù)值e，輸出積分結(jié)果對(duì)于精確的的誤差的對(duì)數(shù)值r。 function r = Arec(n,e)temp0=0;temp1=1; %temp0為T(mén)（h），Temp1為T(mén)（h/2）while(log10(abs(temp0-temp1)>e) temp0=temp1; n=n*2; temp1=rec(n);endr=log10(abs(temp1-pi);end （2）基于梯形法 輸入初始的取點(diǎn)數(shù)n（h=1/n）和所需要的精度的對(duì)數(shù)值efunction r = Alad( n,e)temp0=

6、0;temp1=1;while(log10(abs(temp0-temp1)>e) temp0=temp1; n=n*2; temp1=lad(n);endr=log10(abs(temp1-pi);end（3）基于Simpson法 輸入初始的取點(diǎn)數(shù)n（h=1/n）和所需要的精度的對(duì)數(shù)值efunction r = Asimp( n,e)temp0=0;temp1=1;while(log10(abs(temp0-temp1)>e) temp0=temp1; n=n*2; temp1=simp(n);endr=log10(abs(temp1-pi);end（4）基于Romberg法 輸

7、入初始的取點(diǎn)數(shù)n（h=1/n）和所需要的精度的對(duì)數(shù)值e，取遞推步數(shù)k=5function r = Arom( n,e)temp0=0;temp1=1;while(log10(abs(temp0-temp1)>e) temp0=temp1; n=n*2; temp1=rom(n,5);endr=log10(abs(temp1-pi);end三、結(jié)果與分析使用之前的積分函數(shù)進(jìn)行運(yùn)算?？紤]到h的變化范圍比較大，將h等比例地從1取到10(-7),取200個(gè)點(diǎn)。最后將h和誤差按對(duì)數(shù)輸出成圖像。1.矩形積分法運(yùn)行代碼： for i=1:200 h0(i)=-i*7/200; n=fix(10(-h0

8、(i); r Rrec(i)=rec(n); endplot(h0,Rrec,'b');分析：求積方法的最高精度可以達(dá)到10(-14)到10(-15)之間，而當(dāng)h取到10(-6)時(shí)，繼續(xù)提高取點(diǎn)密度對(duì)算法精度沒(méi)有改進(jìn)。而在此之前，精度和h的關(guān)系幾乎是線性的。2.梯形求積法運(yùn)行代碼：for i=1:200 h0(i)=-i*7/200; n=fix(10(-h0(i); r Rlad(i)=lad(n);endplot(h0,Rlad,'b');分析：求積方法的最高精度可以達(dá)到10(-14)到10(-15)之間，而當(dāng)h取到10(-6)時(shí)，精度隨h取小開(kāi)始波動(dòng)，繼續(xù)

9、減小h并不能保證精度的增加。而在此之前，精度和h的關(guān)系幾乎是線性的。3.Simpson求積法運(yùn)行代碼：for i=1:200 h0(i)=-i*7/200; n=fix(10(-h0(i) r Rsimp(i)=simp(n);endplot(h0,Rsimp,'b');分析：求積方法的最高精度可以達(dá)到10(-15)到10(-16)之間，而當(dāng)h取到10(-2)時(shí)，精度隨h取小開(kāi)始波動(dòng)，繼續(xù)減小h并不能保證精度的增加。而在此之前，精度和h的關(guān)系幾乎是線性的。圖像上存在間斷點(diǎn)，實(shí)際上是由于求積方法的精度已經(jīng)超過(guò)MATLAB所給出的值，計(jì)算機(jī)得到的誤差為零，求對(duì)數(shù)運(yùn)算無(wú)意義，故出現(xiàn)間

10、斷點(diǎn)。4.三種算法的精度比較在進(jìn)行前三個(gè)積分法的運(yùn)算的時(shí)候，已經(jīng)分別記錄了結(jié)果，直接輸出就行。運(yùn)行代碼：plot(h0,Rrec,'b');hold on;plot(h0,Rlad,'r');hold on;plot(h0,Rsimp,'k');hold on;分析：對(duì)于相同的h值，在很長(zhǎng)的范圍內(nèi)，Simpson的精度明顯高于另兩種算法，而矩形法又比梯形法略?xún)?yōu)。三種算法的極限精度以Simpson為最好，可以達(dá)到10（-15）到10（-16）之間。5Romberg積分法（固定遞推步數(shù)為5，基于梯形算法）運(yùn)行代碼：k=5for i=1:200 h0(

11、i)=-i*7/200; n=fix(10(-h0(i) r Rrom(i)=rom(n,k);endplot(h0,Rrom,'b');分析：可以看出，積分精度隨h減小提高地非?？?。在（-1,0）這個(gè)區(qū)間里，便出現(xiàn)了間斷點(diǎn)，說(shuō)明算法的精度已經(jīng)達(dá)到MATLAB所能達(dá)到的精度，從而誤差變?yōu)榱?，求?duì)數(shù)運(yùn)算無(wú)意義，出現(xiàn)間斷點(diǎn)。縱觀函數(shù)，最高精度應(yīng)在10(-15)到10(-16)之間。6.自適應(yīng)算法（取自適應(yīng)精度為10(-10)）：（1）基于矩形算法 運(yùn)行代碼：for i=1:200 h0(i)=-i*7/200; n=fix(10(-h0(i) Rarec(i)=Arec(n,-10

12、);endplot(h0,Rarec,'b');(2)基于梯形算法：運(yùn)行代碼：for i=1:200 h0(i)=-i*7/200; n=fix(10(-h0(i) Ralad(i)=Alad(n,-10);endplot(h0,Ralad,'b');（3）基于Simpson算法 運(yùn)行代碼：for i=1:200 h0(i)=-i*7/200; n=fix(10(-h0(i) Rasimp(i)=Asimp(n,-10);endplot(h0,Rasimp,'b');（4）基于Romberg算法（ 固定遞推步數(shù)為五步） 運(yùn)行代碼：for i=1:

13、200 h0(i)=-i*7/200; n=fix(10(-h0(i) Rarom(i)=Arom(n,-10);endplot(h0,Rarom,'b');分析：當(dāng)本來(lái)的算法精度低于所給出的自適應(yīng)精度時(shí)，自適應(yīng)算法對(duì)本來(lái)算法的精度是有改進(jìn)作用的，并使算法精度略高于所給出的自適應(yīng)精度。但是當(dāng)算法對(duì)于具體的h得到的結(jié)果的精度高于自適應(yīng)給出的精度時(shí)，自適應(yīng)算法則對(duì)原本的算法無(wú)改進(jìn)作用。而且有些算法的精度隨h的減小提高得特別快（比如說(shuō)Romberg算法），使用自適應(yīng)算法后算法精度會(huì)明顯高于給出的自適應(yīng)精度。四、討論 算法的精度主要取決于算法本身、h的大小和機(jī)器誤差。幾種算法的精確程度不同，但都要受制于機(jī)器本身的誤差，不可能無(wú)限地提高精度。當(dāng)然，并不是說(shuō)