數(shù)值積分算法與MATLAB實(shí)現(xiàn)陳悅

1、東北大學(xué)秦皇島分校數(shù)值計(jì)算課程設(shè)計(jì)報(bào)告數(shù)值積分算法及MATLAB實(shí)現(xiàn)學(xué) 院數(shù)學(xué)與統(tǒng)計(jì)學(xué)院專 業(yè)信息與計(jì)算科學(xué)學(xué) 號(hào)5133201姓 名陳悅指導(dǎo)教師姜玉山 張建波成 績(jī)教師評(píng)語：指導(dǎo)教師簽字： 2015年07月14日1 緒論數(shù)值分析是計(jì)算數(shù)學(xué)的一個(gè)主要部分，計(jì)算數(shù)學(xué)是數(shù)學(xué)科學(xué)的一個(gè)分支，它研究用計(jì)算機(jī)求解各種數(shù)學(xué)問題的數(shù)值檢索方其理論與軟件的實(shí)現(xiàn).而數(shù)值分析主要研究數(shù)值計(jì)算.現(xiàn)科學(xué)技術(shù)的發(fā)展與進(jìn)步提出了越來越多的復(fù)雜的數(shù)值計(jì)算問題，這些問題的圓滿解決已遠(yuǎn)人工手算所能勝任，必須依靠電子計(jì)算機(jī)快速準(zhǔn)確的數(shù)據(jù)處理能力.這種用計(jì)算機(jī)處理數(shù)值問題的方法，成為科學(xué)計(jì)算.今天，科學(xué)計(jì)算的應(yīng)用范圍非常廣泛，天氣

2、預(yù)報(bào)、工程設(shè)計(jì)、流體計(jì)算、經(jīng)濟(jì)規(guī)劃和預(yù)測(cè)以及國防尖端的一些科研項(xiàng)目，如核武器的研制、導(dǎo)彈和火箭的發(fā)射等，始終是科學(xué)計(jì)算最為活躍的領(lǐng)域.1.1 數(shù)值積分介紹數(shù)值積分是數(shù)值分析的重要環(huán)節(jié)，實(shí)際問題當(dāng)中常常需要計(jì)算積分，有些數(shù)值方法，如微分方程和積分方程的求解，也都和積分計(jì)算相聯(lián)系.求某函數(shù)的定積分時(shí)，在多數(shù)情況下，被積函數(shù)的原函數(shù)很難用初等函數(shù)表達(dá)出來，因此能夠借助微積分學(xué)的牛頓-萊布尼茲公式計(jì)算定積分的機(jī)會(huì)是不多的.另外，許多實(shí)際問題中的被積函數(shù)往往是列表函數(shù)或其他形式的非連續(xù)函數(shù)，對(duì)這類函數(shù)的定積分，也不能用不定積分方法求解.由于以上原因，數(shù)值積分的理論與方法一直是計(jì)算數(shù)學(xué)研究的基本課題.對(duì)微

3、積分學(xué)做出杰出貢獻(xiàn)的數(shù)學(xué)大師，如I.牛頓、L.歐拉、C.F.高斯、拉格朗日等人都在數(shù)值積分這個(gè)領(lǐng)域作出了各自的貢獻(xiàn)，并奠定了這個(gè)分支的理論基礎(chǔ).構(gòu)造數(shù)值積分公式最通常的方法是用積分區(qū)間上的n 次插值多項(xiàng)式代替被積函數(shù)，由此導(dǎo)出的求積公式稱為插值型求積公式.特別在節(jié)點(diǎn)分布等距的情形稱為牛頓-科特斯公式，例如梯形公式(Trapezoidal Approximations)與拋物線公式(Approximations Using Parabolas)就是最基本的近似公式.但它們的精度較差.龍貝格算法是在區(qū)間逐次分半過程中，對(duì)梯形公式的近似值進(jìn)行加權(quán)平均獲得準(zhǔn)確程度較高的積分近似值的一種方法，它具有公式

