2014數(shù)學(xué)建模之計(jì)算機(jī)仿真

1、數(shù)學(xué)建模之,計(jì)算機(jī)仿真,天目學(xué)院基礎(chǔ)部 帥昌浩 tel:18067679183,China Undergraduate Mathematical Contest in Modeling,概述,計(jì)算機(jī)科學(xué)技術(shù)的迅猛發(fā)展，給許多學(xué)科帶來(lái)了巨大的影響．計(jì)算機(jī)不但使問(wèn)題的求解變得更加方便、快捷和精確，而且使得解決實(shí)際問(wèn)題的領(lǐng)域更加廣泛．計(jì)算機(jī)適合于解決那些規(guī)模大、難以解析化以及不確定的數(shù)學(xué)模型．例如對(duì)于一些帶隨機(jī)因素的復(fù)

2、雜系統(tǒng)，用分析方法建模常常需要作許多簡(jiǎn)化假設(shè)，與面臨的實(shí)際問(wèn)題可能相差甚遠(yuǎn)，以致解答根本無(wú)法應(yīng)用，這時(shí)仿真幾乎成為人們的唯一選擇．在歷屆的美國(guó)和中國(guó)大學(xué)生的數(shù)學(xué)建模競(jìng)賽（MCM）中，學(xué)生們經(jīng)常用到計(jì)算機(jī)仿真方法去求解、檢驗(yàn)等．計(jì)算機(jī)仿真(computer simulation)是建模過(guò)程中較為重要的一類方法．,計(jì)算機(jī)仿真的基本概念,計(jì)算機(jī)仿真，是根據(jù)已知的信息和知識(shí)（如數(shù)學(xué)、物理規(guī)律），利用計(jì)算機(jī)模擬現(xiàn)實(shí)情況和系統(tǒng)的演變過(guò)程，發(fā)現(xiàn)新的知

3、識(shí)或規(guī)律，從而解決問(wèn)題的一種方法,計(jì)算機(jī)仿真的基本概念,獨(dú)立于理論研究與實(shí)驗(yàn)研究的認(rèn)識(shí)世界的第三中方法,計(jì)算機(jī)仿真的基本概念,計(jì)算機(jī)仿真的特點(diǎn),代價(jià)小、時(shí)間短、可重復(fù)、參數(shù)設(shè)置靈活,計(jì)算機(jī)仿真可以解決以下5類問(wèn)題:,(1) 難以用數(shù)學(xué)公式表示的系統(tǒng),或者沒(méi)有建立和求解的有效方法.(2)雖然可以用解析的方法解決問(wèn)題,但數(shù)學(xué)的分析與計(jì)算過(guò)于復(fù)雜,這時(shí)計(jì)算機(jī)仿真可能提供簡(jiǎn)單可行的求解方法.(3)希望能在較短的時(shí)間內(nèi)觀察到系統(tǒng)發(fā)展的全過(guò)程

4、,以估計(jì)某些參數(shù)對(duì)系統(tǒng)行為的影響.(4)難以在實(shí)際環(huán)境中進(jìn)行試驗(yàn)和觀察時(shí),計(jì)算機(jī)仿真是唯一可行的方法,如太空飛行的研究.(5)需要對(duì)系統(tǒng)或過(guò)程進(jìn)行長(zhǎng)期運(yùn)行比較,從大量方案中尋找最優(yōu)方案.,計(jì)算機(jī)仿真的分類,計(jì)算機(jī)仿真在計(jì)算機(jī)中運(yùn)行實(shí)現(xiàn),不怕破壞,易修改,可重用,安全經(jīng)濟(jì),不受外界條件和場(chǎng)地空間的限制.仿真分為靜態(tài)仿真(static simulation)和動(dòng)態(tài)仿真(dynamic simulation).數(shù)值積分中的蒙特卡洛方法（

