計(jì)算機(jī)模擬課件

1、數(shù)學(xué)建模培訓(xùn),———計(jì)算機(jī)模擬,模擬的概念,模擬就是利用物理的、數(shù)學(xué)的模型來類比、模仿現(xiàn)實(shí)系統(tǒng)及其演變過程，以尋求過程規(guī)律的一種方法．,模擬的基本思想是建立一個(gè)試驗(yàn)的模型，這個(gè)模型包含所研究系統(tǒng)的主要特點(diǎn)．通過對這個(gè)實(shí)驗(yàn)?zāi)Ｐ偷倪\(yùn)行，獲得所要研究系統(tǒng)的必要信息.,對實(shí)際系統(tǒng)及其過程用功能相似的實(shí)物系統(tǒng)去模仿．例如，軍事演習(xí)、船艇實(shí)驗(yàn)、沙盤作業(yè)等．,物理模擬通?；ㄙM(fèi)較大、周期較長，且在物理模型上改變系統(tǒng)結(jié)構(gòu)和系數(shù)都較困難．而且，許多系統(tǒng)無法

2、進(jìn)行物理模擬，如社會(huì)經(jīng)濟(jì)系統(tǒng)、生態(tài)系統(tǒng)等．,模擬的方法,物理模擬,,,在一定的假設(shè)條件下，運(yùn)用數(shù)學(xué)運(yùn)算模擬系統(tǒng)的運(yùn)行，稱為數(shù)學(xué)模擬．現(xiàn)代的數(shù)學(xué)模擬都是在計(jì)算機(jī)上進(jìn)行的，稱為計(jì)算機(jī)模擬．,數(shù)學(xué)模擬,計(jì)算機(jī)模擬可以反復(fù)進(jìn)行，改變系統(tǒng)的結(jié)構(gòu)和系數(shù)都比較容易．,在實(shí)際問題中，面對一些帶隨機(jī)因素的復(fù)雜系統(tǒng)，用分析方法建模常常需要作許多簡化假設(shè)，與面臨的實(shí)際問題可能相差甚遠(yuǎn)，以致解答根本無法應(yīng)用．這時(shí)，計(jì)算機(jī)模擬幾乎成為唯一的選擇．計(jì)算機(jī)模擬與數(shù)學(xué)

3、建模的關(guān)系如圖所示：,案例一 圓周率的近似計(jì)算 ——蒙特卡羅法,在正方形中隨機(jī)撒一把豆子，通過考察落在其內(nèi)切圓內(nèi)黃豆的數(shù)目，用隨機(jī)模擬的方法計(jì)算圓周率P的近似值（如圖）,具體方案如下：①利用計(jì)算機(jī)產(chǎn)生兩組[0，1]上的均勻隨機(jī)數(shù)，x=RAND，y=RAND；②統(tǒng)計(jì)試驗(yàn)總次數(shù)N和落在內(nèi)切圓內(nèi)的點(diǎn)數(shù)N1（滿足條件x2+y2≤1的點(diǎn)（x，y）的個(gè)數(shù)）；③計(jì)算頻率N1

4、/N ，即為點(diǎn)落在圓內(nèi)的概率的近似值；④設(shè)圓的面積為S，由幾何概率公式得點(diǎn)落在內(nèi)切圓部分的概率為P=S/ 4 ．∴S/ 4 = N1 /N ．∴S≈4N1/N ，即為圓的面積的近似值． 又S圓=πr2=π，∴π=S≈4N1 /N ，即為圓周率的近似值．,模擬過程流程圖,cs=0n=500for i=1:n a=rand(1,2); if a(1)^2+a(2)^2<=1 cs=cs

5、+1 endend4*cs/n,Matlab程序：,從計(jì)算結(jié)果看：這種數(shù)據(jù)模擬算法，在實(shí)驗(yàn)次數(shù)較少時(shí)，算的近似值距離真實(shí)值誤差較大，但這種數(shù)據(jù)模擬方法簡單易行。在精度要求不高的情況下，這種取隨機(jī)數(shù)進(jìn)行數(shù)據(jù)模擬的方法還是有一定的實(shí)用價(jià)值，通過簡單編程我們可以模擬出許多數(shù)學(xué)現(xiàn)象。,依次取n=500,1000,100000取得的近似值分別為： 3.0960；3.0680；3.140

