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

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

一個問題我們做一個實(shí)驗(yàn):把一個硬幣擲一萬次,統(tǒng)計兩個面出現(xiàn)的次數(shù)。這樣做很簡單但卻需要大師時間，有沒有一咱較快的辦法把這個實(shí)驗(yàn)完成呢？擲硬幣仿真流程圖概述

計算機(jī)科學(xué)技術(shù)的迅猛發(fā)展，給許多學(xué)科帶來了巨大的影響．計算機(jī)不但使問題的求解變得更加方便、快捷和精確，而且使得解決實(shí)際問題的領(lǐng)域更加廣泛．計算機(jī)適合于解決那些規(guī)模大、難以解析化以及不確定的數(shù)學(xué)模型．例如對于一些帶隨機(jī)因素的復(fù)雜系統(tǒng)，用分析方法建模常常需要作許多簡化假設(shè)，與面臨的實(shí)際問題可能相差甚遠(yuǎn)，以致解答根本無法應(yīng)用，這時仿真幾乎成為人們的唯一選擇．在歷屆的美國和中國大學(xué)生的數(shù)學(xué)建模競賽（MCM）中，學(xué)生們經(jīng)常用到計算機(jī)仿真方法去求解、檢驗(yàn)等．計算機(jī)仿真(computersimulation)是建模過程中較為重要的一類方法．計算機(jī)仿真可以解決以下5類問題:(1)難以用數(shù)學(xué)公式表示的系統(tǒng),或者沒有建立和求解的有效方法.(2)雖然可以用解析的方法解決問題,但數(shù)學(xué)的分析與計算過于復(fù)雜,這時計算機(jī)仿真可能提供簡單可行的求解方法.(3)希望能在較短的時間內(nèi)觀察到系統(tǒng)發(fā)展的全過程,以估計某些參數(shù)對系統(tǒng)行為的影響.(4)難以在實(shí)際環(huán)境中進(jìn)行試驗(yàn)和觀察時,計算機(jī)仿真是唯一可行的方法,如太空飛行的研究.(5)需要對系統(tǒng)或過程進(jìn)行長期運(yùn)行比較,從大量方案中尋找最優(yōu)方案.計算機(jī)仿真案例1模型建立：由于本題要求使從攪拌中心到各個工地運(yùn)輸混凝土的總的噸公里數(shù)最少，所以，該問題的目標(biāo)函數(shù)是求解方法：1、高數(shù)中的方法2、數(shù)值計算方法3、計算機(jī)仿真：離散化，遍歷！計算機(jī)仿真案例2例2（趕火車過程仿真）一列火車從A站經(jīng)過B站開往C站，某人每天趕往B站乘這趟火車。已知火車從A站到B站的運(yùn)行時間是均值為30min、標(biāo)準(zhǔn)差為2min的正態(tài)隨機(jī)變量?；疖嚧蠹s在下午1點(diǎn)離開Ａ站?；疖囯x開時刻的頻率分布和這個人到達(dá)Ｂ站時刻的頻率分布如下表所示。問他能趕上火車的概率有多大？