5、統(tǒng)計(jì)模擬方法）是典型的靜態(tài)仿真．動(dòng)態(tài)仿真又分為連續(xù)系統(tǒng)仿真和離散系統(tǒng)仿真．連續(xù)系統(tǒng)是指狀態(tài)變量隨著時(shí)間連續(xù)變化的系統(tǒng)，例如傳染病的檢測(cè)與預(yù)報(bào)系統(tǒng).離散系統(tǒng)是指系統(tǒng)狀態(tài)變量只在有限的時(shí)間點(diǎn)或可數(shù)的時(shí)間點(diǎn)上發(fā)生變化的系統(tǒng),例如排隊(duì)系統(tǒng).,仿真系統(tǒng)，必須設(shè)置一個(gè)仿真時(shí)鐘(simulate clock)，它能將時(shí)間從一個(gè)時(shí)刻向另一個(gè)時(shí)刻進(jìn)行推進(jìn)，并且能隨時(shí)反映系統(tǒng)時(shí)間的當(dāng)前值．其中，模擬時(shí)間推進(jìn)方式有兩種:時(shí)間步長(zhǎng)法(均勻間隔時(shí)間推進(jìn)法,連續(xù)系

6、統(tǒng)常用)和事件步長(zhǎng)法(下次事件推進(jìn)法,離散系統(tǒng)常用) ．,主要內(nèi)容,一: 準(zhǔn)備知識(shí):隨機(jī)數(shù)的產(chǎn)生二:隨機(jī)變量的模擬三:連續(xù)系統(tǒng)的模擬-時(shí)間步長(zhǎng)法四: 離散系統(tǒng)的模擬-事件步長(zhǎng)法●五：蒙特卡洛方法,一: 準(zhǔn)備知識(shí):隨機(jī)數(shù)的產(chǎn)生,由于仿真研究的實(shí)際系統(tǒng)要受到多種隨機(jī)因素的作用和影響,在仿真過(guò)程中必須處理大量的隨機(jī)因素.要解決此問(wèn)題的前提是確定隨機(jī)變量的類型和選擇合適的隨機(jī)數(shù)產(chǎn)生的方法.對(duì)隨機(jī)現(xiàn)象進(jìn)行模擬,實(shí)質(zhì)是要給出隨機(jī)變量的模擬

7、,也就是說(shuō)要利用計(jì)算機(jī)隨機(jī)產(chǎn)生一系列數(shù)值,使它們服從一定的概率分布,稱這些數(shù)值為隨機(jī)數(shù).最基本,最常用的是(0,1)區(qū)間內(nèi)均勻分布的隨機(jī)數(shù).其他分布的隨機(jī)數(shù)均可利用它來(lái)產(chǎn)生.,1:產(chǎn)生模擬隨機(jī)數(shù)的計(jì)算機(jī)命令,在MATLAB中,可以直接產(chǎn)生滿足各種分布的隨機(jī)數(shù),命令如下:常見(jiàn)的分布函數(shù) MATLAB語(yǔ)句 均勻分布U[0,1] R=rand(m,n) 均勻分布U[a,b]

8、 R=unifrnd(a,b,m,n) 指數(shù)分布E(λ) R = exprnd(λ, m , n)正態(tài)分布N(mu,sigma) R=normrnd(mu,sigma,m,n) 標(biāo)準(zhǔn)正態(tài)分布N(0,1) R=randn(m,n) 二項(xiàng)分布B(n,p) R=binornd(n,p,m,n) 泊松分布 P(λ ) R=poi

9、ssrnd (λ, m , n)以上語(yǔ)句均產(chǎn)生m× n 的矩陣．,2:案例分析,例1: unifrnd(2,3)unifrnd(1,32,1,4)normrnd(1,2) normrnd(1,2,2,3)rand(2,3)randn(2,3),2:案例分析,ans =2.8132ans =1.3057 5.3056 7.2857 7.1604ans = 0.2527ans =

10、 2.7429 0.0219 2.7759 2.2756 0.0992 -0.9560ans = 0.6038 0.1988 0.7468 0.2722 0.0153 0.4451ans = -0.0945 -1.3089 -0.2440 -0.2141 0.8248 -0.1778>>,2:案例分析,例2:敵空戰(zhàn)部

11、隊(duì)對(duì)我方港口進(jìn)行空襲,其到達(dá)規(guī)律服從泊松分布,平均每分鐘到達(dá)4架飛機(jī).(1)模擬敵機(jī)在3分鐘內(nèi)到達(dá)目標(biāo)區(qū)域的數(shù)量,以及在第1,2,3分鐘內(nèi)各到達(dá)幾架飛機(jī);(2)模擬在3分鐘內(nèi)每架飛機(jī)的到達(dá)時(shí)刻.分析:(1) n1=poissrnd(4), n2=poissrnd(4), n3=poissrnd(4),n=n1+n2+n3 (2) 由排隊(duì)論知識(shí),敵機(jī)到達(dá)規(guī)律服從泊松分布等價(jià)于敵機(jī)到達(dá)港口的間隔時(shí)間服從參數(shù)為1/4的指數(shù)分布,故