6、5,案例二 隨機(jī)模擬拋硬幣問題,硬幣有正反兩面，拋擲的均勻的硬幣，正面和反面朝上的機(jī)會(huì)都是二分之一，我們利用matlab函數(shù)rand產(chǎn)生在（0,1）區(qū)間內(nèi)均與分布的偽隨機(jī)數(shù)的功能，設(shè)計(jì)算法，編寫程序，模擬拋擲的均與的硬幣這種隨機(jī)行為。,,算法： 設(shè)x是（0，1）區(qū)間分布的隨機(jī)數(shù)，定義x>=0.5表示正面朝上，x=0.5，k=k+1;否則k不變； 第5步 計(jì)算頻率p=k/n,拋硬幣的matlab程序,k=

7、0n=10for i=1:n x=rand; if x>=0.5 k=k+1; endendkp=k/n,案例三 隨機(jī)模擬擲骰子問題,骰子呈正方體，有6面，點(diǎn)數(shù)分別是1、2、3、4、5、6，如果骰子是質(zhì)地均勻的，那么投擲骰子擲出每一個(gè)點(diǎn)數(shù)的機(jī)會(huì)都是六分之一。我們要設(shè)計(jì)算法，編寫程序，模擬投擲骰子的這種隨機(jī)行為。,算法： 設(shè)x是（0，1）區(qū)間分布的隨

8、機(jī)數(shù)，定義0<x<1/6表示擲出1點(diǎn)， 1/6<x<1/3表示擲出2點(diǎn),….., 5/6<x<1表示擲出6點(diǎn). 輸入 投骰子的總次數(shù)n 輸出 投出每一個(gè)點(diǎn)數(shù)的頻率p（有六個(gè)分量） 第1步 初始化計(jì)數(shù)器（有六個(gè)分量）； 第2步 對i=1,2,….,n,循環(huán)進(jìn)行第3,4步； 第3步 x=rand; 第4步 如果0<x<

9、1/6 ，k(1)=k(1)+1; 如果1/6<x<1/3，k(2)=k(2)+1; ……. 如果5/6<x<1，k(6)=k(6)+1; 第5步 計(jì)算頻率p=k/n,n=100;k=zeros(1,6);for i=1:n x=rand; if x<1/6

10、 k(1)=k(1)+1; elseif x<2/6 k(2)=k(2)+1; elseif x<3/6 k(3)=k(3)+1; elseif x<4/6 k(4)=k(4)+1; elseif x<5/6 k(5)=k(5)+1; else k(6)=k(6)+1; endend

11、kp=k/n,投擲骰子的matlab程序,產(chǎn)生模擬隨機(jī)數(shù)的計(jì)算機(jī)命令,,在MATLAB軟件中，可以直接產(chǎn)生滿足各種分布的隨機(jī)數(shù)，命令如下：,2．產(chǎn)生m×n階［０，１］均勻分布的隨機(jī)數(shù)矩陣： rand (m, n) 產(chǎn)生一個(gè)［０，１］均勻分布的隨機(jī)數(shù)：rand,1．產(chǎn)生m×n階［a，b］上均勻分布U（a，b）的隨機(jī)數(shù)矩陣： unifrnd (a,b,m, n) 產(chǎn)生一個(gè)［a，b］

12、均勻分布的隨機(jī)數(shù)：unifrnd (a,b),當(dāng)只知道一個(gè)隨機(jī)變量取值在（a，b）內(nèi)，但不知道（也沒理由假設(shè)）它在何處取值的概率大，在何處取值的概率小，就只好用U（a，b）來模擬它．,例 1的計(jì)算機(jī)模擬,To MATLAB（rnd）,當(dāng)研究對象視為大量相互獨(dú)立的隨機(jī)變量之和，且其中每一種變量對總和的影響都很小時(shí)，可以認(rèn)為該對象服從正態(tài)分布．,機(jī)械加工得到的零件尺寸的偏差、射擊命中點(diǎn)與目標(biāo)的偏差、各種測量誤差、人的身高、體重等，都可近似