出發(fā)時刻1:001:051:10到達(dá)時刻1:281:301:321:34頻率0.70.20.1頻率0.30.40.20.1仿真過程：1、生成火車的發(fā)車時間、運(yùn)行時間，從而達(dá)得到其到達(dá)B站的時間。2、生成此人達(dá)到B站的時間。3、如果此人到達(dá)Ｂ站的時間早于火車到達(dá)時間，則算趕上火車一次。4、將上述過程重復(fù)一萬次，統(tǒng)計趕上火車的頻率作為所求概率。分析：這個問題用概率論的方法求解十分困難，它涉及此人到達(dá)時刻、火車離開A站的時刻、火車運(yùn)行時間幾個隨機(jī)變量。我們可以用計算機(jī)仿真的方法來解決。概述計算機(jī)仿真在計算機(jī)中運(yùn)行實(shí)現(xiàn),不怕破壞,易修改,可重用,安全經(jīng)濟(jì),不受外界條件和場地空間的限制.仿真分為靜態(tài)仿真(staticsimulation)和動態(tài)仿真(dynamicsimulation).數(shù)值積分中的蒙特卡洛方法是典型的靜態(tài)仿真．動態(tài)仿真又分為連續(xù)系統(tǒng)仿真和離散系統(tǒng)仿真．連續(xù)系統(tǒng)是指狀態(tài)變量隨著時間連續(xù)變化的系統(tǒng)，例如傳染病的檢測與預(yù)報系統(tǒng).離散系統(tǒng)是指系統(tǒng)狀態(tài)變量只在有限的時間點(diǎn)或可數(shù)的時間點(diǎn)上發(fā)生變化的系統(tǒng),例如排隊(duì)系統(tǒng).概述仿真系統(tǒng)，必須設(shè)置一個仿真時鐘(simulateclock)，它能將時間從一個時刻向另一個時刻進(jìn)行推進(jìn)，并且能隨時反映系統(tǒng)時間的當(dāng)前值．其中，模擬時間推進(jìn)方式有兩種:時間步長法(均勻間隔時間推進(jìn)法,連續(xù)系統(tǒng)常用)和事件步長法(下次事件推進(jìn)法,離散系統(tǒng)常用)．主要內(nèi)容一:準(zhǔn)備知識:隨機(jī)數(shù)的產(chǎn)生二:隨機(jī)變量的模擬三:連續(xù)系統(tǒng)的模擬-時間步長法四:離散系統(tǒng)的模擬-事件步長法五：蒙特卡洛方法一:準(zhǔn)備知識:隨機(jī)數(shù)的產(chǎn)生由于仿真研究的實(shí)際系統(tǒng)要受到多種隨機(jī)因素的作用和影響,在仿真過程中必須處理大量的隨機(jī)因素.要解決此問題的前提是確定隨機(jī)變量的類型和選擇合適的隨機(jī)數(shù)產(chǎn)生的方法.對隨機(jī)現(xiàn)象進(jìn)行模擬,實(shí)質(zhì)是要給出隨機(jī)變量的模擬,也就是說要利用計算機(jī)隨機(jī)產(chǎn)生一系列數(shù)值,使它們服從一定的概率分布,稱這些數(shù)值為隨機(jī)數(shù).最基本,最常用的是(0,1)區(qū)間內(nèi)均勻分布的隨機(jī)數(shù).其他分布的隨機(jī)數(shù)均可利用它來產(chǎn)生.1:產(chǎn)生模擬隨機(jī)數(shù)的計算機(jī)命令在MATLAB中,可以直接產(chǎn)生滿足各種分布的隨機(jī)數(shù),命令如下:常見的分布函數(shù)MATLAB語句均勻分布U[0,1]R=rand(m,n)均勻分布U[a,b]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,n1)泊松分布P(λ)R=poissrnd(λ,m,n)以上語句均產(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)ans=2.8132ans=1.30575.30567.28577.1604ans=0.2527ans=2.74290.02192.77592.27560.0992-0.9560ans=0.60380.19880.74680.27220.01530.4451ans=-0.0945-1.3089-0.2440-0.21410.8248-0.1778>>2:案例分析例2:敵空戰(zhàn)部隊(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á)時刻.分析:(1)n1=poissrnd(4),n2=poissrnd(4),n3=poissrnd(4),n=n1+n2+n3(2)由排隊(duì)論知識,敵機(jī)到達(dá)規(guī)律服從泊松分布等價于敵機(jī)到達(dá)港口的間隔時間服從參數(shù)為1/4的指數(shù)分布,故可由指數(shù)分布模擬每架飛機(jī)的到達(dá)時刻.注：如果單位時間發(fā)生的次數(shù)（如到達(dá)的人數(shù)）服從參數(shù)為r的泊松分布,則任連續(xù)發(fā)生的兩次時間的間隔時間序列服從參數(shù)為r的指數(shù)分布!泊松分布的期望是λ，根據(jù)到泊松分布和指數(shù)分布的關(guān)系，可以推出指數(shù)分布的期望是1/λ。

2:案例分析cleart=0;j=0;%到達(dá)的飛機(jī)數(shù)whilet<3j=j+1t=t+exprnd(1/4)end二:隨機(jī)變量的模擬在很多實(shí)際問題中，我們需要模擬服從一定分布的隨機(jī)變量,來進(jìn)行計算和預(yù)測。利用均勻分布的隨機(jī)數(shù)可以產(chǎn)生具有任意分布的隨機(jī)變量的樣本,從而可以對隨機(jī)變量的取值情況進(jìn)行模擬.1連續(xù)型隨機(jī)變量的模擬具有給定分布的連續(xù)型隨機(jī)變量可以利用在區(qū)間(0,1)上均勻分布的隨機(jī)數(shù)來模擬,最常用的方法是逆變換法.結(jié)論:若隨機(jī)變量Y有連續(xù)的分布函數(shù)F(y),