12、可由指數(shù)分布模擬每架飛機(jī)的到達(dá)時(shí)刻.,2:案例分析,cleart=0;j=0; %到達(dá)的飛機(jī)數(shù) while t<3 j=j+1 t=t+exprnd(1/4)end,二:隨機(jī)變量的模擬,利用均勻分布的隨機(jī)數(shù)可以產(chǎn)生具有任意分布的隨機(jī)變量的樣本,從而可以對(duì)隨機(jī)變量的取值情況進(jìn)行模擬.1 連續(xù)型隨機(jī)變量的模擬具有給定分布的連續(xù)型隨機(jī)變量可以利用在區(qū)間(0

13、,1)上均勻分布的隨機(jī)數(shù)來(lái)模擬,最常用的方法是逆變換法.結(jié)論:若隨機(jī)變量Y有連續(xù)的分布函數(shù)F(y), 則Z與Y有相同的分布.,1 連續(xù)型隨機(jī)變量的模擬,若已知Y的概率密度為如果給定區(qū)間(0,1)上均勻分布的隨機(jī)數(shù) ,則具有給定分布Y的隨機(jī)數(shù) 可由方程 解出.例 :模擬服從參數(shù)為 的指數(shù)分布時(shí),由

14、 可得,2 離散型隨機(jī)變量的模擬,設(shè)隨機(jī)變量X的分布律為:,有相同的發(fā)生的概率.因此我們可以用隨機(jī)變量R落在小區(qū)間內(nèi)的情況來(lái)模擬離散的隨機(jī)變量X的取值情況.,2 離散型隨機(jī)變量的模擬,例 3 :隨機(jī)變量 表示每分鐘到達(dá)銀行柜臺(tái)的顧客數(shù).X的分布列見(jiàn)下表,試模擬10分鐘內(nèi)顧客到達(dá)柜臺(tái)的情況. 表1

15、10分鐘內(nèi)顧客到達(dá)柜臺(tái)的情況 Xk 0 1 2 pk 0.4 0.3 0.3分析:因?yàn)槊糠昼姷竭_(dá)柜臺(tái)的人數(shù)是隨機(jī)的,所以可用計(jì)算機(jī)隨機(jī)生成一組(0,1)的數(shù)據(jù),由X的概率分布情況,可認(rèn)為隨機(jī)數(shù)在(0,0.4)范圍內(nèi)時(shí)沒(méi)有顧客光顧,在[0.4,0.7)時(shí),有一個(gè)顧客光顧,在[0.7,1)時(shí),有兩個(gè)顧客光顧. 從而有MATLAB程序:,2 離散型隨機(jī)變量的模擬,r

16、=rand(1,10);for i=1:10; if r(i)<0.4 n(i)=0; elseif 0.4<=r(i)&r(i)<0.7 n(i)=1; else n(i)=2; end;

17、 end r n,三:連續(xù)系統(tǒng)的模擬-時(shí)間步長(zhǎng)法,對(duì)連續(xù)系統(tǒng)的計(jì)算機(jī)模擬是近似地獲取系統(tǒng)狀態(tài)在一些離散時(shí)刻點(diǎn)上的數(shù)值．在一定假設(shè)條件下，利用數(shù)學(xué)運(yùn)算模擬系統(tǒng)的運(yùn)行過(guò)程．連續(xù)系統(tǒng)模型一般是微分方程，它在數(shù)值模擬中最基本的算法是數(shù)值積分算法．例如有一系統(tǒng)可用微分方程來(lái)描述: 已知輸出量y的初始條件 ，現(xiàn)在要求出輸出量y隨時(shí)間變化的過(guò)程y(

18、t)。,1 理論介紹,最直觀的想法是：首先將時(shí)間離散化，令 ，稱為第k步的計(jì)算步距（一般是等間距的），然后按以下算法計(jì)算狀態(tài)變量在各時(shí)刻 上的近似值: 其中初始點(diǎn) 按照這種作法即可求出整個(gè)的曲線．這種最簡(jiǎn)單的數(shù)值積分算法稱為歐拉法．除此之外，還有其他一些算法． 􀀢,1 理

19、論介紹,因此，連續(xù)系統(tǒng)模擬方法是：首先確定系統(tǒng)的連續(xù)狀態(tài)變量，然后將它在時(shí)間上進(jìn)行離散化處理，并由此模擬系統(tǒng)的運(yùn)行狀態(tài)．模擬過(guò)程分為許多相等的時(shí)間間隔,時(shí)間步長(zhǎng)的長(zhǎng)度可以根據(jù)實(shí)際問(wèn)題分別取為秒,分,小時(shí),天等.仿真時(shí)鐘按時(shí)間步長(zhǎng)等距推進(jìn),每次推進(jìn)都要掃描系統(tǒng)中所有活動(dòng),按預(yù)定計(jì)劃和目標(biāo)進(jìn)行分析,計(jì)算,記錄系統(tǒng)狀態(tài)的變化,直到滿足某個(gè)終止條件為止.,,例:追逐問(wèn)題1．問(wèn)題提出假設(shè)正方形ABCD的四個(gè)頂點(diǎn)處各站一人．在某一時(shí)刻，四人同時(shí)

