1、海南大學(xué)數(shù)理統(tǒng)計課程設(shè)計題目：基于蒙特卡洛思想的定積分?jǐn)?shù)值解法班級：信息與計算科學(xué)姓名：體貼的瑾色學(xué)號：指導(dǎo)教師：日期：2017.06目錄基于蒙特卡洛方法的定積分?jǐn)?shù)值解法3摘要3Abstract3一、前言4二、蒙特卡洛求解定積分法42.1 隨機(jī)投點法42.2 平均值法62.3 兩種方法的比較7三、蒙特卡洛計算二重積分8四、結(jié)語10參考文獻(xiàn):10附錄：MATLAB®序10基于蒙特卡洛方法的定積分?jǐn)?shù)值解法摘要微積分是現(xiàn)代數(shù)學(xué)和現(xiàn)代物理的一個重要基礎(chǔ)，而定積分在生活生產(chǎn)上也具有十分廣泛的應(yīng)用，然而，定積分的計算特別是涉及到較復(fù)雜的定積分的話便變的十分困難.本文介紹了蒙特卡洛方法解決定積分求

2、解的兩種思路：隨機(jī)取點法和平均值法，并給出了用matlab實現(xiàn)兩種方法的算法程序，最后用兩種方法分別計算了幾個實例.關(guān)鍵詞：蒙特卡洛定積分matlabAbstractCalculusisanimportantfoundationofmodernmathematicsandmodernphysics,anddefiniteintegralarealsowidelyusedinlifeproduction.However,itisverydifficulttocalculatethedefiniteintegral,especiallywhenitcomestomorecomplex.Thispa

3、perintroducestwomethodsofMonteCarlomethodtosolvedefiniteintegralsolution:randompointmethodandmeanmethod,andgivesthealgorithmtorealizetwomethodsbymatlab.Finally,weusetwomethodstocalculateseveralexamplesrespectively.Keys:MonteCarlodefiniteintegralmatlab一、前言蒙特卡洛方法(英語：MonteCarlomethod),也稱統(tǒng)計模擬方法，是1940年代中

4、期由于科學(xué)技術(shù)的發(fā)展和電子計算機(jī)的發(fā)明，而提出的一種以概率統(tǒng)計理論為指導(dǎo)的數(shù)值計算方法.是指使用隨機(jī)數(shù)(或更常見的偽隨機(jī)數(shù))來解決很多計算問題的方法.蒙特卡洛方法是屬于統(tǒng)計實驗的一種方法，主要通過統(tǒng)計抽樣實驗為各種各樣的數(shù)學(xué)問題提供近似解，所以也被稱為隨機(jī)抽樣技術(shù).而正是因為蒙特卡洛方法的原理使得這種方法得到的結(jié)果基本上都會存在誤差，不過可喜的是這種誤差隨著取的隨機(jī)數(shù)的數(shù)量的增大而減小.所以，隨著計算機(jī)的發(fā)展，蒙特卡洛算法的精確度一直在提高.定積分的求解是蒙特卡洛方法解決的最好的問題之一，特別是在涉及到重積分的情況下，這種方法往往是科研工作者比較喜歡的方法之一.二、蒙特卡洛求解定積分法2.1隨