則Z與Y有相同的分布.1連續(xù)型隨機(jī)變量的模擬一般說來，具有給定分布的連續(xù)型隨機(jī)變量可以利用在區(qū)間(0,1)上均勻分布的隨機(jī)數(shù)來模擬．最常用的方法是反函數(shù)法．由概率論的理論可以證明，若隨機(jī)變量Y有連續(xù)的分布函數(shù)F(y)，而X是區(qū)間(0,1)上均勻分布的隨機(jī)變量，令

，則Z與Y有相同的分布．由此，若已知Y的概率密度為

，由

可得反函數(shù)法如果給定區(qū)間(0,1)上均勻分布的隨機(jī)數(shù),則具有給定分布Y的隨機(jī)數(shù)可由方程解出.例:模擬服從參數(shù)為的指數(shù)分布時,由可得注：指數(shù)分布的密度函數(shù)2離散型隨機(jī)變量的模擬設(shè)隨機(jī)變量X的分布律為:有相同的發(fā)生的概率.因此我們可以用隨機(jī)變量R落在小區(qū)間內(nèi)的情況來模擬離散的隨機(jī)變量X的取值情況.因此我們可以用隨機(jī)變量R落在小區(qū)間內(nèi)的情況來模擬離散的隨機(jī)變量X的取值情況．具體執(zhí)行的過程是：產(chǎn)生一個(0,1)上均勻分布的隨機(jī)數(shù)r(簡稱隨機(jī)數(shù))，若p(n-1)<r≤p(n)則理解為發(fā)生事件(X=xn)．于是就可以模擬隨機(jī)變量的取值情況．2離散型隨機(jī)變量的模擬例3:隨機(jī)變量表示每分鐘到達(dá)銀行柜臺的顧客數(shù).X的分布列見下表,試模擬10分鐘內(nèi)顧客到達(dá)柜臺的情況.

表110分鐘內(nèi)顧客到達(dá)柜臺的情況

Xk012

pk0.40.30.3分析:因?yàn)槊糠昼姷竭_(dá)柜臺的人數(shù)是隨機(jī)的,所以可用計算機(jī)隨機(jī)生成一組(0,1)的數(shù)據(jù),由X的概率分布情況,可認(rèn)為隨機(jī)數(shù)在(0,0.4)范圍內(nèi)時沒有顧客光顧,在[0.4,0.7)時,有一個顧客光顧,在[0.7,1)時,有兩個顧客光顧.

從而有MATLAB程序:2離散型隨機(jī)變量的模擬r=rand(1,10);fori=1:10;ifr(i)<0.4

n(i)=0;

elseif0.4<=r(i)&r(i)<0.7

n(i)=1;elsen(i)=2;end;endrn三:連續(xù)系統(tǒng)的模擬-時間步長法對連續(xù)系統(tǒng)的計算機(jī)模擬是近似地獲取系統(tǒng)狀態(tài)在一些離散時刻點(diǎn)上的數(shù)值．在一定假設(shè)條件下，利用數(shù)學(xué)運(yùn)算模擬系統(tǒng)的運(yùn)行過程．連續(xù)系統(tǒng)模型一般是微分方程，它在數(shù)值模擬中最基本的算法是數(shù)值積分算法．例如有一系統(tǒng)可用微分方程來描述:已知輸出量y的初始條件，現(xiàn)在要求出輸出量y隨時間變化的過程y(t)。1理論介紹最直觀的想法是：首先將時間離散化，令，稱為第k步的計算步距（一般是等間距的），然后按以下算法計算狀態(tài)變量在各時刻上的近似值:

其中初始點(diǎn)按照這種作法即可求出整個的曲線．這種最簡單的數(shù)值積分算法稱為歐拉法．除此之外，還有其他一些算法．1理論介紹因此，連續(xù)系統(tǒng)模擬方法是：首先確定系統(tǒng)的連續(xù)狀態(tài)變量，然后將它在時間上進(jìn)行離散化處理，并由此模擬系統(tǒng)的運(yùn)行狀態(tài)．模擬過程分為許多相等的時間間隔,時間步長的長度可以根據(jù)實(shí)際問題分別取為秒,分,小時,天等.仿真時鐘按時間步長等距推進(jìn),每次推進(jìn)都要掃描系統(tǒng)中所有活動,按預(yù)定計劃和目標(biāo)進(jìn)行分析,計算,記錄系統(tǒng)狀態(tài)的變化,直到滿足某個終止條件為止.計算機(jī)仿真舉例例4：追擊問題

