基于蒙特卡洛思想的定積分數(shù)值解法

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

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

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

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

5、機投點法s伯努利大數(shù)定理說明：隨著試驗次數(shù)n的增大，事件A發(fā)生的頻率Sn與其頻率p的偏n差Sn-p大于預定給定的精度名的可能性愈來愈小，要多小就多小，這就是說當試驗次數(shù)n足夠多的話，頻率穩(wěn)定于概率.所以，隨機投點法的方法的原理就是重復進行大量實驗，然i后用觀測到的頻率代替概率作為所求結果.比如求定積分J=(f(x)dx，其中0Ef(x)E1,可設二維隨機變量(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"內的隨機投的點，用隨機點在區(qū)域y<f(x)中的頻率作為定積分的值.算法流程是:(1) .先用計算機產生(0,1上均勻分布的2n個隨機數(shù)，組成n對隨機數(shù)(x,yi),i=1,2,.n,其中n應該足夠大(2) .對這n對隨機數(shù)(x,y),i=1,2,.n,記錄滿足yWf(Xj)的次數(shù)，由此可得到事件SSA發(fā)生的頻率,則J=五.nn,b注意對于一般區(qū)間a,b上的定積分J=fg(x)dx,其中g(x)為可積函數(shù)，文獻1,a給出一種解法是做線性變換將J轉換成J的形式.本文給出另外一種

7、算法來計算J,我們知道一重定積分的幾何意義是求曲線和坐標軸所包圍的面積，但是這里面臨一個問題就是包圍的面積可負可正，另外此時的總區(qū)域面積也不再是1,所以如果按照上述流程(2)顯然會出現(xiàn)問題，所以，本文對(2)進行改良后的新的第二步是：(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，,運行100次后得到的結果為口用機超與法巧用辟。泯口魄Q甌D吃DE4EDl=4i100模擬后得到的值的分析為(誤差計算公式為工(模擬值M確值)2，下同)：t曾方法平均值力差誤差精確解隨機投點法0.650153.733E-043.830E-040.64984當n=1000時的投點圖為2.2平均值法1計算定積分J=If(x)dx,其中f(x)可積.設隨機變量X服從(0,1摩的均勻分布，則Y=f(x)的數(shù)學期望為iE(f(x)=J。f(x)dx=J.所以估計期望就是估計所求值，由辛欽大數(shù)定理，可以用f(X)觀察值的平均去估計f(X)的期望值.具體做法如下：先用計

9、算機產生n個在(0,1)上均勻分布的隨機數(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,運行100次后得到的結果為835精而貂+攫報也1m1p1口通加4U如E07的nIUJ模擬后得到的值的分析為方法平均值力差誤差精確解平均值法1.912942.378E-032.381E-031.913122.3兩種方法的比較本為將通過計算幾個例子來比較兩種方法的優(yōu)劣性比

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

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

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

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

14、在計算一重還是二重定積分上，平均值法還是比隨機投點法性能要更好一些參考文獻：1其詩松，程依明.概率論與數(shù)理統(tǒng)計M.高等教育出版社.20112朱陸陸.蒙特卡洛方法及應用D.華中師范大學,2014.附錄：MATLABfunctiongailv(time,num,dlim,ulim)%time次數(shù)%num生成的點數(shù)%min最大值%max最小值%隨機投點法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(,隨機投點法計算定積分結果,);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(,隨機投點法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最小值%隨機投點法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);