第6講計算機模擬

1、數(shù)學(xué)實驗數(shù)學(xué)實驗 計算機模擬計算機模擬模擬的概念模擬的概念 模擬就是利用物理的、數(shù)學(xué)的模型來類比、模仿現(xiàn)實系統(tǒng)及其演變過程，以尋求過程規(guī)律的一種方法 模擬的基本思想是建立一個試驗的模型，這個模型包含所研究系統(tǒng)的主要特點通過對這個實驗?zāi)Ｐ偷倪\行，獲得所要研究系統(tǒng)的必要信息. 在實際問題中，面對一些帶隨機因素的復(fù)雜系統(tǒng)，用分析方法建模常常需要作許多簡化假設(shè)，與面臨的實際問題可能相差甚遠，以致解答根本無法應(yīng)用這時，計算機模擬幾乎成為唯一的選擇 在一定的假設(shè)條件下，運用數(shù)學(xué)運算模擬系統(tǒng)的運行，稱為數(shù)學(xué)模擬現(xiàn)代的數(shù)學(xué)模擬都是在計算機上進行的，稱為計算機模擬模擬的方法模擬的方法: 數(shù)學(xué)模擬 計算機模擬可以

2、反復(fù)進行，改變系統(tǒng)的結(jié)構(gòu)和系數(shù)都比較容易 蒙特卡羅（蒙特卡羅（Monte CarloMonte Carlo）方法）方法是一種應(yīng)用隨機數(shù)來進行計算機模擬的方法此方法對研究的系統(tǒng)進行隨機觀察抽樣，通過對樣本值的觀察統(tǒng)計，求得所研究系統(tǒng)的某些參數(shù)例例1 1在我方某前沿防守地域，敵人以一個炮排（含兩門火炮）為單位對我方進行干擾和破壞為躲避我方打擊，敵方對其陣地進行了偽裝并經(jīng)常變換射擊地點 經(jīng)過長期觀察發(fā)現(xiàn)，我方指揮所對敵方目標(biāo)的指示有50是準(zhǔn)確的，而我方火力單位，在指示正確時，有1/3的射擊效果能毀傷敵人一門火炮，有1/6的射擊效果能全部消滅敵人 現(xiàn)在希望能用某種方式把我方將要對敵人實施的20次打擊結(jié)

3、果顯現(xiàn)出來，確定有效射擊的比率及毀傷敵方火炮的平均值分析分析：這是一個概率問題，可以通過理論計算得到相應(yīng)的概率和期望值.但這樣只能給出作戰(zhàn)行動的最終靜態(tài)結(jié)果，而顯示不出作戰(zhàn)行動的動態(tài)過程. 為了能顯示我方20次射擊的過程，現(xiàn)采用模擬的方式 需要模擬出以下兩件事： 1. 問題分析問題分析2 2 當(dāng)指示正確時，我方火力單位的射擊結(jié)果情況當(dāng)指示正確時，我方火力單位的射擊結(jié)果情況1 1 觀察所對目標(biāo)的指示正確與否觀察所對目標(biāo)的指示正確與否模擬試驗有兩種結(jié)果，每種結(jié)果出現(xiàn)的概率都是1/2 因此，可用投擲可用投擲1 1枚硬幣的方式予以確定枚硬幣的方式予以確定，當(dāng)硬幣出現(xiàn)正面時為指示正確，反之為不正確 模擬

4、試驗有三種結(jié)果：毀傷1門火炮的可能性為1/3(即2/6)，毀傷兩門的可能性為1/6，沒能毀傷敵火炮的可能性為1/2(即3/6) 這時可用投擲骰子的方法來確定可用投擲骰子的方法來確定：如果出現(xiàn)的是、點：則認為沒能擊中敵人；如果出現(xiàn)的是、點：則認為毀傷敵人一門火炮；若出現(xiàn)的是點：則認為毀傷敵人兩門火炮2. 符號假設(shè)符號假設(shè)i：要模擬的打擊次數(shù)； k1：沒擊中敵人火炮的射擊總數(shù)； k2：擊中敵人一門火炮的射擊總數(shù)；k3：擊中敵人兩門火炮的射擊總數(shù)E：有效射擊比率； E1：20次射擊平均每次毀傷敵人的火炮數(shù)3. 模擬框圖模擬框圖初始化:i=0, k1=0, k2=0, k3=0i=i+1骰子點數(shù)?k1