4、簡(jiǎn)練、計(jì)算結(jié)果準(zhǔn)確、使用方便、穩(wěn)定性好等優(yōu)點(diǎn)，因此在等距情形宜采用龍貝格求積公式(Rhomberg Integration).當(dāng)用不等距節(jié)點(diǎn)進(jìn)行計(jì)算時(shí)，常用高斯型求積公式計(jì)算，它在節(jié)點(diǎn)數(shù)目相同情況下，準(zhǔn)確程度較高，穩(wěn)定性好，而且還可以計(jì)算無窮積分.數(shù)值積分還是微分方程數(shù)值解法的重要依據(jù).許多重要公式都可以用數(shù)值積分方程導(dǎo)出.現(xiàn)探討數(shù)值積分算法以及運(yùn)用MATLAB軟件的具體實(shí)現(xiàn)1.2 MATLAB軟件MATLAB是美國MathWorks公司出品的商業(yè)數(shù)學(xué)軟件，用于算法開發(fā)、數(shù)據(jù)可視化、數(shù)據(jù)分析以及數(shù)值計(jì)算的高級(jí)技術(shù)計(jì)算語言和交互式環(huán)境，主要包括MATLAB和Simulink兩大部分.MATLAB

5、是matrix&laboratory兩個(gè)詞的組合，意為矩陣工廠（矩陣實(shí)驗(yàn)室）.是由美國mathworks公司發(fā)布的主要面對(duì)科學(xué)計(jì)算、可視化以及交互式程序設(shè)計(jì)的高科技計(jì)算環(huán)境.它將數(shù)值分析、矩陣計(jì)算、科學(xué)數(shù)據(jù)可視化以及非線性動(dòng)態(tài)系統(tǒng)的建模和仿真等諸多強(qiáng)大功能集成在一個(gè)易于使用的視窗環(huán)境中，為科學(xué)研究、工程設(shè)計(jì)以及必須進(jìn)行有效數(shù)值計(jì)算的眾多科學(xué)領(lǐng)域提供了一種全面的解決方案，并在很大程度上擺脫了傳統(tǒng)非交互式程序設(shè)計(jì)語言（如C、Fortran）的編輯模式，代表了當(dāng)今國際科學(xué)計(jì)算軟件的先進(jìn)水平.MATLAB和Mathematica、Maple并稱為三大數(shù)學(xué)軟件.它在數(shù)學(xué)類科技應(yīng)用軟件中在數(shù)值計(jì)算

6、方面首屈一指.MATLAB可以進(jìn)行矩陣運(yùn)算、繪制函數(shù)和數(shù)據(jù)、實(shí)現(xiàn)算法、創(chuàng)建用戶界面、連接其他編程語言的程序等，主要應(yīng)用于工程計(jì)算、控制設(shè)計(jì)、信號(hào)處理與通訊、圖像處理、信號(hào)檢測(cè)、金融建模設(shè)計(jì)與分析等領(lǐng)域.MATLAB的基本數(shù)據(jù)單位是矩陣，它的指令表達(dá)式與數(shù)學(xué)、工程中常用的形式十分相似，故用MATLAB來解算問題要比用C，F(xiàn)ORTRAN等語言完成相同的事情簡(jiǎn)捷得多，并且MATLAB也吸收了像Maple等軟件的優(yōu)點(diǎn)，使MATLAB成為一個(gè)強(qiáng)大的數(shù)學(xué)軟件.在新的版本中也加入了對(duì)C，F(xiàn)ORTRAN，C+，JAVA的支持.2 數(shù)值積分的基本概念 一般的，我們可以在區(qū)間上適當(dāng)選取某些節(jié)點(diǎn)，然后用的加權(quán)平均得

7、到平均高度的近似值，這樣構(gòu)造出的求積公式具有下列形式：，式中稱為求積節(jié)點(diǎn)；稱為求積系數(shù)，亦稱伴隨節(jié)點(diǎn)的權(quán).權(quán)僅僅與節(jié)點(diǎn)的選取有關(guān)，而不依賴于被積函數(shù)的具體形式.2.1 代數(shù)精度的概念如果某個(gè)求積公式對(duì)于次數(shù)不超過的多項(xiàng)式均能準(zhǔn)確的成立，但對(duì)于次多項(xiàng)式就不準(zhǔn)確成立，則稱該求積公式具有次代數(shù)精度（或代數(shù)精確值）一般地，欲使求積公式具有次代數(shù)精度，只要令它對(duì)于都能準(zhǔn)確成立，這就要求:2.2 求積公式的余項(xiàng)令求積公式的余項(xiàng)為，其中.；區(qū)間可以是有限的或無限的.構(gòu)造求積公式的問題就是確定和使得在某種意義下盡可能地小.3 數(shù)值積分方法及MATLAB實(shí)現(xiàn)3.1 復(fù)合辛普森公式 插值型求積公式用插值多項(xiàng)式替換

