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

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

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

3、做出杰出貢獻(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)出的求積公式稱(chēng)為插值型求積公式特別在節(jié)點(diǎn)分布等距的情形稱(chēng)為牛頓-科特斯 公式，例如梯形公式 (Trapezoidal Approximate ns )與拋物線(xiàn)公式(Approximate nsUsing Parabolas ）就是最基本的近似公式.但它們的精度較差.龍貝格算法是在區(qū)間逐次 分半過(guò)程中，對(duì)梯形公式的近似值進(jìn)行加權(quán)平均獲得準(zhǔn)確程度較高的積分近似值的一種 方法，它具有公式簡(jiǎn)練

4、、計(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ì)算無(wú)窮積分.數(shù)值積分還是微分方程數(shù)值解法的重要依據(jù).許多重要公式都可以用數(shù)值 積分方程導(dǎo)出.現(xiàn)探討數(shù)值積分算法以及運(yùn)用 MATLAB軟件的具體實(shí)現(xiàn)1.2 MATLAB 軟件MATLAB是美國(guó)MathWorks 公司出品的商業(yè)數(shù)學(xué)軟件，用于算法開(kāi)發(fā)、數(shù)據(jù)可 視化、數(shù)據(jù)分析以及數(shù)值計(jì)算的高級(jí)技術(shù)計(jì)算語(yǔ)言和交互式環(huán)境，主要包括MATLAB和Simulink兩大部分

5、.MATLAB是matrix&laboratory兩個(gè)詞的組合，意為矩陣工廠(chǎng)（矩陣實(shí)驗(yàn)室）.是由美國(guó)mathworks公司發(fā)布的主要面對(duì)科學(xué)計(jì)算、可視化以及交互式程序設(shè)計(jì)的高科 技計(jì)算環(huán)境.它將數(shù)值分析、矩陣計(jì)算、科學(xué)數(shù)據(jù)可視化以及非線(xiàn)性動(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ì)語(yǔ)言（如 C、Fortran ）的編輯模式，代表了當(dāng)今國(guó)際科學(xué)計(jì)算 軟件的先進(jìn)水平.MATLAB和Mathematica、Maple并稱(chēng)為三大數(shù)學(xué)軟件.它在數(shù)學(xué)

6、類(lèi)科技應(yīng)用軟 件中在數(shù)值計(jì)算方面首屈一指.MATLAB可以進(jìn)行矩陣運(yùn)算、繪制函數(shù)和數(shù)據(jù)、實(shí)現(xiàn)算 法、創(chuàng)建用戶(hù)界面、連接其他編程語(yǔ)言的程序等，主要應(yīng)用于工程計(jì)算、控制設(shè)計(jì)、信 號(hào)處理與通訊、圖像處理、信號(hào)檢測(cè)、金融建模設(shè)計(jì)與分析等領(lǐng)域.MATLAB的基本數(shù)據(jù)單位是矩陣，它的指令表達(dá)式與數(shù)學(xué)、工程中常用的形式十 分相似，故用MATLAB來(lái)解算問(wèn)題要比用C，F(xiàn)ORTRAN等語(yǔ)言完成相同的事情簡(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ū)間a,b

7、上適當(dāng)選取某些節(jié)點(diǎn)Xk，然后用f Xk的加權(quán)平均得到平均高度f(wàn) 的近似值，這樣構(gòu)造出的求積公式具有下列形式：anb f x dxAk f Xk ，式中Xk稱(chēng)為求積節(jié)點(diǎn)；Ak稱(chēng)為求積系數(shù)，亦稱(chēng)伴隨節(jié)點(diǎn)Xk的k 0權(quán).權(quán)Ak僅僅與節(jié)點(diǎn)xk的選取有關(guān)，而不依賴(lài)于被積函數(shù) f x的具體形式.2.1代數(shù)精度的概念如果某個(gè)求積公式對(duì)于次數(shù)不超過(guò) m的多項(xiàng)式均能準(zhǔn)確的成立，但對(duì)于 m 1次多項(xiàng)式就不準(zhǔn)確成立，則稱(chēng)該求積公式具有 m次代數(shù)精度（或代數(shù)精確值）一般地，欲使求積公式具有 m次代數(shù)精度，只要令它對(duì)于 f x 1,x,L ,xm都能準(zhǔn)確成立，這就要求：Ak=b-a,. 1,22AkXk=b -a2M1

8、m+1m+1m+1b -a令求積公式的余項(xiàng)為bR f，其中 R f = a2.2 求積公式的余項(xiàng)nX dxAkf Xk .；區(qū)間a,b可以是k 0有限的或無(wú)限的.構(gòu)造求積公式的問(wèn)題就是確定 Xj和Aj使得R f在某種意義下盡可能 地小.3 數(shù)值積分方法及MATLAB實(shí)現(xiàn)3.1復(fù)合辛普森公式3.1.1 插值型求積公式用插值多項(xiàng)式Ln(X)替換積分I*x dx中的被積函數(shù)f然后計(jì)算bILnax dx作為積分的近似值，這樣建立的求積公式稱(chēng)為插值型求積公式用插值多n項(xiàng)式Ln(x)的表達(dá)式Ln(x)lk Xk 0X k ，代入得InAk f Xk，其中:k 0bAk(X Xo)L (x Xk i)(x

9、Xk i)L (x Xn)a (Xk X°)L (Xk Xk 1)(Xk Xk1)Ldx .(Xk Xn)(n1L XdX.1)!取等距節(jié)點(diǎn)，把積分區(qū)間a, b剖分成n等分.令步長(zhǎng)h j，并記X0na,Xn b，則n 1個(gè)節(jié)點(diǎn)為Xk X0 kh,k0,1L ,n，代入得:(1)n khn(n k)!k! 0t(t 1)L (t k 1)(t k1)L(t n)dt.*b其余項(xiàng)為:R f I Ia (n3.1.2 牛頓-科特斯公式介紹這種等距節(jié)點(diǎn)的插值型求積公式通常稱(chēng)為牛頓 -科特斯公式.3.1.3 辛普森公式利用牛頓-科特斯公式,b a取n=2，此時(shí)為S w f(a)4f(節(jié))f(b

10、)，即為辛普森公式，其余項(xiàng)為Rs f*(b a)5f(4)().9023.1.4復(fù)合辛普森公式將積分區(qū)間a,b分成n等分，分點(diǎn)為Xk a kh(k 0,1,L , n),區(qū)間Xk,Xk 1的中點(diǎn)為X 1，在每個(gè)小區(qū)間Xk12Xk , Xk 1上用辛普森公式,其中h 引，記n貝U得到所謂的復(fù)合辛普森公式:Snn1hk 0 6xk4fxk22xk 1Xkxk22余項(xiàng)為R Sn18刖廠(chǎng)(a,b).3.1.5 復(fù)合辛普森公式的MATLAB實(shí)現(xiàn)代碼如下:function s=x in puse n(fun, a,b ,n) h=(b-a)/n s1=0;s2=0; for k=0:( n-1)x=a+h

11、*k;s1=s1+feval(f un, x);end for k=0:( n-1)x=a+h*(k+1/2); s2=s2+feval(fu n,x);ends=h/6*(feval(fu n, a)+feval(fu n,b)+2*s1+4*s2);3.2 龍貝格公式3.2.1 梯形法的遞推化將區(qū)間a,b分成n等份，共有n 1個(gè)分點(diǎn)，如果將求積區(qū)間再二分一次，則分點(diǎn)增Xk2f Xk 1f(Xk1)，2至2n 1個(gè)，用復(fù)合梯形公式求得該子區(qū)間上的積分值為-4將每個(gè)子區(qū)間上的積分值相加得T2n = -f Xkf ( Xk i)n 1hkof(Xk1).21得到遞推公式：T2n-Tnxk丄23.

12、2.2 龍貝格算法公式當(dāng)f x在a,b上充分光滑時(shí),可證用(h)逼近I的截?cái)嗾`差是:I T|(h)= a/2 a2h4 L a2k L按理查森外推法:R(h)F h，F(xiàn)m1(h) Fm(qh) qpmFm(h),(m 1,2,L )Pm其中,q為滿(mǎn)足1 qPm用Tm11 q0(m 1,2,L )的適當(dāng)正數(shù).取序列：mh4 TmTn hAm, m 1,2,L.41，這種算法就是龍貝格算法.Tm1 hh來(lái)逼近I的誤差為0 h2 m 1龍貝格算法MATLAB實(shí)現(xiàn)代碼如下:function s=lon gbeige(fu n, a,b,tol) if nargin<4 ，tol=1e-4;end

13、i=1;j=1;h=b-a;T(1,1)=h*(feval(fu n,a)+feval(fu n,b)/2;T(i+1,j)=T(l,j)/2+sum(feval(fu n,a+h/2:h:b-h/2)*h/2;T(i+1,j+1)=(4Aj*T(i+1,j)-T(i.j)/(4Aj-1);while (abs(T(i+1,i+1)-T(l,i)>tol)i=i+1;h=h/2;T(i+1,1)=T(I,1)/2+sum(feval(fu n,a+h/2:h:b-h/2)*h/2;for j=1:iT(i+1,j+1)=(4Aj*T(i+1,j)-T(I,j)/(4Aj-1);enden

14、dTs=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ū)間的積分方式.自適應(yīng)積分法的基本步驟如下：11(1) 將積分區(qū)間a,b分成兩個(gè)相等的1級(jí)子區(qū)間a,a -h和a丄h,a h，且22h b a ；1(2) 在上述兩個(gè) 1級(jí)子區(qū)間上用辛普森積分得到積分I a,a丄h八1和a -h,a h A 12將子區(qū)間分成兩個(gè)相等的2級(jí)子區(qū)間 a,

15、 a2丄h,a22h ；1a, a h2a, a11 2a h, a2比較I _ a,a1h a2，如果 | I _ a, a 1h a 1a, a 丄 h a 22|< 10*1gpsi，其中epsi為整體積分所需要的精度，則認(rèn)為子區(qū)1間a，a 2h上的積分已達(dá)到所需精度，不需要再細(xì)分；否則就需要代碼如下：fun ctio n g=gaosi(f name,a,b ,n,m)switch mcase 1t=0;A=1;case 2t=-1/sqrt(3),1/sqrt(3);A=1,1;case 3t=-sqrt(0.6),0.0,sqrt(0.6);A=5/9,8/9,5/9;cas

16、e 4t=-0.8612136,-0.339981,0.339981,0.861136;A=0.347855,0.652145,0.652145,0.347855;case 5t=-0.906180,-0.538469,0.0,0.538469,0.906180;A=0.236927,0.478629,0.568889,0.478629,0.236927;case 6t=-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;

17、 otherwiseerrorendx=li nspace(a,b, n+1);g=o；for i=1: ng=g+gsi nt(fname,x(i),x(i+1),A,t);end function g=gsint(fname,a,b,A,t)g=(b-a)/2*sum(A.*feval(fname,(b-a)/2*t+(a+b)/2);3.5 多重積分法考慮二重積分f(x,y)dA，它是曲面z f (x,y)與平面區(qū)域R圍成的體積，對(duì)于矩R形區(qū)域 R (x,y) a x b,c y d .可將它寫(xiě)成累次積分f (x, y)dxRba(d將a,b,c,d分成N , M等份，步長(zhǎng)a,kd cd

18、百，先對(duì)積分c f(X,y)dy應(yīng)用復(fù)合辛普森公式，令yc ik， y 1 ci -2(i扣.dc f(x,y)dyk f (x,yo)6M 14f(x,yii oii)2M 12 f (x, yi)f (x, yM ).i 1f (x, y)dy)dx若用復(fù)合辛普森公式，可分別 c從而得:kf (x, y)dydx -6ba f(x,yo)M 1b4 af(x,y.1)ai i 0i 2bba f (x, yi)af(x,yM)dx .aa4數(shù)值積分方法比較例：用數(shù)值積分方法計(jì)算方程4.1復(fù)合辛普森公式求解 >> fun=i nli ne('4丿(l+x.T)')

19、;>> xin puse n(fun, 0,1,10) ans =3.2749259863031184.2 龍貝格算法求解Ion gbeige(i nlin e('4./(1+x.A2)'),0,1,1e-6) ans =3.10003.13333.13123.14163.14213.13903.14163.14163.14163.14093.14163.14163.14163.14163.14143.14163.14163.14163.14163.00003.14164.3高斯算法求解>> gaosi(i nlin e('4./(1+x.A2)

20、'),0,1,2,3)ans =>> gaosi(i nlin e('4./(1+x.A2)'),0,1,4,4)ans =4.4三種方法比較分析結(jié)果顯示每一個(gè)算法都接近真實(shí)值，但龍貝格算法相比較復(fù)合辛普森算法，高斯算法來(lái)說(shuō)更加接近.對(duì)于代數(shù)精度來(lái)說(shuō)，復(fù)合辛普森的代數(shù)精度為11，龍貝格代數(shù)精度為11，高斯代數(shù)精度為11.可見(jiàn)代數(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)算中繁瑣的問(wèn)題，實(shí)現(xiàn)學(xué)習(xí)的簡(jiǎn)單，快捷化.同時(shí) 意識(shí)到用MATLAB編程時(shí)，要實(shí)現(xiàn)代碼的層次性，做到有規(guī)有矩，那樣才能把 MATLAB 運(yùn)用自如.