20、以勻速v沿順時(shí)針?lè)较蜃分鹣乱粋€(gè)人，并且在任意時(shí)刻他們始終保持追逐的方向是對(duì)準(zhǔn)追逐目標(biāo)，例如，A追逐B(yǎng)，任意時(shí)刻A始終向著B(niǎo)追．可以證明四人的運(yùn)動(dòng)軌跡將按螺旋曲線狀匯合于中心O． 怎樣證明呢？有兩種證明方法．一是分別求出四人的運(yùn)動(dòng)軌跡曲線解析式，求證四條曲線在某時(shí)刻相交于一點(diǎn)．另一方法則是用計(jì)算機(jī)模擬將四人的運(yùn)動(dòng)軌跡直觀地表示在圖形上.,2 應(yīng)用舉例,2．建立模型及模擬方法,模擬步驟： 1）建立平面直角坐標(biāo)系． 2）以時(shí)間間隔Δt進(jìn)行

21、采樣，在每一時(shí)t計(jì)算每個(gè)人在下一時(shí)t+Δt時(shí)的坐標(biāo)． 3）不妨設(shè)甲的追逐對(duì)象是乙，在時(shí)間t時(shí)，甲,2．建立模型及模擬方法,4)選取足夠小的 ，模擬到 時(shí)為止5）連接四人在各時(shí)刻的位置，就得到所求的軌跡．根據(jù)以上模擬步驟，編出MATLAB程序(simu2.m)如下：,%取v=1,t=12,A,B,C,D點(diǎn)的坐標(biāo)分另為（0，10），（10，10），（10，0），(0, 0)v=1;dt=0.05;

22、d=20;x=[0 0 0 10 10 10 10 0];x(9)=x(1);x(10)=x(2);holdaxis('equal')axis([0 10 0 10]);,3 MATLAB實(shí)現(xiàn),for k=1:2:7plot(x(k),x(k+1),'.' )endwhile(d>0.1)for i=1:2:7d=sqrt((x(i)-x(i+1))^2+(x(i+1)-x(i

23、+3))^2);x(i)=x(i)+v*dt*(x(i+2)-x(i))/d;x(i+1)=x(i+1)+v*dt*(x(i+3)-x(i+1))/d;plot(x(i),x(i+1),'.')endx(9)= x(1);x(10)= x(2);endhold,3 MATLAB實(shí)現(xiàn),3 MATLAB實(shí)現(xiàn),例9.6 水池含鹽量問(wèn)題,某水池有2000m3水,其中含鹽2kg,以6m3/min的速率向水池內(nèi)注入含鹽

24、為0.5kg/m3的鹽水,同時(shí)又以4m3/min的速率從水池流出攪拌均勻的鹽水.試用計(jì)算機(jī)仿真該水池內(nèi)鹽水的變化過(guò)程,并每隔10min計(jì)算水池中水的體積,含鹽量,含鹽率.欲使池中鹽水含鹽率達(dá)到0.2kg/m3,需經(jīng)過(guò)多長(zhǎng)時(shí)間?分析:這是一個(gè)連續(xù)系統(tǒng),首先要將系統(tǒng)離散化,在一些離散點(diǎn)上進(jìn)行考察,這些離散點(diǎn)的間隔就是時(shí)間步長(zhǎng).可取步長(zhǎng)為1min,即隔1min考察一次系統(tǒng)的狀態(tài),并相應(yīng)地記錄和分析.在注入和流出的作用下,池中水的體積與含

25、鹽量,含鹽率均隨時(shí)間變化,初始時(shí)刻含鹽率為0.001kg/m3,以后每分鐘注入含鹽率為0.5kg/m3的水6m3,流出混合均勻的鹽水為4m3,當(dāng)池中水的含鹽率達(dá)到0.2kg/m3時(shí),仿真過(guò)程結(jié)束.,例9.6 水池含鹽量問(wèn)題,記T時(shí)刻的體積為w(m3)，水的含鹽量為s(kg)，水的含鹽率為r=s/w(kg/m3)，每隔1min池水的動(dòng)態(tài)變化過(guò)程如下：每分鐘水的體積增加6-4=2(m3)；每分鐘向池內(nèi)注入鹽6×0.5=3(kg)；

26、每分鐘向池外流出鹽4×r(kg)；每分鐘池內(nèi)增加鹽3-4×r(kg).本例還可以用微分方程建立數(shù)學(xué)模型，并求出它的解析解，這個(gè)解析解就是問(wèn)題的精確解，有興趣的讀者可以按照這個(gè)思路求出該問(wèn)題的精確解，考察相應(yīng)時(shí)刻精確解與仿真解的差異，還可以進(jìn)一步調(diào)整仿真過(guò)程的時(shí)間步長(zhǎng)，通過(guò)與精確解的比較來(lái)研究時(shí)間步長(zhǎng)的大小對(duì)仿真度的影響。,MATLAB實(shí)現(xiàn),clearh=1; %時(shí)間步長(zhǎng)為1s0=2; %初始