8、積分中的被積函數(shù)，然后計(jì)算作為積分的近似值，這樣建立的求積公式稱為插值型求積公式.用插值多項(xiàng)式的表達(dá)式，代入得，其中: .其余項(xiàng)為:.3.1.2 牛頓-科特斯公式介紹 取等距節(jié)點(diǎn)，把積分區(qū)間剖分成等分.令步長(zhǎng)，并記，則個(gè)節(jié)點(diǎn)為，代入得:=.這種等距節(jié)點(diǎn)的插值型求積公式通常稱為牛頓-科特斯公式. 辛普森公式利用牛頓-科特斯公式，取=2，此時(shí)為，即為辛普森公式，其余項(xiàng)為.3.1.4 復(fù)合辛普森公式將積分區(qū)間分成等分，分點(diǎn)為，其中，記區(qū)間的中點(diǎn)為，在每個(gè)小區(qū)間上用辛普森公式，則得到所謂的復(fù)合辛普森公式：.余項(xiàng)為，.3.1.5 復(fù)合辛普森公式的MATLAB實(shí)現(xiàn)代碼如下：function s=xinpu

9、sen(fun,a,b,n)h=(b-a)/ns1=0;s2=0;for k=0:(n-1) x=a+h*k; s1=s1+feval(fun,x);endfor k=0:(n-1) x=a+h*(k+1/2); s2=s2+feval(fun,x);ends=h/6*(feval(fun,a)+feval(fun,b)+2*s1+4*s2);3.2 龍貝格公式3.2.1 梯形法的遞推化將區(qū)間分成等份，共有個(gè)分點(diǎn)，如果將求積區(qū)間再二分一次，則分點(diǎn)增至個(gè)，用復(fù)合梯形公式求得該子區(qū)間上的積分值為，將每個(gè)子區(qū)間上的積分值相加得=.得到遞推公式：.3.2.2 龍貝格算法公式當(dāng)在上充分光滑時(shí)， 可證用逼