我輯私艦雷達(dá)發(fā)現(xiàn)距c里處有一艘走私船正以勻速度a沿直線行駛,輯私艦立即以最大的速度b追趕,若用雷達(dá)進(jìn)行跟蹤,保持船的瞬時速度方向始終指向走私船,試求輯私艦追逐路線和追上的時間,并且給出緝私艦和走私船的路線圖。

選取坐標(biāo)系如圖:假設(shè)緝私艇初始位置在點(diǎn)(c,0),走私船初始位置在點(diǎn)(0,0),走私船的行駛方向?yàn)閥方向.Oxy(c,0)D(x,y).R(0,at)設(shè)緝私艇為動點(diǎn)D,走私船為動點(diǎn)R.在時刻t,緝私艇的位置是D(x,y),走私船的位置是R(0,at).取時間間隔(步長)為,則在時刻t+,D的位置是,追趕方向(D的運(yùn)動方向)為DR.用方向余弦表示為:(*)計算機(jī)仿真舉例算法:賦初值:初始時刻t0,時間步長,速度a,b,初始位置c循環(huán):由公式(*)計算動點(diǎn)D在各時刻的坐標(biāo)

D

計算動點(diǎn)R在各時刻的坐標(biāo)R.

終止:

當(dāng)D,R的距離小于給定值時終止仿真算法：

第二步：計算動點(diǎn)緝私艇D在時刻

時的坐標(biāo)

計算走私船R在時刻

時的坐標(biāo)第一步：設(shè)置時間步長,速度a,b及初始位置第三步：計算緝私艇與走私船這兩個動點(diǎn)之間的距離：

根據(jù)事先給定的距離，判斷緝私艇是否已經(jīng)追上了走私船，從而判斷退出循環(huán)還是讓時間產(chǎn)生一個步長，返回到第二步繼續(xù)進(jìn)入下一次循環(huán)；

第四步：當(dāng)從上述循環(huán)退出后，由點(diǎn)列和可分別繪制成兩條曲線即為緝私艇和走私船走過的軌跡曲線。

試求緝私艦追逐路線和追上的時間。取c=3千米，a=0.4千米/分鐘，b=0.8千米/分鐘計算機(jī)仿真舉例將上述算法中求出的點(diǎn)D用直線段連接起來，就得到追趕曲線。循環(huán)終止時的即為追趕時間?？梢钥闯觯介L選取得越小，所求的曲線越精確。

例5:追逐問題1．問題提出假設(shè)正方形ABCD的四個頂點(diǎn)處各站一人．在某一時刻，四人同時以勻速v沿順時針方向追逐下一個人，并且在任意時刻他們始終保持追逐的方向是對準(zhǔn)追逐目標(biāo)，例如，A追逐B(yǎng)，任意時刻A始終向著B追．可以證明四人的運(yùn)動軌跡將按螺旋曲線狀匯合于中心O．怎樣證明呢？有兩種證明方法．一是分別求出四人的運(yùn)動軌跡曲線解析式，求證四條曲線在某時刻相交于一點(diǎn)．另一方法則是用計算機(jī)模擬將四人的運(yùn)動軌跡直觀地表示在圖形上.2應(yīng)用舉例2．建立模型及模擬方法模擬步驟：1）建立平面直角坐標(biāo)系．

2）以時間間隔tΔ進(jìn)行采樣，在每一時t計算每個人在下一時t+tΔ時的坐標(biāo)．

