龍貝格積分的程序實現(共3頁)

1、精選優(yōu)質文檔-傾情為你奉上計算方法實驗報告3 【課題名稱】龍貝格積分的程序實現【目的和意義】函數變化有急有緩，為了照顧變化劇烈部分的誤差，我們需要加密格點。對于變化緩慢的部分，加密格點會造成計算的浪費。以此我們介紹一種算法，可以自動在變化劇烈的地方加密格點計算，而變化緩慢的地方，則取稀疏的格點。實際計算中，由于要事先給出一個合適的步長往往很困難，所以我們往往采用變步長的計算方案，即在步長逐步分半的過程中，反復利用復化求積公式進行計算，直到所求得的積分值滿足精度要求為止。在步長逐步分半過程中將粗糙的積分值逐步加工為精度較高的積分值，或者說將收斂緩慢的梯形值序列加工成收斂迅速的積分值序列。這種加速

2、方法稱為龍貝格算法?！居嬎愎健吭O表示復化梯形求得的積分值，其下標是等分數，由此則有遞推公式其中 ，其中由復化梯形公式的截斷誤差公式可得 ， 。由此可知， 。這樣導出的加速公式是辛普森公式：同理可得 。由此便可得加速的算法：龍貝格算法?！君堌惛袂蠓e算法流程圖】定義被積函數f，積分上下限a，b和精度c定義一個15×4的零矩陣，用于存放t值先按照公式算出t(1,1)至t(1,3)，以便循環(huán)按公式計算其余t值abs(t(k,4)-t(k-1,4)<c并且 k>6否是否是輸出“不收斂”輸出近似值k>=15【龍貝格求積算法Matlab主程序】functiont=rbg(f,a

3、,b,c) %定義龍貝格積分函數，f為待積函數，a與b為積分上下限，c為精度控制； t=zeros(15,4); %生成一零矩陣，用于存放t值； t(1,1)=(b-a)/2*(f(a)+f(b); %由于矩陣行列值均從1開始，所以將原本的t(0,0)記為t(1,1),行列均加1； for k=2:4 %先算出第一列的4個（包括t(1,1)）值，以便程后面可以直接循環(huán)計算； sum=0; for i=1:2(k-2) sum=sum+f(a+(2*i-1)*(b-a)/2(k-1); end t(k,1)=0.5*t(k-1,1)+(b-a)/2(k-1)*sum; for i=2:k t(k

4、,i)=(4(i-1)*t(k,i-1)-t(k-1,i-1)/(4(i-1)-1); end end for k=5:15 %循環(huán)按照公式計算出t值； sum=0; for i=1:2(k-2) sum=sum+f(a+(2*i-1)*(b-a)/2(k-1); end t(k,1)=0.5*t(k-1,1)+(b-a)/2(k-1)*sum; for i=2:4 t(k,i)=(4(i-1)*t(k,i-1)-t(k-1,i-1)/(4(i-1)-1); end if k>6 %可知最小二分次數，防止假收斂； if abs(t(k,4)-t(k-1,4)<c %若此時t值滿足精度，則輸出積分值； disp('答案 ',num2str(t(k,4); break; end end end if k>=15 disp('不收斂'); %二分次數達15次仍不收斂；end 【調用函數解題