計算機模擬課件

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

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

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

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

6、5,案例二 隨機模擬拋硬幣問題,硬幣有正反兩面，拋擲的均勻的硬幣，正面和反面朝上的機會都是二分之一，我們利用matlab函數rand產生在（0,1）區(qū)間內均與分布的偽隨機數的功能，設計算法，編寫程序，模擬拋擲的均與的硬幣這種隨機行為。,,算法： 設x是（0，1）區(qū)間分布的隨機數，定義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面，點數分別是1、2、3、4、5、6，如果骰子是質地均勻的，那么投擲骰子擲出每一個點數的機會都是六分之一。我們要設計算法，編寫程序，模擬投擲骰子的這種隨機行為。,算法： 設x是（0，1）區(qū)間分布的隨

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

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

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

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

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

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

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

18、天的庫存量;S1 :方案甲前1 天的銷售量;S21 :方案乙前1 天的銷售量;S22 :方案乙前2 天的銷售量;S3 :方案甲當天實際銷售量;S4 :方案乙當天實際銷售量;L1 :方案甲當天的利潤;L2 :方案乙當天的利潤;TL1 :方案甲累計總利潤;TL2 :方案甲累計總利潤;T :預定模擬天數.,引入下列符號:,模擬過程的流程圖,【模型的求解】,利用Matlab 編程來實現這一過程, 建立如下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、為躲避我方打擊，敵方對其陣地進行了偽裝并經常變換射擊地點．,經過長期觀察發(fā)現，我方指揮所對敵方目標的指示有50％是準確的，而我方火力單位，在指示正確時，有1/3的射擊效果能毀傷敵人一門火炮，有1/6的射擊效果能全部消滅敵人．,現在希望能用某種方式把我方將要對敵人實施的20次打擊結果顯現出來，確定有效射擊的比率及毀傷敵方火炮的平均值．,分析：這是一個概率問題，可以通過理論計算得到相應的概率和期望值.但這樣只能給出作戰(zhàn)行動的最終靜態(tài)結果，而

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

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

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

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

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

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

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