13、看成服從正態(tài)分布．,,若連續(xù)型隨機(jī)變量X的概率密度函數(shù)為 其中 >0為常數(shù)，則稱X服從參數(shù)為 的指數(shù)分布．,指數(shù)分布的期望值為,指數(shù)分布可以用來表示獨(dú)立隨機(jī)事件發(fā)生的時(shí)間間隔，比如旅客進(jìn)機(jī)場的時(shí)間間隔，排隊(duì)服務(wù)系統(tǒng)中顧客到達(dá)率為常數(shù)時(shí)的到達(dá)間隔、故障率為常數(shù)時(shí)零件的壽命都服從指數(shù)分布．,指數(shù)分布在排隊(duì)論、可靠性分析中有廣泛應(yīng)用．,注意：MATLAB中，產(chǎn)生參數(shù)為 的指數(shù)分布的命令為exprnd( ),例 顧客到達(dá)某商

14、店的間隔時(shí)間服從參數(shù)為0.1的指數(shù)分布,指數(shù)分布的均值為1/0.1=10． 指兩個(gè)顧客到達(dá)商店的平均間隔時(shí)間是10個(gè)單位時(shí)間.即平均10個(gè)單位時(shí)間到達(dá)1個(gè)顧客. 顧客到達(dá)的間隔時(shí)間可用exprnd(10)模擬．,設(shè)離散型隨機(jī)變量X的所有可能取值為0,1,2,…,且取各個(gè)值的概率為其中 >0為常數(shù)，則稱X服從參數(shù)為 的泊松分布．,泊松分布在排隊(duì)系統(tǒng)、產(chǎn)品檢驗(yàn)、天文、物理等領(lǐng)域有廣泛應(yīng)用．,泊松分布的期望值為,案例四　敵坦克

15、分隊(duì)對我方陣地實(shí)施突襲，其到達(dá)規(guī)律服從泊松分布，平均每分鐘到達(dá)４輛．（1）模擬敵坦克在３分鐘內(nèi)到達(dá)目標(biāo)區(qū)的數(shù)量，以及在第１、２、３分鐘內(nèi)各到達(dá)幾輛坦克．（2）模擬在3分鐘內(nèi)每輛敵坦克的到達(dá)時(shí)刻．,（1）用poissrnd(4)進(jìn)行模擬．,To MATLAB（poiss）,（2）坦克到達(dá)的間隔時(shí)間應(yīng)服從參數(shù)為4的負(fù)指數(shù)分布，用exprnd（1/4）模擬．,To MATLAB（time）,案例四 企業(yè)庫存系統(tǒng)的計(jì)算機(jī)模擬方法,在銷售部門、工

16、廠等領(lǐng)域中都存在庫存問題, 庫存太多造成浪費(fèi)以及資金積壓,庫存少了不能滿足需求也造成損失. 部門的工作人員需決定何時(shí)進(jìn)貨,進(jìn)多少,使得所花費(fèi)的平均費(fèi)用最少,而收益最大,這就是庫存問題. 某企業(yè)生產(chǎn)易變質(zhì)的產(chǎn)品. 當(dāng)天生產(chǎn)的產(chǎn)品必須售出,否則就會(huì)變質(zhì). 該產(chǎn)品單位成本為2.5 元,單位產(chǎn)品售價(jià)為5 元. 企業(yè)為避免存貨過多而造成損失,擬從以下2 種庫存方案中選出一個(gè)較優(yōu)的方案: 方案甲:按前1 天的銷售量作為當(dāng)天的庫存量

17、; 方案乙:按前2 天的平均銷售量作為當(dāng)天的庫存量.,【問題的提出】,【模型建立】,假定市場對該產(chǎn)品的每天需求量是一個(gè)隨機(jī)變量,但從以往的統(tǒng)計(jì)分析得知它服從正態(tài)分布: N(135 , 22. 4) .,第一：獲得市場對該產(chǎn)品的需求量的數(shù)據(jù);,第二：計(jì)算出按照2 種不同方案經(jīng)T 天后企業(yè)的利潤值；,第三：比較大小,從中選出一個(gè)更優(yōu)的方案.,計(jì)算機(jī)模擬的基本思路:,D :每天需求量;Q1 :方案甲當(dāng)天的庫存量;Q2 :方案甲當(dāng)