10、近的截?cái)嗾`差是：=按理查森外推法：，其中，為滿足的適當(dāng)正數(shù) .取序列：.用來逼近的誤差為，這種算法就是龍貝格算法.3.2.3 龍貝格算法MATLAB實(shí)現(xiàn)代碼如下：function s=longbeige(fun,a,b,tol)if nargin<4，tol=1e-4;endi=1;j=1;h=b-a;T(1,1)=h*(feval(fun,a)+feval(fun,b)/2;T(i+1,j)=T(I,j)/2+sum(feval(fun,a+h/2:h:b-h/2)*h/2;T(i+1,j+1)=(4j*T(i+1,j)-T(i.j)/(4j-1);while (abs(T(i+1,i

11、+1)-T(I,i)>tol) i=i+1;h=h/2; T(i+1,1)=T(I,1)/2+sum(feval(fun,a+h/2:h:b-h/2)*h/2; for j=1:i T(i+1,j+1)=(4j*T(i+1,j)-T(I,j)/(4j-1); endendTs=T(i+1,j+1);3.3 自適應(yīng)法自適應(yīng)積分法是一種比較經(jīng)濟(jì)而且快速的求積分的方法.他能自動(dòng)地在被積函數(shù)變化劇烈的區(qū)域增多節(jié)點(diǎn)，而在被積函數(shù)變化平緩的地方減少節(jié)點(diǎn).因此它是一種不均勻區(qū)間的積分方法.按照子區(qū)間上的積分方式它可以分為自適應(yīng)辛普森積分法和自適應(yīng)梯形積分法.通常是采用自適應(yīng)辛普森積分法作為子區(qū)間的積分

12、方式. 自適應(yīng)積分法的基本步驟如下：(1) 將積分區(qū)間分成兩個(gè)相等的1級(jí)子區(qū)間和，且；(2) 在上述兩個(gè)1級(jí)子區(qū)間上用辛普森積分得到積分和;(3) 將子區(qū)間分成兩個(gè)相等的2級(jí)子區(qū)間和;(4) ;(5) 比較和， 如果| - |<，其中為整體積分所需要的精度，則認(rèn)為子區(qū)間上的積分已達(dá)到所需精度，不需要再細(xì)分;否則就需要再細(xì)分，對(duì)每個(gè)2級(jí)子區(qū)間做同樣的判斷1級(jí)子區(qū)間的操作過程完全與上面相同.3.4 高斯法對(duì)已知求積公式可以討論它的代數(shù)精度，反之也可以按照代數(shù)精度要求導(dǎo)出求積公式.對(duì)于求積公式，當(dāng)求積節(jié)點(diǎn)（）固定時(shí)，公式有個(gè)待定參數(shù)，故此時(shí)可要求滿足對(duì)“準(zhǔn)確”這樣個(gè)約束條件，從而使之至少具有次

13、代數(shù)精度.進(jìn)一步，可考慮將也視為待定參數(shù)，這樣公式的待定參數(shù)就有個(gè)，從而可望公式的代數(shù)精度達(dá)到.此類高精度的求積公式稱為高斯型公式，而對(duì)應(yīng)的節(jié)點(diǎn)稱為區(qū)間上的高斯點(diǎn).3.4.1 高斯法的MATLAB實(shí)現(xiàn)代碼如下：function g=gaosi(fname,a,b,n,m)switch m case 1 t=0;A=1; case 2 t=-1/sqrt(3),1/sqrt(3);A=1,1; case 3 t=-sqrt(0.6),0.0,sqrt(0.6);A=5/9,8/9,5/9; case 4 t=-0.8612136,-0.339981,0.339981,0.861136; A=0.

14、347855,0.652145,0.652145,0.347855; case 5 t=-0.906180,-0.538469,0.0,0.538469,0.906180; A=0.236927,0.478629,0.568889,0.478629,0.236927; case 6 t=-0.932470,-0.661209,-0.238619,0.238619,0.661209,0.932470; A=0.171325,0.360762,0.467914,0.467914,0.360762,0.171325; otherwise errorendx=linspace(a,b,n+1);g=0

15、;for i=1:n g=g+gsint(fname,x(i),x(i+1),A,t);endfunction g=gsint(fname,a,b,A,t)g=(b-a)/2*sum(A.*feval(fname,(b-a)/2*t+(a+b)/2);3.5 多重積分法考慮二重積分，它是曲面與平面區(qū)域圍成的體積，對(duì)于矩形區(qū)域.可將它寫成累次積分若用復(fù)合辛普森公式，可分別將分成N，M等份，步長(zhǎng)，先對(duì)積分應(yīng)用復(fù)合辛普森公式，令，.從而得：.4 數(shù)值積分方法比較例：用數(shù)值積分方法計(jì)算方程的值.4.1 復(fù)合辛普森公式求解>> fun=inline('4./(1+x.2)')

16、;>> xinpusen(fun,0,1,10)ans = 3.2749259863031184.2 龍貝格算法求解longbeige(inline('4./(1+x.2)'),0,1,1e-6)T = 3.0000 3.1000 3.1333 3.1312 3.1416 3.1421 3.1390 3.1416 3.1416 3.1416 3.1409 3.1416 3.1416 3.1416 3.1416 3.1414 3.1416 3.1416 3.1416 3.1416 3.1416ans =4.3 高斯算法求解>> gaosi(inline(&

17、#39;4./(1+x.2)'),0,1,2,3)ans =>> gaosi(inline('4./(1+x.2)'),0,1,4,4)ans =4.4 三種方法比較分析結(jié)果顯示每一個(gè)算法都接近真實(shí)值，但龍貝格算法相比較復(fù)合辛普森算法，高斯算法來說更加接近.對(duì)于代數(shù)精度來說，復(fù)合辛普森的代數(shù)精度為11，龍貝格代數(shù)精度為11，高斯代數(shù)精度為11.可見代數(shù)精度相同時(shí)，龍貝格的求積精度最小，所以相同條件下龍貝格求積公式最能接近準(zhǔn)確值.總結(jié)隨著數(shù)學(xué)實(shí)驗(yàn)的興起， 對(duì)整個(gè)數(shù)學(xué)課程教學(xué)改革起到了積極的推動(dòng)作用， 我們要熟悉的運(yùn)用各種數(shù)學(xué)軟件，解決數(shù)學(xué)運(yùn)算中繁瑣的問題，實(shí)現(xiàn)學(xué)習(xí)的簡(jiǎn)單，快捷化.同時(shí)意識(shí)到用MATLAB編程時(shí)，要實(shí)現(xiàn)代碼的層次性，做到有規(guī)有矩，那樣才能把MATLAB運(yùn)用自如.這次課程設(shè)計(jì)，用MATLAB實(shí)驗(yàn)對(duì)數(shù)值積分進(jìn)行了實(shí)現(xiàn)，簡(jiǎn)介了5種不同的數(shù)值積分的方法，并且實(shí)現(xiàn)了其中的3中方法，實(shí)現(xiàn)過程中發(fā)現(xiàn)了各種方法之間的區(qū)別和