3）不妨設(shè)甲的追逐對象是乙，在時間t時，甲2．建立模型及模擬方法4)選取足夠小的，模擬到時為止(即已追到）5）連接四人在各時刻的位置，就得到所求的軌跡．根據(jù)以上模擬步驟，編出MATLAB程序如下：%取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;d=20;x=[0001010

10

100];%DABC點(diǎn)x,y坐標(biāo)x(9)=x(1);x(10)=x(2);%D點(diǎn)holdaxis('equal')%各坐標(biāo)軸刻度增量相同axis([010010]);%坐標(biāo)軸范圍3MATLAB實(shí)現(xiàn)fork=1:2:7plot(x(k),x(k+1),'.')endwhile(d>0.1)fori=1:2:7%循環(huán)求每個點(diǎn)的下一個時間點(diǎn)的坐標(biāo)d=sqrt((x(i)-x(i+2))^2+(x(i+1)-x(i+3))^2);%與追擊點(diǎn)的距離x(i)=x(i)+v*dt*(x(i+2)-x(i))/d;%下一個時間點(diǎn)的x坐標(biāo)x(i+1)=x(i+1)+v*dt*(x(i+3)-x(i+1))/d;%下一個時間點(diǎn)的y坐標(biāo)plot(x(i),x(i+1),'.')endx(9)=x(1);x(10)=x(2);endhold3MATLAB實(shí)現(xiàn)3MATLAB實(shí)現(xiàn)例6水池含鹽量問題某水池有2000m3水,其中含鹽2kg,以6m3/min的速率向水池內(nèi)注入含鹽為0.5kg/m3的鹽水,同時又以4m3/min的速率從水池流出攪拌均勻的鹽水.試用計算機(jī)仿真該水池內(nèi)鹽水的變化過程,并每隔10min計算水池中水的體積,含鹽量,含鹽率.欲使池中鹽水含鹽率達(dá)到0.2kg/m3,需經(jīng)過多長時間?分析:這是一個連續(xù)系統(tǒng),首先要將系統(tǒng)離散化,在一些離散點(diǎn)上進(jìn)行考察,這些離散點(diǎn)的間隔就是時間步長.可取步長為1min,即隔1min考察一次系統(tǒng)的狀態(tài),并相應(yīng)地記錄和分析.在注入和流出的作用下,池中水的體積與含鹽量,含鹽率均隨時間變化,初始時刻含鹽率為0.001kg/m3,以后每分鐘注入含鹽率為0.5kg/m3的水6m3,流出混合均勻的鹽水為4m3,當(dāng)池中水的含鹽率達(dá)到0.2kg/m3時,仿真過程結(jié)束.例6水池含鹽量問題記T時刻的體積為w(m3)，水的含鹽量為s(kg)，水的含鹽率為r=s/w(kg/m3)，每隔1min池水的動態(tài)變化過程如下：每分鐘水的體積增加6-4=2(m3)；每分鐘向池內(nèi)注入鹽6×0.5=3(kg)；每分鐘向池外流出鹽4×r(kg)；每分鐘池內(nèi)增加鹽3-4×r(kg).本例還可以用微分方程建立數(shù)學(xué)模型，并求出它的解析解，這個解析解就是問題的精確解，有興趣的讀者可以按照這個思路求出該問題的精確解，考察相應(yīng)時刻精確解與仿真解的差異，還可以進(jìn)一步調(diào)整仿真過程的時間步長，通過與精確解的比較來研究時間步長的大小對仿真度的影響。MATLAB實(shí)現(xiàn)clearh=1;%時間步長為1s0=2;%初始含鹽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)^2;fori=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+t(i)^3)/(1000+t(i))^2;m=floor(i/10);MATLAB實(shí)現(xiàn)

ifi/10-m<0.1

tm(m)=m;

wm(m)=w(i);

sm(m)=s(i);

rm(m)=r(i);endifr(i)>0.2%若第i步后的鹽水濃度大于0.2t02=i*h;r02=r(i);breakendend[t02,r02][10*tm',sm',rm']%'表示逆subplot(1,2,1),plot(t,s,'blue');holdonsubplot(1,2,2),plot(t,y,'red');應(yīng)用舉例-面積計算例7:面積計算

求由曲線x2+y2=16,x2/36+y2=1,以及(x-2)2+(y+1)2=9所圍成圖形的面積。

用Matlab語句作出圖形得區(qū)域圖t=0:0.01:2*pi;x=sin(t);y=cos(t);plot(4*x,4*y,6*x,y,2+3*x,3*y-1)axis([-6,6,-6,6])應(yīng)用舉例-面積計算應(yīng)用舉例-面積計算此區(qū)域形狀復(fù)雜，理論分析困難，可以用計算機(jī)仿真實(shí)現(xiàn)。將可能的區(qū)域等分，考察每個小區(qū)域是否在此區(qū)域中，將在此區(qū)域中的小面積相加即可其仿真圖如下：應(yīng)用舉例-面積計算Matlab程序?yàn)閤=-2:0.01:6;y=-2:0.01:2;s=0;h=0.01;fori=1:800forj=1:400xx=-2+i*h;

yy=-2+j*h;ifxx^2+yy^2<=16ifxx^2/36+yy^2<=1if(xx-2)^2+(yy+1)^2<=9s=s+h^2;endendendendends運(yùn)行后給出面積的值

8.8310。

四:離散系統(tǒng)的模擬-事件步長法