18、天的庫存量;S1 :方案甲前1 天的銷售量;S21 :方案乙前1 天的銷售量;S22 :方案乙前2 天的銷售量;S3 :方案甲當(dāng)天實(shí)際銷售量;S4 :方案乙當(dāng)天實(shí)際銷售量;L1 :方案甲當(dāng)天的利潤;L2 :方案乙當(dāng)天的利潤;TL1 :方案甲累計(jì)總利潤;TL2 :方案甲累計(jì)總利潤;T :預(yù)定模擬天數(shù).,引入下列符號:,模擬過程的流程圖,【模型的求解】,利用Matlab 編程來實(shí)現(xiàn)這一過程, 建立如下M文件:,funct

19、ion[ TL1 , TL2 ] = kucun( T , S1 , S21 , S22)TL1 = 0; TL2 = 0 ; k = 1;while k < T　Q1 = S1 ; Q2 = ( S21 + S22) / 2 ;　D = normrnd(135 ,22. 4) ;　if D < Q1　　S3 = Q1 ;　else　　S3 = D;　end　if D < Q2　　S4 = Q2

20、;　else　　S4 = D;　end　L1 = 5 *S3 - 2. 5 * Q1 ; L2 = 5 *S4 - 2. 5 *Q2 ;　TL1 = TL1 + L1 ; TL2 = TL2 + L2 ;　k = k + 1 ;endS1 = S3 ; S22 = S21 ; S21 = S4 ;,案例五連續(xù)系統(tǒng)模擬實(shí)例: 追逐問題,狀態(tài)隨時(shí)間連續(xù)變化的系統(tǒng)稱為連續(xù)系統(tǒng)．對連續(xù)系統(tǒng)的計(jì)算機(jī)模擬只能是近似的，只要這種近似

21、達(dá)到一定的精度，也就可以滿足要求．,追逐問題: 如圖,正方形ABCD的4個(gè)頂點(diǎn)各有1人.在某一時(shí)刻,4人同時(shí)出發(fā)以勻速v=1m/s按順時(shí)針方向追逐下一人,如果他們始終保持對準(zhǔn)目標(biāo),則最終按螺旋狀曲線于中心點(diǎn)O.試求出這種情況下每個(gè)人的行進(jìn)軌跡.,1. 建立平面直角坐標(biāo)系: A(x1,y1), B(x2,y2), C(x3,y3), D(x4,y4).,2. 取時(shí)間間隔為Δt,計(jì)算每一點(diǎn)在各個(gè)時(shí)刻的坐標(biāo).,4. 對每一個(gè)點(diǎn)，連接它在各時(shí)

22、刻的位置,即得所求運(yùn)動(dòng)軌跡.,,求解過程:,To MATLAB(chase),v=1;dt=0.05;x=[0 0 10 10];y=[0 10 10 0];for i=1:4 plot(x(i),y(i),'.'),hold onendd=20;while(d>0.1) x(5)=x(1);y(5)=y(1); for i=1:4 d=sqrt((x(i+1)-x(

23、i))^2+(y(i+1)-y(i))^2); x(i)=x(i)+v*dt*(x(i+1)-x(i))/d; y(i)=y(i)+v*dt*(y(i+1)-y(i))/d; plot(x(i),y(i),'.'),hold on end end,計(jì)算程序:,To MATLAB(chase),案例六 在我方某前沿防守地域，敵人以一個(gè)炮排（含兩門火炮）為單位對我方進(jìn)行干擾和破壞．