27、含鹽2kgw0=2000; %初始水池有水2000m3r0=s0/w0; %初始濃度s(1)=s0+0.5*6*h-4*h*r0; %一分鐘后的含鹽量w(1)=w0+2*h; %一分鐘后水池中的鹽水體積r(1)=s(1)/w(1); %一分鐘后的濃度 t(1)=h;y(1)=(2000000+3000000*h+3000*h^2+h^3)/(1000+h)

28、^2;for i=2:200 t(i)=i*h; s(i)=s(i-1)+0.5*6*h-4*h*r(i-1); %第i步后的含鹽量 w(i)=w(i-1)+2*h; %第i步后的鹽水體積 r(i)=s(i)/w(i); %第i步后的鹽水濃度 y(i)=(2000000+3000000*t(i)+3000*t(i)^2+

29、t(i)^3)/(1000+t(i))^2; m=floor(i/10);,MATLAB實(shí)現(xiàn),if i/10-m0.2 %若第i步后的鹽水濃度大于0.2 t02=i*h; r02=r(i); break endend[t02,r02][10*tm',sm',rm'] % '表示逆subplot(1,2,

30、1),plot(t,s,'blue');hold onsubplot(1,2,2),plot(t,y,'red');,四: 離散系統(tǒng)的模擬-事件步長(zhǎng)法,離散系統(tǒng)(discrete system)是指系統(tǒng)狀態(tài)只在有限的時(shí)間點(diǎn)或可數(shù)的時(shí)間點(diǎn)上有隨機(jī)事件驅(qū)動(dòng)的系統(tǒng)．例如排隊(duì)系統(tǒng)（queue system），顯然狀態(tài)量的變化只是在離散的隨機(jī)時(shí)間點(diǎn)上發(fā)生．假設(shè)離散系統(tǒng)狀態(tài)的變化是在一個(gè)時(shí)間點(diǎn)上瞬間完成的． 常

31、用的是事件步長(zhǎng)法（下次事件推進(jìn)法）．其過(guò)程是：置模擬時(shí)鐘的初值為0，跳到第一個(gè)事件發(fā)生的時(shí)刻，計(jì)算系統(tǒng)的狀態(tài)，產(chǎn)生未來(lái)事件并加入到隊(duì)列中去，跳到下一事件，計(jì)算系統(tǒng)狀態(tài)，……，重復(fù)這一過(guò)程直到滿足某個(gè)終止條件為止．,例9.7 收款臺(tái)前的排隊(duì)過(guò)程,假設(shè):(1) 顧客到達(dá)收款臺(tái)是隨機(jī)的,平均時(shí)間間隔為0.5min,即間隔時(shí)間服從lamda=2的指數(shù)分布.(2) 對(duì)不同的顧客收款和裝袋的時(shí)間服從正態(tài)分布N(1,1/3).試模擬20位顧客到

32、收款臺(tái)前的排隊(duì)情況,我們關(guān)心的問(wèn)題是每個(gè)顧客的平均等待時(shí)間,隊(duì)長(zhǎng)及服務(wù)員的工作效率,例9.7 收款臺(tái)前的排隊(duì)過(guò)程,分析:單服務(wù)臺(tái)結(jié)構(gòu)的排隊(duì)系統(tǒng)有兩類原發(fā)事件即到來(lái)和離去,顧客到來(lái)的后繼事件是顧客接受服務(wù),顧客離去的后繼事件是服務(wù)臺(tái)尋找服務(wù),這四類事件各自的子程序框圖如圖1所示.假設(shè):t(i)為第i位顧客到達(dá)時(shí)刻;t2(i)為第i位顧客受到服務(wù)的時(shí)間(隨機(jī)變量);T(i)為第i位顧客離去時(shí)刻.將第i位顧客到達(dá)作為件事發(fā)生:t(i+1)

33、-t(i)=r(i)(隨機(jī)變量);平衡關(guān)系:當(dāng)t(i+1) ≥ T(i)時(shí), T(i+1)= t(i+1) + t2(i+1);否則, T(i+1)= T(i) + t2(i+1).,圖1 到來(lái)事件子程序,,,,系統(tǒng)人數(shù)+1,產(chǎn)生下一個(gè)顧客到來(lái)時(shí)刻,調(diào)用接受服務(wù)事件的程序,,圖 2 離去事件子程序,,系統(tǒng)人數(shù)-1,,,置服務(wù)臺(tái)”閑”,已服務(wù)人數(shù)+1,調(diào)用尋找服務(wù)事件子程序,,圖3 :接受服務(wù)事件子程序,,服務(wù)臺(tái)空,,,置服務(wù)臺(tái)”忙”,

34、產(chǎn)生服務(wù)結(jié)束時(shí)刻,登記到事件表,,,,排隊(duì)人數(shù)+1,否,是,,,圖4 :尋找服務(wù)事件子程序,,排隊(duì)人數(shù)≥1,,,置服務(wù)臺(tái)”忙”,產(chǎn)生服務(wù)結(jié)束時(shí)刻,排隊(duì)人數(shù)-1,,,,排隊(duì)人數(shù)+1,否,是,,,五：蒙特卡洛方法,在用傳統(tǒng)方法難以解決的問(wèn)題中，有很大一部分可以用概率模型進(jìn)行描述．由于這類模型含有不確定的隨機(jī)因素，分析起來(lái)通常比確定性的模型困難．有的模型難以作定量分析，得不到解析的結(jié)果，或者是雖有解析結(jié)果，但計(jì)算代價(jià)太大以至不能使用．在這種情

35、況下，可以考慮采用Monte Carlo方法。Monte Carlo方法是計(jì)算機(jī)模擬的基礎(chǔ)，它的名字來(lái)源于世界著名的賭城——摩納哥的蒙特卡洛，其歷史起源于1777年法國(guó)科學(xué)家蒲豐提出的一種計(jì)算圓周π的方法——隨機(jī)投針?lè)?，即著名的蒲豐投針問(wèn)題。,五：蒙特卡洛方法,Monte Carlo方法的基本思想是首先建立一個(gè)概率模型，使所求問(wèn)題的解正好是該模型的參數(shù)或其他有關(guān)的特征量．然后通過(guò)模擬一統(tǒng)計(jì)試驗(yàn)，即多次隨機(jī)抽樣試驗(yàn)（確定m和n），統(tǒng)計(jì)出某

