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

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

2、( )(a,b)12具有一次的代數(shù)精度，很明顯，這樣的近似求解精度很難滿足計(jì)算的要求，因而，可以采用將積分區(qū)間不停地對(duì)分，當(dāng)區(qū)間足夠小的時(shí)候，利用梯形公式求解每一個(gè)小區(qū)間的積分近似值，然后將所有的區(qū)間加起來(lái)，作為被求函數(shù)的積分， 可以根據(jù)計(jì)算精度的要求，劃分對(duì)分的區(qū)間個(gè)數(shù)，得到復(fù)化梯形公式：nI = " f (x)dx (b a)f(a) f(b) f(a 血包)a2nkdn其截?cái)嗾`差為：r = (b a”2 f“- (a,b)122. Romberg數(shù)值積分算法使用復(fù)化的梯形公式計(jì)算的數(shù)值積分，其收斂速度比減慢，為此，采用Romberg數(shù)值積分。其思想主要是，根據(jù)I的近似值T2n加

3、上I與T2n的近似誤差,作為新的I的近視，反復(fù)迭代，求出滿足計(jì)算精度的近似解。用T2n近似I所產(chǎn)生的誤差可用下式進(jìn)行估算：1 : = I - Tgn(Tn - T?n l)3新的I的近似值：T2k=mT2nj = (0 1 2 .)Romberg數(shù)值積分算法計(jì)算順序i=0(1)T；i=1(2)T20(3)Ti=2T20(5)T：(6)T22i=3(7)T(8)T212(9)T22(10)T23i=4(11)T20(12)T；(13)T22(14)T23其中，第一列是二階收斂的，第二列是四階收斂的，第三列是六階收斂的,第四列是八階收斂的，即Romberg序列。三、復(fù)化梯形法以及Romberg算法

4、程序流程圖圖1復(fù)化梯形法程序流程圖輸入被積函數(shù)、積分區(qū)間、收斂條件土差值是否滿足收斂條件構(gòu)建并存儲(chǔ) T1i+i i是否大于等于v差值是否滿足收斂條件構(gòu)建并存儲(chǔ) T 2ii是否大于等于3差值是否滿足收斂條件3構(gòu)建并存儲(chǔ)T 3ii是否大于等于差值是否滿足收斂條件輸岀結(jié)果圖2 Romberg算法程序流程圖四、計(jì)算實(shí)例依據(jù)上文所述的流程圖，編寫(xiě)復(fù)化梯形程序以及 Romberg算法程序，并且利 用實(shí)例驗(yàn)證程序的正確性，示例如下（計(jì)算精度）：兀2dx1 x2表2計(jì)算結(jié)果計(jì)算精度0.5 X 10A-50.5 X 10A-70.5 X 10A-9復(fù)化梯形時(shí)間0.0698263946330.2166358023

5、043.459824945493算法近似值3.1415901104583.1415926138533.141592653434Romberg時(shí)間0.0456873297100.0433617263570.044913907518算法近似值3.1415925024583.1415926512243.141592653552從上表中可以看出，當(dāng)要求的計(jì)算精度不高時(shí)，復(fù)化梯形算法與Romberg算法計(jì)算時(shí)間相差不太大，但是Romberg算法是要快于復(fù)化梯形算法的；當(dāng)要求 的計(jì)算精度更高的時(shí)候，Romberg算法是明顯快于復(fù)化梯形算法。本文所編寫(xiě)的程序適用于多項(xiàng)式的數(shù)值積分，且對(duì)于積分區(qū)間內(nèi)，被積函數(shù)

6、在每一點(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,f=poly2sym(a)/poly2sym(b) %利用梯形公式計(jì)算此數(shù)值積分disp('利用梯形公式計(jì)算數(shù)值積分的結(jié)果')kk=zeros();kk(1,1)=1/2*(M-m)/1*(subs(f,'x',m)+subs(f,'x',M) % for i=1:1:2A30t=0;for j=0:1:2A(i-1)-1 v=m+(2*j+1)*(M-m)/(2Ai)

7、 vv=polyval(a,v)/polyval(b,v);t=t+(M-m)/(2Ai)*vvendy=1/2*kk(i,1)+tkk(i+1,1)=y f=i+1;if(1/3*(kk(i+1,1)-kk(i,1)<=d) break;endend time=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，f=poly2sym(a)/poly2sym(b)disp('利用梯形公式計(jì)算數(shù)值

8、積分的結(jié)果')kk=zeros();kk(1,1)=1/2*(M-m)/1*(subs(f,'x',m)+subs(f,'x',M); % for i=1:1:2A40t=0;for j=0:1:2A(i-1)-1 v=m+(2*j+1)*(M-m)/(2Ai);vv=polyval(a,v)/polyval(b,v);積分下限m,積分上限M,以及計(jì)算精度d')%用于給用戶顯示被積函數(shù)的形式%用于存放結(jié)果先存儲(chǔ)首項(xiàng)%通項(xiàng)公式計(jì)算各項(xiàng)值%存儲(chǔ)其他項(xiàng)%記錄符合條件的值的下標(biāo)積分下限m,積分上限M,以及計(jì)算精度d')%用于給用戶顯示被積函數(shù)的形

9、式%用于存放結(jié)果先存儲(chǔ)首項(xiàng)endt=t+(M-m)/(29)*vv;y=1/2*kk(i,1)+t;kk(i+1,1)=y ;if(abs(1/3*(kk(i+1,1)-kk(i,1)<=d)disp('The result is:')kk()kk(i+1,1)break;elses=(4*kk(i+1,1)-kk(i,1)/(4-1);kk(i+1,2)=sif(i+1>=3)%通項(xiàng)公式計(jì)算各項(xiàng)值%存儲(chǔ)其他項(xiàng)%判斷梯形公式值是否達(dá)到要求%梯形值滿足要求，輸出結(jié)果%構(gòu)造 simpson 各項(xiàng)%存儲(chǔ)if(i+1>=3 & abs(1/15*(kk(i+1

10、,2)-kk(i,2)<=d) kk() disp('The result is:')%simpson 值滿足要求，輸出結(jié)果kk(i+1,2)pan1=0;break;elsec=(4A2*kk(i+1,2)-kk(i,2)/(4A2-1);% 構(gòu)造 cotes 值 kk(i+1,3)=c%存儲(chǔ) cotes 值if(i+1>=4)if(i+1>=4 & abs(1/63*(kk(i+1,3)-kk(i,3)<=d) disp('The result is:') kk(i+1,3) break;elser=(4A3*kk(i+1,3)-kk(i,3)/(4A3-1)% 構(gòu)造 romberg 值 kk(i+1,4)=r%存儲(chǔ) r