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

利用Matlab實(shí)現(xiàn)Romberg數(shù)值積分算法----系統(tǒng)建模與仿真結(jié)課作業(yè)利用Matlab實(shí)現(xiàn)Romberg數(shù)值積分算法----系統(tǒng)建模與仿真結(jié)課作業(yè)利用Matlab實(shí)現(xiàn)Romberg數(shù)值積分算法----系統(tǒng)建模與仿真結(jié)課作業(yè)xxx公司利用Matlab實(shí)現(xiàn)Romberg數(shù)值積分算法----系統(tǒng)建模與仿真結(jié)課作業(yè)文件編號：文件日期：修訂次數(shù)：第1.0次更改批準(zhǔn)審核制定方案設(shè)計，管理制度利用Matlab實(shí)現(xiàn)Romberg數(shù)值積分算法一、內(nèi)容摘要針對于某些多項(xiàng)式積分，利用Newton—Leibniz積分公式求解時有困難，可以采用數(shù)值積分的方法，求解指定精度的近似解，本文利用Matlab中的.m文件編寫了復(fù)化梯形公式與Romberg的數(shù)值積分算法的程序，求解多項(xiàng)式的數(shù)值積分，比較兩者的收斂速度。二、數(shù)值積分公式1.復(fù)化梯形公式求解數(shù)值積分的基礎(chǔ)是將區(qū)間一等分時的Newton—Cotes求積公式:=其幾何意義是，利用區(qū)間端點(diǎn)的函數(shù)值、與端點(diǎn)構(gòu)成的梯形面積來近似在區(qū)間[a,b]上的積分值，截斷誤差為：具有一次的代數(shù)精度，很明顯，這樣的近似求解精度很難滿足計算的要求，因而，可以采用將積分區(qū)間不停地對分，當(dāng)區(qū)間足夠小的時候，利用梯形公式求解每一個小區(qū)間的積分近似值，然后將所有的區(qū)間加起來，作為被求函數(shù)的積分，可以根據(jù)計算精度的要求，劃分對分的區(qū)間個數(shù)，得到復(fù)化梯形公式：=其截斷誤差為：數(shù)值積分算法使用復(fù)化的梯形公式計算的數(shù)值積分，其收斂速度比減慢，為此，采用Romberg數(shù)值積分。其思想主要是，根據(jù)的近似值加上與的近似誤差，作為新的的近視，反復(fù)迭代，求出滿足計算精度的近似解。用近似所產(chǎn)生的誤差可用下式進(jìn)行估算：新的的近似值：=（012….）Romberg數(shù)值積分算法計算順序i=0(1)i=1(2)(3)i=2(4)(5)(6)i=3(7)(8)(9)(10)i=4(11)(12)(13)(14)…………其中，第一列是二階收斂的，第二列是四階收斂的，第三列是六階收斂的，第四列是八階收斂的，即Romberg序列。三、復(fù)化梯形法以及Romberg算法程序流程圖圖1復(fù)化梯形法程序流程圖圖2Romberg算法程序流程圖四、計算實(shí)例依據(jù)上文所述的流程圖，編寫復(fù)化梯形程序以及Romberg算法程序，并且利用實(shí)例驗(yàn)證程序的正確性，示例如下(計算精度)：表2計算結(jié)果計算精度×10^-5×10^-7×10^-9復(fù)化梯形算法時間近似值3.3.3.Romberg算法時間近似值3.3.3.從上表中可以看出，當(dāng)要求的計算精度不高時，復(fù)化梯形算法與Romberg算法計算時間相差不太大，但是Romberg算法是要快于復(fù)化梯形算法的；當(dāng)要求的計算精度更高的時候，Romberg算法是明顯快于復(fù)化梯形算法。本文所編寫的程序適用于多項(xiàng)式的數(shù)值積分，且對于積分區(qū)間內(nèi)，被積函數(shù)在每一點(diǎn)必須有定義，在以后的學(xué)習(xí)中進(jìn)一步改進(jìn)。附錄：1.復(fù)化梯形算法程序function[]=sf(a,b,m,M,d)ticdisp('請輸入分子多項(xiàng)式a，分母多項(xiàng)式b，積分下限m,積分上限M,以及計算精度d')f=poly2sym(a)/poly2sym(b)%用于給用戶顯示被積函數(shù)的形式%利用梯形公式計算此數(shù)值積分disp('利用梯形公式計算數(shù)值積分的結(jié)果')kk=zeros();%用于存放結(jié)果kk(1,1)=1/2*(M-m)/1*(subs(f,'x',m)+subs(f,'x',M))%先存儲首項(xiàng)fori=1:1:2^30t=0;forj=0:1:2^(i-1)-1v=m+(2*j+1)*(M-m)/(2^i)vv=polyval(a,v)/polyval(b,v);t=t+(M-m)/(2^i)*vvendy=1/2*kk(i,1)+t%通項(xiàng)公式計算各項(xiàng)值kk(i+1,1)=y%存儲其他項(xiàng)f=i+1;%記錄符合條件的值的下標(biāo)if(1/3*(kk(i+1,1)-kk(i,1))<=d)break;endendtime=tocfprintf('Theresultis%f\n',kk(f,1))算法程序function[]=romberg(a,b,m,M,d)ticdisp('請輸入分子多項(xiàng)式a，分母多項(xiàng)式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));%先存儲首項(xiàng)fori=1:1:2^40t=0;forj=0:1:2^(i-1)-1v=m+(2*j+1)*(M-m)/(2^i);vv=polyval(a,v)/polyval(b,v);t=t+(M-m)/(2^i)*vv;endy=1/2*kk(i,1)+t;%通項(xiàng)公式計算各項(xiàng)值kk(i+1,1)=y;%存儲其他項(xiàng)if(abs(1/3*(kk(i+1,1)-kk(i,1)))<=d)%判斷梯形公式值是否達(dá)到要求disp('Theresultis:')kk()kk(i+1,1)%梯形值滿足要求，輸出結(jié)果break;elses=(4*kk(i+1,1)-kk(i,1))/(4-1);%構(gòu)造simpson各項(xiàng)kk(i+1,2)=s%存儲if(i+1>=3)if(i+1>=3&abs(1/15*(kk(i+1,2)-kk(i,2)))<=d)kk()disp('Theresultis:')kk(i+1,2)%simpson值滿足要求，輸出結(jié)果pan1=0;break;elsec=(4^2*kk(i+1,2)-kk(i,2))/(4^2-1);%構(gòu)造cotes值kk(i+1,3)=c%存儲cotes值if(i+1>=4)if(i+1>=4&abs(1/63*(kk(i+1,3)-kk(i,3)))<=d)disp('Theresultis:')kk(i+1,3)break;elser=(4^3*kk(i+1,3)-kk(i,3))/(4^3-1)%構(gòu)造romberg值kk(i+1,4)=r%存儲romberg值if(i+1>=5)if(i+1>=5&abs(1/127*(kk(i+1,4)-kk(i,4)))<=d)disp('Theresultis:')