5、機(jī)投點法s伯努利大數(shù)定理說明：隨著試驗次數(shù)n的增大，事件A發(fā)生的頻率Sn與其頻率p的偏n差Sn-p大于預(yù)定給定的精度名的可能性愈來愈小，要多小就多小，這就是說當(dāng)試驗次數(shù)n足夠多的話，頻率穩(wěn)定于概率.所以，隨機(jī)投點法的方法的原理就是重復(fù)進(jìn)行大量實驗，然i后用觀測到的頻率代替概率作為所求結(jié)果.比如求定積分J=(f(x)dx，其中0Ef(x)E1,可設(shè)二維隨機(jī)變量(X,Y)服從正方形0<x<1,0<y<1上的均勻分布，則可知X服從0,1上的均勻分布，Y也服從0,1上的均勻分布且X和Y相互獨立，又記事件A=丫三f(X),則A的概率為1 f(x)1p=P(Y_f(X)=00dyd

6、x=0f(x)dx=J我們只需要找到A事件出現(xiàn)的頻率即可，即將(X丫)看成是向正方形0WxW1,0WyW"內(nèi)的隨機(jī)投的點，用隨機(jī)點在區(qū)域y<f(x)中的頻率作為定積分的值.算法流程是:(1) .先用計算機(jī)產(chǎn)生(0,1上均勻分布的2n個隨機(jī)數(shù)，組成n對隨機(jī)數(shù)(x,yi),i=1,2,.n,其中n應(yīng)該足夠大(2) .對這n對隨機(jī)數(shù)(x,y),i=1,2,.n,記錄滿足yWf(Xj)的次數(shù)，由此可得到事件SSA發(fā)生的頻率,則J=五.nn,b注意對于一般區(qū)間a,b上的定積分J=fg(x)dx,其中g(shù)(x)為可積函數(shù)，文獻(xiàn)1,a給出一種解法是做線性變換將J轉(zhuǎn)換成J的形式.本文給出另外一種

7、算法來計算J,我們知道一重定積分的幾何意義是求曲線和坐標(biāo)軸所包圍的面積，但是這里面臨一個問題就是包圍的面積可負(fù)可正，另外此時的總區(qū)域面積也不再是1,所以如果按照上述流程(2)顯然會出現(xiàn)問題，所以，本文對(2)進(jìn)行改良后的新的第二步是：(2)記c=max(f(x),d=mm(f(x),記錄滿足0Wy<f(x)出現(xiàn)的次數(shù)k1,和滿kko一足f(x)WyM0的次數(shù)k2,算得新的頻率為.右CA0且d<0,則nS=(cd)*(b-a),或者dA0則S=c*(ba)或者cc0,則S=d*(ba).最后得至Uj=knk*s.n2 金12譬如計算J=Je/，2ndx，其精確值由matlab函數(shù)in

8、tegral給出(下同)，取n=10，,運(yùn)行100次后得到的結(jié)果為口用機(jī)超與法巧用辟。泯口魄Q甌D吃DE4EDl=4i100模擬后得到的值的分析為(誤差計算公式為工(模擬值M確值)2，下同)：t曾方法平均值力差誤差精確解隨機(jī)投點法0.650153.733E-043.830E-040.64984當(dāng)n=1000時的投點圖為2.2平均值法1計算定積分J=If(x)dx,其中f(x)可積.設(shè)隨機(jī)變量X服從(0,1摩的均勻分布，則Y=f(x)的數(shù)學(xué)期望為iE(f(x)=J。f(x)dx=J.所以估計期望就是估計所求值，由辛欽大數(shù)定理，可以用f(X)觀察值的平均去估計f(X)的期望值.具體做法如下：先用計

9、算機(jī)產(chǎn)生n個在(0,1)上均勻分布的隨機(jī)數(shù)x，i=1,2,.n,然后對每個x計算f(X),最后得到J的估計值為1 nbJ=Zf(xj.至于一般區(qū)間(a,b)上定積分J=1f(x)dx的計算，按照上述流程因為ni1ab11b-anE(f(x)=-f(x)dxJ,所以Jf(x).ab-ab-ani1譬如計算J2x/2dx6取n=10,運(yùn)行100次后得到的結(jié)果為835精而貂+攫報也1m1p1口通加4U如E07的nIUJ模擬后得到的值的分析為方法平均值力差誤差精確解平均值法1.912942.378E-032.381E-031.913122.3兩種方法的比較本為將通過計算幾個例子來比較兩種方法的優(yōu)劣性比

10、如計算J=ed2x/、Wdx2前林方法的比和結(jié)果I1I-看存價十平歸值去018.JS16.41B3,B1C.3T8219.15T611805C1口3030如505口m90制HDD次數(shù)綜合的的比較為：方法平均值、.、.廣.力左精確解隨機(jī)投點法18.250364.732E-014.732E-0118.25077平均值法18.253138.831E-028.886E-0218.25077又比如計算J=yX2-xdx結(jié)果為：綜合比較結(jié)果為:方法平均值、.、.廣.力左精確解隨機(jī)投點法9.164791.898E-011.901E-019.16667平均值法9.165898.122E-028.128E-02

11、9.16667二匚*1C20前叩國溝703：90兩種方法的匕軸造果顯然在一重定積分的計算上，雖然兩種方法平均值都接近于精確解，但是不管是平均值法方差還是誤差都較之隨機(jī)投點法更小，所以平均值法性能更好一此三、蒙特卡洛計算二重積分計算二重積分其實和計算一重積分差距并不是很大。在隨機(jī)投點法中與只需要將生成的隨機(jī)數(shù)變成三維便，y,z)并把投點區(qū)域換成一個正方體即可，比如求區(qū)間在ki-k?(a,b)*(e,f)上的二重積分，計算公式應(yīng)該為J一2*S*(f-e)其中k1,k2,n,S與n上文含義一致.而平均值法除了生成的隨機(jī)數(shù)也要變成三維(x,yz)之外，最大的不同在于計算一般積分區(qū)域(a,b)*(e,f

12、)上的積分，計算公式應(yīng)該變成j=(b-a)*(f-e)£f(為).另外值得一提的是，因為matlab在最優(yōu)化上效果并不是十nii分的出色，所以本文在隨機(jī)投點法中為了求投點區(qū)域，本文采用lingo求函數(shù)最大值和最小值.譬如求2X-3.5110210y5dxdy取n=106,使用兩種方法得到的結(jié)果為黃種吃小一叟.士勉令總工S2W£O70CDydxdy"他項設(shè)下上情扁M平均F徜方法平均值力差誤差精確解隨機(jī)投點法5.203144.70790E-044.71132E-045.20320平均值法5.203212.14826E-052.15035E-055.20320又比如求積

13、分兩種方法計用二歪雙分的轉(zhuǎn)果圖-0.&270.626Q蹈-0.S3.的0&32-0&33-0.63-D&35OCTcr+*D10劉3D409Q0D7DBD90100模擬后得到的結(jié)果為:O隨機(jī)投點法精踴解平均值法方法平均值力差誤差精確解隨機(jī)投點法5.208363.64836E-013.67495E-015.20320平均值法5.202682.12750E-022.13024E-025.20320-&3E次封四、結(jié)語通過上文給出的結(jié)果，可以看到蒙特卡洛方法用來求定積分的誤差十分的小，每次得到的結(jié)果也十分的穩(wěn)定，所以是求解定積分的一種很好的方法.但是，顯然不管

14、在計算一重還是二重定積分上，平均值法還是比隨機(jī)投點法性能要更好一些參考文獻(xiàn)：1其詩松，程依明.概率論與數(shù)理統(tǒng)計M.高等教育出版社.20112朱陸陸.蒙特卡洛方法及應(yīng)用D.華中師范大學(xué),2014.附錄：MATLABfunctiongailv(time,num,dlim,ulim)%time次數(shù)%num生成的點數(shù)%min最大值%max最小值%隨機(jī)投點法fun=(x)exp(x.A2)-x;gun=(x)-(exp(x.A2)-x);val=integral(fun,dlim,ulim)I=zeros(time,1);,min=fminbnd(fun,dlim,ulim);,max=fminbnd(

15、gun,dlim,ulim);max=-max;fori=1:timei1=0;i2=0;xi=unifrnd(dlim,ulim,num,1);ifmax<0yi=min+(-min)*rand(num,1);i2=1;elseifmin>0yi=max*rand(num,1);i1=1;elseyi=min+(max-min)*rand(num,1);endendy=fun(xi);k=sum(yi<=y&yi>0&y>0);p=sum(yi>=y&yi<0&y<0);ifi=100a=(yi<=y&am

16、p;yi>0&y>0);b=(yi>=y&yi<0&y<0);c=(a+b);holdonplot(xi(a),yi(a),'ro')plot(xi(b),yi(b),'bo')plot(xi(c),yi(c),'co')%fplot('(exp(-x.A3)-x)',dlim,ulim)holdoffendI(i,1)=(k-p)/num*(max-min)+i1*min-i2*max)*(ulim-dlim)*(ulim-dlim);ave_r=sum(I)/time;wu

17、cha_r=sum(I-val).A2);var_r=sum(I-ave_r).A2);endfigureplot(1:time,I,'bo');%¥均值法fori=1:timexi=dlim+(ulim-dlim)*rand(num,1);y=fun(xi);I(i,1)=sum(y)/num;endholdonplot(1:time,zeros(1,time)+val,'c')title(,隨機(jī)投點法計算定積分結(jié)果,);plot(1:time,I,'ro');holdoffave_a=sum(I)/time;var_a=sum(I-

18、ave_a).A2);wucha_a=sum(I-val).A2);fprintf('方法ttt平均值ttt方差tttt誤差n')fprintf(,隨機(jī)投點法t%1.5ftt%1.5Ett%1.5En',ave_r,var_rwucha_r);fprintf('平均值法tt%1.5ftt%1.5Ett%1.5En',ave_a,var_awucha_a);二重定積分程序：functiongailv2(time,num,dlim,ulim,ydlim,yulim)%time次數(shù)%num生成的點數(shù)%min最大值%max最小值%隨機(jī)投點法fun=(x,y)(-

19、x.A2+y+5);val=dblquad(fun,dlim,ulim,ydlim,yulim);I=zeros(time,1);max=6.3;min=4.06;fori=1:timei1=0;i2=0;xi=dlim+(ulim-dlim)*rand(num,1);yi=ydlim+(yulim-ydlim)*rand(num,1);ifmax<0zi=min+(-min)*rand(num,1);i2=1;elseifmin>0zi=max*rand(num,1);i1=1;elsezi=min+(max-min)*rand(num,1);endendy=fun(xi,yi);k=sum(zi<=y&zi>0&y>=0);p=sum(zi>=y&zi<0&y<0);plot3(xi,yi,y,'b');holdoffendI(i,1)=(k-p)/num*(max-min)+i1*min-i2*max)*(ulim-dlim)*(yulim-ydlim);endave_r=sum(I)/time;wucha_r=sum(I-val).A2);var_r=sum(I-ave_r).A2);