24、為躲避我方打擊，敵方對其陣地進(jìn)行了偽裝并經(jīng)常變換射擊地點(diǎn)．,經(jīng)過長期觀察發(fā)現(xiàn)，我方指揮所對敵方目標(biāo)的指示有50％是準(zhǔn)確的，而我方火力單位，在指示正確時(shí)，有1/3的射擊效果能毀傷敵人一門火炮，有1/6的射擊效果能全部消滅敵人．,現(xiàn)在希望能用某種方式把我方將要對敵人實(shí)施的20次打擊結(jié)果顯現(xiàn)出來，確定有效射擊的比率及毀傷敵方火炮的平均值．,分析：這是一個(gè)概率問題，可以通過理論計(jì)算得到相應(yīng)的概率和期望值.但這樣只能給出作戰(zhàn)行動(dòng)的最終靜態(tài)結(jié)果，而

25、顯示不出作戰(zhàn)行動(dòng)的動(dòng)態(tài)過程.,為了能顯示我方20次射擊的過程，現(xiàn)采用模擬的方式．,需要模擬出以下兩件事：,1. 問題分析,[2] 當(dāng)指示正確時(shí)，我方火力單位的射擊結(jié)果情況,[1] 觀察所對目標(biāo)的指示正確與否,模擬試驗(yàn)有兩種結(jié)果，每種結(jié)果出現(xiàn)的概率都是1/2．,因此，可用投擲1枚硬幣的方式予以確定，當(dāng)硬幣出現(xiàn)正面時(shí)為指示正確，反之為不正確．,模擬試驗(yàn)有三種結(jié)果：毀傷1門火炮的可能性為1/3(即2/6)，毀傷兩門的可能性為1/6，沒能毀傷

26、敵火炮的可能性為1/2(即3/6)．,這時(shí)可用投擲骰子的方法來確定：如果出現(xiàn)的是１、２、３點(diǎn)：則認(rèn)為沒能擊中敵人；如果出現(xiàn)的是４、５點(diǎn)：則認(rèn)為毀傷敵人一門火炮；若出現(xiàn)的是６點(diǎn)：則認(rèn)為毀傷敵人兩門火炮．,2. 符號假設(shè),i：要模擬的打擊次數(shù)； k1：沒擊中敵人火炮的射擊總數(shù)； k2：擊中敵人一門火炮的射擊總數(shù)；k3：擊中敵人兩門火炮的射擊總數(shù)．E：有效射擊比率； E1：20次射擊平

27、均每次毀傷敵人的火炮數(shù)．,3. 模擬框圖,4. 模擬結(jié)果,6. 結(jié)果比較,雖然模擬結(jié)果與理論計(jì)算不完全一致，但它卻能更加真實(shí)地表達(dá)實(shí)際戰(zhàn)斗動(dòng)態(tài)過程．,用蒙特卡羅方法進(jìn)行計(jì)算機(jī)模擬的步驟：,[1] 設(shè)計(jì)一個(gè)邏輯框圖，即模擬模型．這個(gè)框圖要正確反映系統(tǒng)各部分運(yùn)行時(shí)的邏輯關(guān)系．[2] 模擬隨機(jī)現(xiàn)象．可通過具有各種概率分布的模擬隨機(jī)數(shù)來模擬隨機(jī)現(xiàn)象．,,cleark1=0;k2=0;k3=0;for i=1:20 R1=rand

28、if R15/6 k3=k3+1; else k2=k2+1; end endelse k1=k1+1;endendiE=(k2+k3)/20E1=(0*k1+1*k2+2*k3)/20,離散系統(tǒng)模擬實(shí)例: 排隊(duì)問題,排隊(duì)論主要研究隨機(jī)服務(wù)系統(tǒng)的工作過程．,在排隊(duì)系統(tǒng)中，服務(wù)對象的到達(dá)時(shí)間和服務(wù)時(shí)間都是隨機(jī)的．排隊(duì)論通過對每個(gè)個(gè)別的隨機(jī)服務(wù)現(xiàn)象的統(tǒng)計(jì)研究