36、事件發(fā)生的百分比．只要試驗(yàn)次數(shù)很大，該百分比便近似于事件發(fā)生的概率．這實(shí)際上就是概率的統(tǒng)計(jì)定義．利用建立的概率模型，求出要估計(jì)的參數(shù)．蒙特卡洛方法屬于試驗(yàn)數(shù)學(xué)的一個(gè)分支．,五：蒙特卡洛方法,蒙特卡洛方法適用范圍很廣泛，它既能求解確定性的問(wèn)題，也能求解隨機(jī)性的問(wèn)題以及科學(xué)研究中的理論問(wèn)題．例如利用蒙特卡洛方法可以近似地計(jì)算定積分，即產(chǎn)生數(shù)值積分問(wèn)題．,蒙特卡洛法求圓周率,clearn=50000X=rand(n,1);Y=rand(

37、n,1); k=0;for i=1:n; if X(i)^2+Y(i)^2<=1 k=k+1; endend4*k/n,普豐投針求圓周率的M程序,function pi_value=pinpi(k,a,l) %k投針次數(shù)%a線距%l針長(zhǎng)（必須小于等于d）%pi_value返回pi值m=0; if l>a; fprintf('error：針長(zhǎng)必須小于等于%d\

38、n',a); fprintf('請(qǐng)重新調(diào)用函數(shù)pinpi(k,d,l) \n'); pi_value=0;else for i=1:k if a*rand(1)<=l*sin(pi*rand(1)) m=m+1; end end p=m/k;pi_value=2*l/(a*p); fprintf('投針?lè)ㄇ蟮胮i=%d\n',pi_valu

39、e);end,,format long g>> pinpi(100000,4,3)投針?lè)ㄇ蟮胮i=3.143797e+000ans = 3.14379728795087,任意曲邊梯形面積的近似計(jì)算,一個(gè)古老的問(wèn)題：用一堆石頭測(cè)量一個(gè)水塘的面積．應(yīng)該怎樣做呢？測(cè)量方法如下：假定水塘位于一塊面積已知的矩形農(nóng)田之中．如圖8．2所示．隨機(jī)地向這塊農(nóng)田扔石頭使得它們都落在農(nóng)田內(nèi)．被扔到農(nóng)田中的石頭可能濺