離散系統(tǒng)(discretesystem)是指系統(tǒng)狀態(tài)只在有限的時間點(diǎn)或可數(shù)的時間點(diǎn)上有隨機(jī)事件驅(qū)動的系統(tǒng)．例如排隊(duì)系統(tǒng)（queuesystem），顯然狀態(tài)量的變化只是在離散的隨機(jī)時間點(diǎn)上發(fā)生．假設(shè)離散系統(tǒng)狀態(tài)的變化是在一個時間點(diǎn)上瞬間完成的．常用的是事件步長法（下次事件推進(jìn)法）．其過程是：置模擬時鐘的初值為0，跳到第一個事件發(fā)生的時刻，計算系統(tǒng)的狀態(tài)，產(chǎn)生未來事件并加入到隊(duì)列中去，跳到下一事件，計算系統(tǒng)狀態(tài)，……，重復(fù)這一過程直到滿足某個終止條件為止．例7收款臺前的排隊(duì)過程假設(shè):

(1)顧客到達(dá)收款臺是隨機(jī)的,平均時間間隔為0.5min,即間隔時間服從lamda=2的指數(shù)分布.(2)對不同的顧客收款和裝袋的時間服從正態(tài)分布N(1,1/3).試模擬20位顧客到收款臺前的排隊(duì)情況,我們關(guān)心的問題是每個顧客的平均等待時間,隊(duì)長及服務(wù)員的工作效率例7收款臺前的排隊(duì)過程分析:單服務(wù)臺結(jié)構(gòu)的排隊(duì)系統(tǒng)有兩類原發(fā)事件即到來和離去,顧客到來的后繼事件是顧客接受服務(wù),顧客離去的后繼事件是服務(wù)臺尋找服務(wù),這四類事件各自的子程序框圖如圖1所示.假設(shè):t(i)為第i位顧客到達(dá)時刻;t2(i)為第i位顧客受到服務(wù)的時間(隨機(jī)變量);T(i)為第i位顧客離去時刻.將第i位顧客到達(dá)作為事件發(fā)生:t(i+1)-t(i)=r(i)(隨機(jī)變量);平衡關(guān)系:當(dāng)t(i+1)≥T(i)時,T(i+1)=t(i+1)+t2(i+1);否則,T(i+1)=T(i)+t2(i+1).圖1到來事件子程序

系統(tǒng)人數(shù)+1產(chǎn)生下一個顧客到來時刻調(diào)用接受服務(wù)事件的程序圖2離去事件子程序

系統(tǒng)人數(shù)-1置服務(wù)臺”閑”已服務(wù)人數(shù)+1調(diào)用尋找服務(wù)事件子程序圖3:接受服務(wù)事件子程序

服務(wù)臺空置服務(wù)臺”忙”產(chǎn)生服務(wù)結(jié)束時刻登記到事件表排隊(duì)人數(shù)+1否是圖4:尋找服務(wù)事件子程序

排隊(duì)人數(shù)≥1置服務(wù)臺”忙”產(chǎn)生服務(wù)結(jié)束時刻排隊(duì)人數(shù)-1排隊(duì)人數(shù)+1否是五：蒙特卡洛方法在用傳統(tǒng)方法難以解決的問題中，有很大一部分可以用概率模型進(jìn)行描述．由于這類模型含有不確定的隨機(jī)因素，分析起來通常比確定性的模型困難．有的模型難以作定量分析，得不到解析的結(jié)果，或者是雖有解析結(jié)果，但計算代價太大以至不能使用．在這種情況下，可以考慮采用MonteCarlo方法。MonteCarlo方法是計算機(jī)模擬的基礎(chǔ)，它的名字來源于世界著名的賭城——摩納哥的蒙特卡洛，其歷史起源于1777年法國科學(xué)家蒲豐提出的一種計算圓周π的方法——隨機(jī)投針法，即著名的蒲豐投針問題。18世紀(jì)，法國數(shù)學(xué)家蒲豐和勒可萊爾提出的“投針問題”，記載于蒲豐1777年出版的著作中：“在平面上畫有一組間距為a的平行線，將一根長度為L（L=a/2）的針任意擲在這個平面上，求此針與平行線中任一條相交的概率?！逼沿S本人證明了，這個概率是：p=2L/（πa）π為圓周率利用這個公式可以用概率的方法得到圓周率的近似值。下面是一些資料試驗(yàn)者時間投擲次數(shù)相交次數(shù)圓周率估計值Wolf1850年500025323.1596Smith1855年32041218.53.1554C.DeMorgan1860年600382.53.137Fox1884年10304893.1595Lazzerini1901年340818083.1415929Reina1925年25208593.1795五：蒙特卡洛方法MonteCarlo方法的基本思想是首先建立一個概率模型，使所求問題的解正好是該模型的參數(shù)或其他有關(guān)的特征量．然后通過模擬一統(tǒng)計試驗(yàn)，即多次隨機(jī)抽樣試驗(yàn)（確定m和n），統(tǒng)計出某事件發(fā)生的百分比．只要試驗(yàn)次數(shù)很大，該百分比便近似于事件發(fā)生的概率．這實(shí)際上就是概率的統(tǒng)計定義．利用建立的概率模型，求出要估計的參數(shù)．蒙特卡洛方法屬于試驗(yàn)數(shù)學(xué)的一個分支．五：蒙特卡洛方法蒙特卡洛方法適用范圍很廣泛，它既能求解確定性的問題，也能求解隨機(jī)性的問題以及科學(xué)研究中的理論問題．例如利用蒙特卡洛方法可以近似地計算定積分，即產(chǎn)生數(shù)值積分問題．蒙特卡洛法求圓周率clearn=50000X=rand(n,1);Y=rand(n,1);k=0;fori=1:n;ifX(i)^2+Y(i)^2<=1k=k+1;endend4*k/n%1/4*PI*r^2==k/n蒲豐投針求圓周率的M程序functionpi_value=pinpi(k,a,l)%k投針次數(shù)%a線距%l針長（必須小于等于a）%pi_value返回pi值m=0;ifl>a;

