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

2、不停地對分，當(dāng)區(qū)間足夠小的時候，利用梯形公式求解每一個小區(qū)間的積分近似值，然后將所有的區(qū)間加起來，作為被求函數(shù)的積分，可以根據(jù)計算精度的要求，劃分對分的區(qū)間個數(shù)，得到復(fù)化梯形公式：=其截斷誤差為： 2.Romberg數(shù)值積分算法使用復(fù)化的梯形公式計算的數(shù)值積分，其收斂速度比減慢，為此，采用Romberg數(shù)值積分。其思想主要是，根據(jù)的近似值加上與的近似誤差，作為新的的近視，反復(fù)迭代，求出滿足計算精度的近似解。用近似所產(chǎn)生的誤差可用下式進(jìn)行估算：新的的近似值： =（0 1 2 .）Romberg數(shù)值積分算法計算順序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ù)上文所述的流程圖，編寫復(fù)化梯形程序以及Romberg算法程序，并且利用實例驗證程序的正確性，示例如下(計算精度)：表2 計算結(jié)果計算精度0.510-50.510-70.510-9復(fù)化梯形算法時間0.30.43.3近似值3.83.33.4Romberg算法時間0.00.70.8近似值3.83.43.2從上表中可以看出，當(dāng)要求的

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

5、分的結(jié)果)kk=zeros(); %用于存放結(jié)果kk(1,1)=1/2*(M-m)/1*(subs(f,x,m)+subs(f,x,M) %先存儲首項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 %通項公式計算各項值 kk(i+1,1)=y %存儲其他項 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(請輸入分子多項式a，分母多項式b，積分下限m,積分上限M,以及計算精度d)f=poly2sym(a)/poly2sym(b) %用于給用戶顯示被積函數(shù)的形式disp(利用梯形公式計算數(shù)值積分的結(jié)果)kk=zeros(); %用于存放結(jié)果kk(1,1)=1/2*(M-m)/1*(subs(f,x,m)+subs(f,x,M); %先存儲首項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; %通項公式計算各項值 kk(i+1,1)=y ; %存儲其他項 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