計算機模擬課件

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

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

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

4、/N ，即為點落在圓內(nèi)的概率的近似值；④設(shè)圓的面積為S，由幾何概率公式得點落在內(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程序：,從計算結(jié)果看：這種數(shù)據(jù)模擬算法，在實驗次數(shù)較少時，算的近似值距離真實值誤差較大，但這種數(shù)據(jù)模擬方法簡單易行。在精度要求不高的情況下，這種取隨機數(shù)進行數(shù)據(jù)模擬的方法還是有一定的實用價值，通過簡單編程我們可以模擬出許多數(shù)學(xué)現(xiàn)象。,依次取n=500,1000,100000取得的近似值分別為： 3.0960；3.0680；3.140

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

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

8、機數(shù)，定義0<x<1/6表示擲出1點， 1/6<x<1/3表示擲出2點,….., 5/6<x<1表示擲出6點. 輸入 投骰子的總次數(shù)n 輸出 投出每一個點數(shù)的頻率p（有六個分量） 第1步 初始化計數(shù)器（有六個分量）； 第2步 對i=1,2,….,n,循環(huá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步 計算頻率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)生模擬隨機數(shù)的計算機命令,,在MATLAB軟件中，可以直接產(chǎn)生滿足各種分布的隨機數(shù)，命令如下：,2．產(chǎn)生m×n階［０，１］均勻分布的隨機數(shù)矩陣： rand (m, n) 產(chǎn)生一個［０，１］均勻分布的隨機數(shù)：rand,1．產(chǎn)生m×n階［a，b］上均勻分布U（a，b）的隨機數(shù)矩陣： unifrnd (a,b,m, n) 產(chǎn)生一個［a，b］

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

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

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

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

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

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

18、天的庫存量;S1 :方案甲前1 天的銷售量;S21 :方案乙前1 天的銷售量;S22 :方案乙前2 天的銷售量;S3 :方案甲當天實際銷售量;S4 :方案乙當天實際銷售量;L1 :方案甲當天的利潤;L2 :方案乙當天的利潤;TL1 :方案甲累計總利潤;TL2 :方案甲累計總利潤;T :預(yù)定模擬天數(shù).,引入下列符號:,模擬過程的流程圖,【模型的求解】,利用Matlab 編程來實現(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)模擬實例: 追逐問題,狀態(tài)隨時間連續(xù)變化的系統(tǒng)稱為連續(xù)系統(tǒng)．對連續(xù)系統(tǒng)的計算機模擬只能是近似的，只要這種近似

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

22、刻的位置,即得所求運動軌跡.,,求解過程:,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,計算程序:,To MATLAB(chase),案例六 在我方某前沿防守地域，敵人以一個炮排（含兩門火炮）為單位對我方進行干擾和破壞．

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

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

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

27、均每次毀傷敵人的火炮數(shù)．,3. 模擬框圖,4. 模擬結(jié)果,6. 結(jié)果比較,雖然模擬結(jié)果與理論計算不完全一致，但它卻能更加真實地表達實際戰(zhàn)斗動態(tài)過程．,用蒙特卡羅方法進行計算機模擬的步驟：,[1] 設(shè)計一個邏輯框圖，即模擬模型．這個框圖要正確反映系統(tǒng)各部分運行時的邏輯關(guān)系．[2] 模擬隨機現(xiàn)象．可通過具有各種概率分布的模擬隨機數(shù)來模擬隨機現(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)模擬實例: 排隊問題,排隊論主要研究隨機服務(wù)系統(tǒng)的工作過程．,在排隊系統(tǒng)中，服務(wù)對象的到達時間和服務(wù)時間都是隨機的．排隊論通過對每個個別的隨機服務(wù)現(xiàn)象的統(tǒng)計研究

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

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

31、務(wù)的個數(shù)及每日顧客的平均等待時間．,[2] 符號說明 w：總等待時間；ci：第i個顧客的到達時刻；　　 bi：第i個顧客開始服務(wù)時刻； ei：第i個顧客服務(wù)結(jié)束時刻． xi：第i-1個顧客與第i個顧客到達之間的時間間隔 yi：對第i個顧客的服務(wù)時間,,,,,,,,,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編一個福利彩票雙色球電腦選號的程序并要求前六個數(shù)輸出結(jié)果按升序排序．,3. 某設(shè)備上安裝有4只型號規(guī)格完全相同的電子管，已知電子管壽命服從100～2

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