fprintf('error：針長必須小于等于%d\n',a);

fprintf('請重新調(diào)用函數(shù)pinpi(k,d,l)\n');

pi_value=0;elsefori=1:ka*rand(1)<=l*sin(pi*rand(1))%判斷投針下沿與最近平行線的距離與投針垂直高度的關(guān)系，距離小于等于高度就相交

ifm=m+1;endendp=m/k;pi_value=2*l/(a*p);fprintf('投針法求得pi=%d\n',pi_value);endformatlongg>>pinpi(100000,4,3)投針法求得pi=3.143797e+000ans=3.14379728795087任意曲邊梯形面積的近似計算

一個古老的問題：用一堆石頭測量一個水塘的面積．應(yīng)該怎樣做呢？測量方法如下：假定水塘位于一塊面積已知的矩形農(nóng)田之中．如圖8．2所示．隨機(jī)地向這塊農(nóng)田扔石頭使得它們都落在農(nóng)田內(nèi)．被扔到農(nóng)田中的石頭可能濺上了水，也可能沒有濺上水，估計被“濺上水的”石頭量占總的石頭量的百分比．試想如何利用這估計的百分比去近似計算該水塘面積？

任意曲邊梯形面積的近似計算

任意曲邊梯形面積的近似計算

結(jié)合圖8．2中的圖形(1)分析，只要已知各種參數(shù)及函數(shù)（a，b，H，f(x)），有以下兩種方法可近似計算水塘面積．1．隨機(jī)投點(diǎn)法1）賦初值：試驗(yàn)次數(shù)n=0，成功次數(shù)m=0；規(guī)定投點(diǎn)試驗(yàn)的總次數(shù)N；2）隨機(jī)選擇m個數(shù)對xi,yi,1<i<m，其中a<xi<b,

0<yi<H，置n＝n＋l；3）判斷n≤N，若是，轉(zhuǎn)4，否則停止計算；4）判斷條件(表示一塊濺水的石頭)是否成立，若成立則置m=m+1，轉(zhuǎn)2，否則轉(zhuǎn)2；5）計算水塘面積的近似值S=H×(b-a)×m/N．任意曲邊梯形面積的近似計算

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

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

表訂貨方案方案重新訂貨點(diǎn)P/輛重新訂貨量Q/輛方案1125150