40、上了水，也可能沒(méi)有濺上水，估計(jì)被“濺上水的”石頭量占總的石頭量的百分比．試想如何利用這估計(jì)的百分比去近似計(jì)算該水塘面積？,任意曲邊梯形面積的近似計(jì)算,任意曲邊梯形面積的近似計(jì)算,結(jié)合圖8．2中的圖形(1)分析，只要已知各種參數(shù)及函數(shù)（a，b，H，f(x)），有以下兩種方法可近似計(jì)算水塘面積． 1．隨機(jī)投點(diǎn)法 1）賦初值：試驗(yàn)次數(shù)n=0，成功次數(shù)m=0；規(guī)定投點(diǎn)試驗(yàn)的總次數(shù)N； 2）隨機(jī)選擇m個(gè)數(shù)對(duì)xi,yi,1<i<m

41、，其中ya<xi<b, 0<yi<H，置 n＝n＋l； 3）判斷n≤N，若是，轉(zhuǎn)4，否則停止計(jì)算；4）判斷條件 (表示一塊濺水的石頭)是否成立，若成立則置m=m+1，轉(zhuǎn)2，否則轉(zhuǎn)2；5）計(jì)算水塘面積的近似值S=H×(b-a)×m/N．,任意曲邊梯形面積的近似計(jì)算,2．平均值估計(jì)法1）產(chǎn)生[a,b]區(qū)間的均勻隨機(jī)數(shù)x,i=1,2…,N􀀢

42、;2) 計(jì)算f(xi）， i=1,2…,N3）計(jì)算該方法的特點(diǎn)是估計(jì)函數(shù)f(x)在[a,b]上的平均值，面積近似等于該平均值乘以(b-a)．,例9 庫(kù)存問(wèn)題,在物資的供應(yīng)過(guò)程中,由于到貨和銷售不可能做到同步同量,故總要保持一定的庫(kù)存儲(chǔ)備.如果庫(kù)存過(guò)多,就會(huì)造成積壓浪費(fèi)以及保管費(fèi)的上升;如果庫(kù)存過(guò)少,會(huì)造成缺貨.如何選擇庫(kù)存和訂貨策略,就是一個(gè)需要研究的問(wèn)題.庫(kù)存問(wèn)題有多種類型,一般都比較復(fù)雜,下面討論一種簡(jiǎn)單的情形. 某電動(dòng)車行

43、的倉(cāng)庫(kù)管理人員采取一種簡(jiǎn)單的訂貨策略,當(dāng)庫(kù)存降低到P輛電動(dòng)車時(shí)就向廠家,例9 庫(kù)存問(wèn)題,訂貨,每次訂貨Q輛,如果某一天的需求量超過(guò)了庫(kù)存量,商店就有銷售損失和信譽(yù)損失,但如果庫(kù)存量過(guò)多,就會(huì)導(dǎo)致資金積壓和保管費(fèi)增加.若現(xiàn)在有如下表所示的兩種庫(kù)存策略,試比較選擇一種策略以使總費(fèi)用最少. 表 訂貨方案 方案 重新訂貨點(diǎn)P/輛 重新訂貨量Q/輛 方案1 125

44、 150 方案2 150 250,例9 庫(kù)存問(wèn)題,這個(gè)問(wèn)題的已知條件是:(1) 從發(fā)出訂貨到收到貨物需隔3天.(2) 每輛電動(dòng)車保管費(fèi)為0.50元/天,每輛電動(dòng)車的缺貨損失為1.60元/天,每次的訂貨費(fèi)為75元.(3)每天電動(dòng)車需求量是0到99之間均勻分布的隨機(jī)數(shù).(4)原始庫(kù)存為110輛,假設(shè)第一天沒(méi)有發(fā)出訂貨.,例9 庫(kù)存問(wèn)題,分析: 這一問(wèn)題用解析法討論比較麻煩,但用計(jì)算

45、機(jī)按天仿真?zhèn)}庫(kù)貨物的變動(dòng)情況卻很方便.我們以30天為例,依次對(duì)這兩種方案進(jìn)行仿真,最后比較各方案的總費(fèi)用,從而就可以作出決策.計(jì)算機(jī)仿真時(shí)的工作流程是早上到貨,全天銷售,晚上訂貨,輸入一些常數(shù)和初始數(shù)據(jù)后,以一天為時(shí)間步長(zhǎng)進(jìn)行仿真.首先檢查這一天是否為預(yù)訂到貨日期,如果是,則原有庫(kù)存量加Q,并把預(yù)定到貨量清為0;如果不是,則庫(kù)存量不變.接著用計(jì)算機(jī)中的隨機(jī)函數(shù)仿真隨機(jī)需求量,若庫(kù)存量大于需求量,則新的庫(kù)存量減去需求量;反之,則新庫(kù)存量