29、，找出反映這些隨機(jī)現(xiàn)象平均特性的規(guī)律，從而為設(shè)計(jì)新的服務(wù)系統(tǒng)和改進(jìn)現(xiàn)有服務(wù)系統(tǒng)的工作提供依據(jù)．,對于排隊(duì)服務(wù)系統(tǒng), 顧客常常注意排隊(duì)的人是否太多, 等候的時(shí)間是否長, 而服務(wù)員則關(guān)心他空閑的時(shí)間是否太短. 于是人們常用排隊(duì)的長度、等待的時(shí)間及服務(wù)利用率等指標(biāo)來衡量系統(tǒng)的性能.,[1] 系統(tǒng)的假設(shè)：（1） 顧客源是無窮的； （2） 排隊(duì)的長度沒有限制；（ 3） 到達(dá)系統(tǒng)的顧客按先后順序依次進(jìn)入服務(wù)， 即“先到先服務(wù)”．,單服務(wù)員的排

30、隊(duì)模型：在某商店有一個(gè)售貨員，顧客陸續(xù)來到，售貨員逐個(gè)地接待顧客．當(dāng)?shù)絹淼念櫩洼^多時(shí)，一部分顧客便須排隊(duì)等待，被接待后的顧客便離開商店．設(shè)： 1．顧客到來間隔時(shí)間服從參數(shù)為0.1的指數(shù)分布．２．對顧客的服務(wù)時(shí)間服從［４，１５］上的均勻分布．３．排隊(duì)按先到先服務(wù)規(guī)則，隊(duì)長無限制．,假定一個(gè)工作日為8小時(shí)，時(shí)間以分鐘為單位．[1]模擬一個(gè)工作日內(nèi)完成服務(wù)的個(gè)數(shù)及顧客平均等待時(shí)間t．[2]模擬100個(gè)工作日，求出平均每日完成服

31、務(wù)的個(gè)數(shù)及每日顧客的平均等待時(shí)間．,[2] 符號說明 w：總等待時(shí)間；ci：第i個(gè)顧客的到達(dá)時(shí)刻；　　 bi：第i個(gè)顧客開始服務(wù)時(shí)刻； ei：第i個(gè)顧客服務(wù)結(jié)束時(shí)刻． xi：第i-1個(gè)顧客與第i個(gè)顧客到達(dá)之間的時(shí)間間隔 yi：對第i個(gè)顧客的服務(wù)時(shí)間,,,,,,,,,c1,,b1,,,c3,c4,c5,c2,e1,b2,e2,b3,e3,b4,,e4,b5,ci=ci-1+

32、xiei=bi+yibi=max(ci,ei-1),t,[1]模擬1日 To MATLAB(simu1),[2]模擬100日To MATLAB(simu2),,cleari=2;w=0;x(i)=exprnd(10);c(i)=x(i);b(i)=x(i);while b(i)<=480 y(i)=unifrnd(4,15); e(i)=b(i)+y(i); w=w+b(i)-c(

33、i); i=i+1; x(i)=exprnd(10); c(i)=c(i-1)+x(i); b(i)=max(c(i),e(i-1)); endi=i-2;t=w/im=i,,作 業(yè),1．用matlab編一個(gè)福利彩票雙色球電腦選號的程序并要求前六個(gè)數(shù)輸出結(jié)果按升序排序．,3. 某設(shè)備上安裝有4只型號規(guī)格完全相同的電子管，已知電子管壽命服從100～2

34、00h之間的均勻分布．電子管損壞時(shí)有兩種維修方案，一是每次更換損壞的那只；二是當(dāng)其中1只損壞時(shí)4只同時(shí)更換．已知更換時(shí)間為換1只時(shí)需1h，4只同時(shí)換為2h．更換時(shí)機(jī)器因停止運(yùn)轉(zhuǎn)每小時(shí)的損失為20元，又每只電子管價(jià)格10元，試用模擬方法確定哪一個(gè)方案經(jīng)濟(jì)合理？,4. 設(shè)位于坐標(biāo)原點(diǎn)的甲艦向位于x軸上點(diǎn)A(1, 0)處的乙艦發(fā)射導(dǎo)彈，導(dǎo)彈頭始終對準(zhǔn)乙艦.如果乙艦以最大的速度(常數(shù))沿平行于y軸的直線行駛，導(dǎo)彈的速度為5. 模擬導(dǎo)彈運(yùn)行的軌跡