方案2150250例9庫存問題這個問題的已知條件是:(1)從發(fā)出訂貨到收到貨物需隔3天.(2)每輛電動車保管費(fèi)為0.50元/天,每輛電動車的缺貨損失為1.60元/天,每次的訂貨費(fèi)為75元.(3)每天電動車需求量是0到99之間均勻分布的隨機(jī)數(shù).(4)原始庫存為110輛,假設(shè)第一天沒有發(fā)出訂貨.例9庫存問題分析:這一問題用解析法討論比較麻煩,但用計算機(jī)按天仿真?zhèn)}庫貨物的變動情況卻很方便.我們以30天為例,依次對這兩種方案進(jìn)行仿真,最后比較各方案的總費(fèi)用,從而就可以作出決策.計算機(jī)仿真時的工作流程是早上到貨,全天銷售,晚上訂貨,輸入一些常數(shù)和初始數(shù)據(jù)后,以一天為時間步長進(jìn)行仿真.首先檢查這一天是否為預(yù)訂到貨日期,如果是,則原有庫存量加Q,并把預(yù)定到貨量清為0;如果不是,則庫存量不變.接著用計算機(jī)中的隨機(jī)函數(shù)仿真隨機(jī)需求量,若庫存量大于需求量,則新的庫存量減去需求量;反之,則新庫存量變?yōu)?,并且要在總費(fèi)用上加缺貨損失,然后檢查實(shí)際例9庫存問題庫存量加上預(yù)定到貨量是否小于重新訂貨點(diǎn)P,如果是,則需要重新訂貨,這時就加一次訂貨費(fèi).如此重復(fù)運(yùn)行30天,即可得所需費(fèi)用總值.由此比較這兩種方案的總費(fèi)用,即得最好方案.cleardays=30;P=[125,150];%重新訂貨點(diǎn)P，單位輛Q=[150,250];%重新訂貨量Q，單位輛cost=[0,0];%花費(fèi)初值arrivalinterval=2;%到貨時間間隔為2天storagefee=0.5;lossfee=1.6;bookfee=75;%各種費(fèi)用storage0=110;booknumber=0;arrivedate=0;nr=rand(days,1);%隨機(jī)產(chǎn)生30天的隨機(jī)數(shù)fori=1:2%兩種方案分別仿真

%下面是第一天的情況

storage(1)=storage0;n=round(99*nr(1));%30天中每天的需求量,0-99之間的均勻分布隨機(jī)數(shù)

sale=n;remain=storage(1)-n;%計算庫存

ifremain<=P(i);%第i種方案，庫存少于重新訂貨量，就需要訂貨

booknumber=Q(i);%訂貨數(shù)量

arrivedate=4;%到貨日期

orderfee=bookfee;%訂單費(fèi)用

else

orderfee=0;%未出現(xiàn)缺貨情況，訂單費(fèi)用為0endstorage(1)=remain;%更新貨存

cost(i)=cost(i)+remain*storagefee+orderfee;%總費(fèi)用累加，第一天計算結(jié)束

forj=2:days%第2天到第30天仿真

dh=j;ifabs(dh-arrivedate)<0.01%訂單到貨

storage(j)=storage(j-1)+booknumber;%貨存增加

booknumber=0;%已到貨訂貨數(shù)量歸0arrivedate=j;%到貨日期

else

storage(j)=storage(j-1);%未到貨，當(dāng)前庫存與前一天相等

end%結(jié)束判斷是否訂單到貨

n=round(99*nr(j));%第j天隨機(jī)需求量nifstorage(j)>=n%庫存足夠

sale=n;%成交量

remain=storage(j)-n;%庫存量

shortagenumber=0;%缺貨數(shù)為0else%否則，庫存不夠的處理

sale=storage(j);%只能銷售已有貨存的量，有缺貨

remain=0;%庫存賣完

shortagenumber=n-storage(j);%缺貨量

end

storage(j)=remain;%更新庫存

ifremain+booknumber<=P(i);%庫存量加上預(yù)定到貨量小于重新訂貨點(diǎn)P，重新訂貨

booknumber=Q(i);%訂貨數(shù)

arrivedate=dh+arrivalinterval;%到貨日期

orderfee=bookfee;%訂單費(fèi)用

else

orderfee=0;%無需重新訂貨，未發(fā)生費(fèi)用

end

cost(i)=cost(i)+remain*storagefee+shortagenumber*lossfee+orderfee;%總費(fèi)用增加

end;

mincost=min(cost);%求兩個方案30天總費(fèi)用最小值

endcost

mincost例10趕火車過程仿真

一列火車從A站經(jīng)B站開往C站,某人每天趕往B站乘這趟火車.已知火車從A站到B站運(yùn)行時間為均值30min,標(biāo)準(zhǔn)差為2min的正態(tài)隨機(jī)變量.火車大約在下午1點(diǎn)離開A站,離開時刻的頻率分布見表1.這個人到達(dá)B站時的頻率分布表見表2.用計算機(jī)仿真火車開出,火車到達(dá)B站,這個人到達(dá)B站的情況,并給出能趕上火車的仿真結(jié)果.

表1離開時刻的頻率分布出發(fā)時刻T1:001:051:10頻率0.70.20.1例10趕火車過程仿真

表2人到達(dá)B站時的頻率分布到達(dá)時刻T1:281:301:321:34

頻率:0.30.40.20.