46、變?yōu)?,并且要在總費(fèi)用上加缺貨損失,然后檢查實(shí)際,例9 庫(kù)存問(wèn)題,庫(kù)存量加上預(yù)定到貨量是否小于重新訂貨點(diǎn)P,如果是,則需要重新訂貨,這時(shí)就加一次訂貨費(fèi).如此重復(fù)運(yùn)行30天,即可得所需費(fèi)用總值.由此比較這兩種方案的總費(fèi)用,即得最好方案.,,cleardays=30;P=[125,150];Q=[150,250];cost=[0,0];arrivalinterval=2;storagefee=0.5;lossfee=1.6;bo

47、okfee=75;storage0=110;booknumber=0;arrivedate=0;nr=rand(days,1);for i=1:2 storage(1)=storage0; n=round(99*nr(1)); sale=n; remain=storage(1)-n; if remain<=P(i); booknumber=Q(i

48、); arrivedate=4; orderfee=bookfee; else orderfee=0; end storage(1)=remain; cost(i)=cost(i)+remain*storagefee+orderfee; for j=2:days dh=j;

49、 if abs(dh-arrivedate)<0.01 storage(j)=storage(j-1)+booknumber; booknumber=0;,,arrivedate=j; else storage(j)=storage(j-1);

50、 end n=round(99*nr(j)); if storage(j)>=n sale=n; remain=storage(j)-n; shortagenumber=0; else

51、 sale=storage(j); remain=0; shortagenumber=n-storage(j); end storage(j)=remain; if remain+booknumber<=P(i);

52、 booknumber=Q(i); arrivedate=dh+arrivalinterval; orderfee=bookfee; else orderfee=0; end cost(i)=cost(i)+remain*sto

53、ragefee+shortagenumber*lossfee+orderfee; end; mincost=min(cost); end cost mincost,例9.10 趕火車過(guò)程仿真,一列火車從A站經(jīng)B站開(kāi)往C站,某人每天趕往B站乘這趟火車.已知火車從A站到B站運(yùn)行時(shí)間為均值30min,標(biāo)準(zhǔn)差為2min的正態(tài)隨機(jī)變量.火車大約在下午1點(diǎn)離開(kāi)A站,離開(kāi)時(shí)刻的頻率分布見(jiàn)表1.這

54、個(gè)人到達(dá)B站時(shí)的頻率分布表見(jiàn)表2.用計(jì)算機(jī)仿真火車開(kāi)出,火車到達(dá)B站,這個(gè)人到達(dá)B站的情況,并給出能趕上火車的仿真結(jié)果. 表1 離開(kāi)時(shí)刻的頻率分布出發(fā)時(shí)刻T 1:00 1:05 1:10頻率 0.7 0.2 0.1,例9.10 趕火車過(guò)程仿真,表2 人到達(dá)B站時(shí)的頻率分布 到達(dá)時(shí)刻T 1:28 1:30 1:32 1:34 頻

55、率: 0.3 0.4 0.2 0.1分析:引入以下變量:T1為火車從A站開(kāi)出的時(shí)刻;T2為火車從A站運(yùn)行到B站所需要的時(shí)間;T3為此人到達(dá)B站的時(shí)刻.則有表3. 表3 T1, T2的頻率分布T1,/min 0 5 10 T2/min 28 30 32 34 p 0.7 0.2 0.1 0.3 0.4 0.2 0.1,例9.10 趕火車過(guò)程

56、仿真,顯然,這位旅客要趕上火車的條件是T3< T1+ T2,可以通過(guò)計(jì)算機(jī)模擬出這三個(gè)時(shí)間,再檢驗(yàn)是否滿足T3< T1+ T2.如果滿足,即能夠趕上火車,若不然,則不能趕上火車.,火車運(yùn)行時(shí)間的仿真程序,x=randn(10000,1);for i=1:10000 y(i)=30+2*x(i);end,開(kāi)車時(shí)間的仿真程序,s1=0; s2=0; s3=0; x=rand(10000,1);

57、 for i=1:10000 if x(i)0.9 s3=s3+1; end end[s1/10000, 1-s1/10000-s3/10000,s3/10000],人到達(dá)時(shí)刻的仿真程序,s1=0; s2=0; s3=0;s4=0; x=rand(10000,1); for i=1:10000 if x(i)<

58、;0.3 s1=s1+1; elseif x(i)<0.7 s2=s2+1; else if x(i)<0.9 s3=s3+1; else s4=s4+1; end

59、 endend[s1/10000, s2/10000,s3/10000,s4/10000],趕上火車的仿真程序,s=0;x1=rand(10000,1);x2=rand(10000,1);x3=randn(10000,1);for i=1:10000if x1(i)<0.7T1=0;elseif x1(i)<0.9T1=5; else T1=10;

60、 endT2=30+2*x3(i);if x2(i)<0.3 T3=28;elseif x2(i)<0.7 T3=30; elseif x2(i)<0.9 T3=32; else T3=34; endif T3<T1+T2