利用Matlab實(shí)現(xiàn)Romberg數(shù)值積分算法----系統(tǒng)建模與仿真結(jié)課作業(yè)

1、利用Matlab實(shí)現(xiàn)Romberg數(shù)值積分算法一、內(nèi)容摘要針對(duì)于某些多項(xiàng)式積分，利用NewtonLeibniz積分公式求解時(shí)有困難，可以采用數(shù)值積分的方法，求解指定精度的近似解，本文利用Matlab中的.m文件編寫(xiě)了復(fù)化梯形公式與Romberg的數(shù)值積分算法的程序，求解多項(xiàng)式的數(shù)值積分，比較兩者的收斂速度。二、數(shù)值積分公式1.復(fù)化梯形公式求解數(shù)值積分的基礎(chǔ)是將區(qū)間一等分時(shí)的NewtonCotes求積公式:=其幾何意義是，利用區(qū)間端點(diǎn)的函數(shù)值、與端點(diǎn)構(gòu)成的梯形面積來(lái)近似在區(qū)間a,b上的積分值，截?cái)嗾`差為： 具有一次的代數(shù)精度，很明顯，這樣的近似求解精度很難滿足計(jì)算的要求，因而，可以采用將積分區(qū)間

2、不停地對(duì)分，當(dāng)區(qū)間足夠小的時(shí)候，利用梯形公式求解每一個(gè)小區(qū)間的積分近似值，然后將所有的區(qū)間加起來(lái)，作為被求函數(shù)的積分，可以根據(jù)計(jì)算精度的要求，劃分對(duì)分的區(qū)間個(gè)數(shù)，得到復(fù)化梯形公式：=其截?cái)嗾`差為： 2.Romberg數(shù)值積分算法使用復(fù)化的梯形公式計(jì)算的數(shù)值積分，其收斂速度比減慢，為此，采用Romberg數(shù)值積分。其思想主要是，根據(jù)的近似值加上與的近似誤差，作為新的的近視，反復(fù)迭代，求出滿足計(jì)算精度的近似解。用近似所產(chǎn)生的誤差可用下式進(jìn)行估算：新的的近似值： =（0 1 2 .）Romberg數(shù)值積分算法計(jì)算順序i=0(1)i=1(2)(3)i=2(4)(5)(6)i=3(7)(8)(9)(10

3、)i=4(11)(12)(13)(14)其中，第一列是二階收斂的，第二列是四階收斂的，第三列是六階收斂的，第四列是八階收斂的，即Romberg序列。 三、復(fù)化梯形法以及Romberg算法程序流程圖 圖1 復(fù)化梯形法程序流程圖圖2 Romberg算法程序流程圖四、計(jì)算實(shí)例依據(jù)上文所述的流程圖，編寫(xiě)復(fù)化梯形程序以及Romberg算法程序，并且利用實(shí)例驗(yàn)證程序的正確性，示例如下(計(jì)算精度)：表2 計(jì)算結(jié)果計(jì)算精度0.510-50.510-70.510-9復(fù)化梯形算法時(shí)間0.30.43.3近似值3.83.33.4Romberg算法時(shí)間0.00.70.8近似值3.83.43.2從上表中可以看出，當(dāng)要求的

4、計(jì)算精度不高時(shí)，復(fù)化梯形算法與Romberg算法計(jì)算時(shí)間相差不太大，但是Romberg算法是要快于復(fù)化梯形算法的；當(dāng)要求的計(jì)算精度更高的時(shí)候，Romberg算法是明顯快于復(fù)化梯形算法。本文所編寫(xiě)的程序適用于多項(xiàng)式的數(shù)值積分，且對(duì)于積分區(qū)間內(nèi)，被積函數(shù)在每一點(diǎn)必須有定義，在以后的學(xué)習(xí)中進(jìn)一步改進(jìn)。附錄：1.復(fù)化梯形算法程序function =sf(a,b,m,M,d)ticdisp(請(qǐng)輸入分子多項(xiàng)式a，分母多項(xiàng)式b，積分下限m,積分上限M,以及計(jì)算精度d)f=poly2sym(a)/poly2sym(b) %用于給用戶顯示被積函數(shù)的形式%利用梯形公式計(jì)算此數(shù)值積分disp(利用梯形公式計(jì)算數(shù)值積

5、分的結(jié)果)kk=zeros(); %用于存放結(jié)果kk(1,1)=1/2*(M-m)/1*(subs(f,x,m)+subs(f,x,M) %先存儲(chǔ)首項(xiàng)for i=1:1:230 t=0; for j=0:1:2(i-1)-1 v=m+(2*j+1)*(M-m)/(2i) vv=polyval(a,v)/polyval(b,v); t=t+(M-m)/(2i)*vv end y=1/2*kk(i,1)+t %通項(xiàng)公式計(jì)算各項(xiàng)值 kk(i+1,1)=y %存儲(chǔ)其他項(xiàng) f=i+1; %記錄符合條件的值的下標(biāo) if(1/3*(kk(i+1,1)-kk(i,1)=d) break; endendtime

6、=tocfprintf(The result is %fn, kk(f,1)2.Romberg算法程序function =romberg(a,b,m,M,d)ticdisp(請(qǐng)輸入分子多項(xiàng)式a，分母多項(xiàng)式b，積分下限m,積分上限M,以及計(jì)算精度d)f=poly2sym(a)/poly2sym(b) %用于給用戶顯示被積函數(shù)的形式disp(利用梯形公式計(jì)算數(shù)值積分的結(jié)果)kk=zeros(); %用于存放結(jié)果kk(1,1)=1/2*(M-m)/1*(subs(f,x,m)+subs(f,x,M); %先存儲(chǔ)首項(xiàng)for i=1:1:240 t=0; for j=0:1:2(i-1)-1 v=m+(2*j+1)*(M-m)/(2i); vv=polyval(a,v)/polyval(b,v); t=t+(M-m)/(2i)*vv; end y=1/2*kk(i,1)+t; %通項(xiàng)公式計(jì)算各項(xiàng)值 kk(i+1,1)=y ; %存儲(chǔ)其他項(xiàng) if(abs(1/3*(kk(i+1,1)-kk(i,1)=3) if(i+1=3 & abs(1/15*(kk(i+1,2)-kk(i,2)=4) if(i+1=4 & abs(1/63*(kk(i+1,3)-kk(i,3)=5) if(i+1=5 & abs