5、=k1+1k2=k2+1k3=k3+1k1=k1+1i20?E= E1= 0 +1 +2 停止硬幣正面?YNNY1，2，34，56120k220k320k23()20kk投擲硬幣的計算機模擬投擲硬幣的計算機模擬1產(chǎn)生服從U（0，1）的隨機數(shù)R12將區(qū)間0，1等分： 若 ，則對應(yīng)硬幣正面 若 ，則對應(yīng)硬幣反面5 . 001R15 . 01R擲骰子的計算機模擬擲骰子的計算機模擬1產(chǎn)生服從U（0，1）的隨機數(shù)R22將區(qū)間0，1六等份： 若 ，則對應(yīng)骰子點數(shù)為1 若 ，則對應(yīng)骰子點數(shù)為2 若 ，則對應(yīng)骰子點數(shù)為3 若 ，則對應(yīng)骰子點數(shù)為4 若 ，則對應(yīng)骰子點數(shù)為5 若 ，則對應(yīng)骰子點數(shù)為66102 R

6、62612 R63622 R65642 R64632 R1652 R4. 模擬結(jié)果模擬結(jié)果消滅敵人火炮數(shù)消滅敵人火炮數(shù)試驗試驗序號序號投硬幣投硬幣結(jié)果結(jié)果指示指示正確正確指示指示不正確不正確擲骰子擲骰子結(jié)果結(jié)果正正反正正反正正反反消滅敵人火炮數(shù)消滅敵人火炮數(shù)試驗試驗序號序號投硬幣投硬幣結(jié)果結(jié)果指示指示正確正確指示指示不正確不正確擲骰子擲骰子結(jié)果結(jié)果正反正反正正正正反正5. 理論計算理論計算6. 結(jié)果比較結(jié)果比較 返回 雖然模擬結(jié)果與理論計算不完全一致，但它卻能更加真實地表雖然模擬結(jié)果與理論計算不完全一致，但它卻能更加真實地表達實際戰(zhàn)斗動態(tài)過程達實際戰(zhàn)斗動態(tài)過程 用蒙特卡羅方法進行計算機模擬的步

7、驟用蒙特卡羅方法進行計算機模擬的步驟：1 設(shè)計一個邏輯框圖，即模擬模型這個框圖要正確反映系統(tǒng)各部分運行時的邏輯關(guān)系2 模擬隨機現(xiàn)象可通過具有各種概率分布的模擬隨機數(shù)來模擬隨機現(xiàn)象產(chǎn)生模擬隨機數(shù)的計算機命令產(chǎn)生模擬隨機數(shù)的計算機命令 在MATLAB軟件中，可以直接產(chǎn)生滿足各種分布的隨機數(shù)，命令如下：2產(chǎn)生mn階，均勻分布的隨機數(shù)矩陣： rand (m, n) 產(chǎn)生一個，均勻分布的隨機數(shù)：rand1產(chǎn)生mn階a，b上均勻分布U（a，b）的隨機數(shù)矩陣： unifrnd (a,b,m, n) 產(chǎn)生一個a，b均勻分布的隨機數(shù)：unifrnd (a,b)例 1的計算機模擬 To MATLAB（rnd）當(dāng)研

8、究對象視為大量相互獨立的隨機變量之和，且其中每一種變量對總和的影響都很小時，可以認為該對象服從正態(tài)分布機械加工得到的零件尺寸的偏差、射擊命中點與目標(biāo)的偏差、各種測量誤差、人的身高、體重等，都可近似看成服從正態(tài)分布若連續(xù)型隨機變量X的概率密度函數(shù)為 其中 0為常數(shù)，則稱X服從參數(shù)為 的指數(shù)分布指數(shù)分布e0( )00txf xx指數(shù)分布的期望值為 1排隊服務(wù)系統(tǒng)中顧客到達率為常數(shù)時的到達間隔、故障率為常數(shù)時零件的壽命都服從指數(shù)分布指數(shù)分布在排隊論、可靠性分析中有廣泛應(yīng)用注意：注意：MATLAB中，產(chǎn)生參數(shù)為 的指數(shù)分布的命令為exprnd( )1例例 顧客到達某商店的間隔時間服從參數(shù)為顧客到達某商

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

10、系統(tǒng)模擬實例: 追逐問題追逐問題 狀態(tài)隨時間連續(xù)變化的系統(tǒng)稱為連續(xù)系統(tǒng)連續(xù)系統(tǒng)對連續(xù)系統(tǒng)的計算機模擬只能是近似的，只要這種近似達到一定的精度，也就可以滿足要求例例 追逐問題追逐問題: 如圖,正方形ABCD的4個頂點各有1人.在某一時刻,4人同時出發(fā)以勻速v=1m/s按順時針方向追逐下一人,如果他們始終保持對準(zhǔn)目標(biāo),則最終按螺旋狀曲線于中心點O.試求出這種情況下每個人的行進軌跡.1. 建立平面直角坐標(biāo)系: A(x1,y1), B(x2,y2), C(x3,y3), D(x4,y4).2. 取時間間隔為t,計算每一點在各個時刻的坐標(biāo).4. 對每一個點，連接它在各時刻的位置,即得所求運動軌跡.求解過程求解過程: To MATLAB(chase)返回OBCDAv=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(d0.1) x(5)=x(1);y(5)=y(1); for