1數學建模與數學實驗

計算機模擬2實驗目的實驗內容學習計算機模擬的根本過程與方法．1．模擬的概念．4．實驗作業．3．計算機模擬實例．2．產生隨機數的計算機命令．3連續系統模擬實例:追逐問題離散系統模擬實例:排隊問題用蒙特卡羅法解非線性規劃問題返回計算機模擬實例4模擬的概念

模擬就是利用物理的、數學的模型來類比、模仿現實系統及其演變過程，以尋求過程規律的一種方法．

模擬的根本思想是建立一個試驗的模型，這個模型包含所研究系統的主要特點．通過對這個實驗模型的運行，獲得所要研究系統的必要信息.5模擬的方法1．物理模擬：對實際系統及其過程用功能相似的實物系統去模仿．例如，軍事演習、船艇實驗、沙盤作業等．

物理模擬通?；ㄙM較大、周期較長，且在物理模型上改變系統結構和系數都較困難．而且，許多系統無法進行物理模擬，如社會經濟系統、生態系統等．6

在實際問題中，面對一些帶隨機因素的復雜系統，用分析方法建模常常需要作許多簡化假設，與面臨的實際問題可能相差甚遠，以致解答根本無法應用．這時，計算機模擬幾乎成為唯一的選擇．

在一定的假設條件下，運用數學運算模擬系統的運行，稱為數學模擬．現代的數學模擬都是在計算機上進行的，稱為計算機模擬．2．數學模擬

計算機模擬可以反復進行，改變系統的結構和系數都比較容易．

蒙特卡羅〔MonteCarlo〕方法是一種應用隨機數來進行計算機模擬的方法．此方法對研究的系統進行隨機觀察抽樣，通過對樣本值的觀察統計，求得所研究系統的某些參數．7例1在我方某前沿防守地域，敵人以一個炮排〔含兩門火炮〕為單位對我方進行干擾和破壞．為躲避我方打擊，敵方對其陣地進行了偽裝并經常變換射擊地點．

經過長期觀察發現，我方指揮所對敵方目標的指示有50％是準確的，而我方火力單位，在指示正確時，有1/3的射擊效果能毀傷敵人一門火炮，有1/6的射擊效果能全部消滅敵人．

現在希望能用某種方式把我方將要對敵人實施的20次打擊結果顯現出來，確定有效射擊的比率及毀傷敵方火炮的平均值．分析：這是一個概率問題，可以通過理論計算得到相應的概率和期望值.但這樣只能給出作戰行動的最終靜態結果，而顯示不出作戰行動的動態過程.

為了能顯示我方20次射擊的過程，現采用模擬的方式．8

需要模擬出以下兩件事：

1.問題分析[2]當指示正確時，我方火力單位的射擊結果情況[1]觀察所對目標的指示正確與否模擬試驗有兩種結果，每種結果出現的概率都是1/2．

因此，可用投擲1枚硬幣的方式予以確定，當硬幣出現正面時為指示正確，反之為不正確．

模擬試驗有三種結果：毀傷1門火炮的可能性為1/3(即2/6)，毀傷兩門的可能性為1/6，沒能毀傷敵火炮的可能性為1/2(即3/6)．

這時可用投擲骰子的方法來確定：如果出現的是１、２、３點：那么認為沒能擊中敵人；如果出現的是４、５點：那么認為毀傷敵人一門火炮；假設出現的是６點：那么認為毀傷敵人兩門火炮．92.符號假設i：要模擬的打擊次數；k1：沒擊中敵人火炮的射擊總數；k2：擊中敵人一門火炮的射擊總數；k3：擊中敵人兩門火炮的射擊總數．E：有效射擊比率；E1：20次射擊平均每次毀傷敵人的火炮數．3.模擬框圖初始化:i=0,k1=0,k2=0,k3=0i=i+1骰子點數?k1=k1+1k2=k2+1k3=k3+1k1=k1+1i＜20?E=E1=0×

+1×

+2×

停止硬幣正面?YNNY1，2，34，56104.模擬結果11125.理論計算136.結果比較

返回

雖然模擬結果與理論計算不完全一致，但它卻能更加真實地表達實際戰斗動態過程．用蒙特卡羅方法進行計算機模擬的步驟：[1]設計一個邏輯框圖，即模擬模型．這個框圖要正確反映系統各局部運行時的邏輯關系．[2]模擬隨機現象．可通過具有各種概率分布的模擬隨機數來模擬隨機現象．14產生模擬隨機數的計算機命令

在MATLAB軟件中，可以直接產生滿足各種分布的隨機數，命令如下：2．產生m×n階［０，１］均勻分布的隨機數矩陣：

rand(m,n)

產生一個［０，１］均勻分布的隨機數：rand1．產生m×n階［a，b］上均勻分布U〔a，b〕的隨機數矩陣：

unifrnd(a,b,m,n)

產生一個［a，b］均勻分布的隨機數：unifrnd(a,b)

當只知道一個隨機變量取值在〔a，b〕內，但不知道〔也沒理由假設〕它在何處取值的概率大，在何處取值的概率小，就只好用U〔a，b〕來模擬它．例1的計算機模擬15當研究對象視為大量相互獨立的隨機變量之和，且其中每一種變量對總和的影響都很小時，可以認為該對象服從正態分布．機械加工得到的零件尺寸的偏差、射擊命中點與目標的偏差、各種測量誤差、人的身高、體重等，都可近似看成服從正態分布．16若連續型隨機變量X的概率密度函數為其中>0為常數，則稱X服從參數為的指數分布．指數分布的期望值為

排隊效勞系統中顧客到達率為常數時的到達間隔、故障率為常數時零件的壽命都服從指數分布．指數分布在排隊論、可靠性分析中有廣泛應用．注意：MATLAB中，產生參數為的指數分布的命令為exprnd()例

指數分布的均值為1/0.1=10．指兩個顧客到達商店的平均間隔時間是10個單位時間.即平均10個單位時間到達1個顧客.顧客到達的間隔時間可用exprnd(10)模擬．17設離散型隨機變量X的所有可能取值為0,1,2,…,且取各個值的概率為其中>0為常數，則稱X服從參數為的泊松分布．泊松分布在排隊系統、產品檢驗、天文、物理等領域有廣泛應用．泊松分布的期望值為18如相繼兩個事件出現的間隔時間服從參數為的指數分布，則在單位時間間隔內事件出現的次數服從參數為的泊松分布．即單位時間內該事件出現k次的概率為：反之亦然．指數分布與泊松分布的關系：(1)指兩個顧客到達商店的平均間隔時間是10個單位時間.即平均10個單位時間到達1個顧客.例(1)顧客到達某商店的間隔時間服從參數為0.1的指數分布

(2)該商店在單位時間內到達的顧客數服從參數為0.1的泊松分布19

返回例2敵坦克分隊對我方陣地實施突襲，其到達規律服從泊松分布，平均每分鐘到達４輛．〔1〕模擬敵坦克在３分鐘內到達目標區的數量，以及在第１、２、３分鐘內各到達幾輛坦克．〔2〕模擬在3分鐘內每輛敵坦克的到達時刻．〔1〕用poissrnd(4)進行模擬．〔2〕坦克到達的間隔時間應服從參數為4的負指數分布，用exprnd〔1/4〕模擬．ToMATLAB〔time〕20連續系統模擬實例:追逐問題

狀態隨時間連續變化的系統稱為連續系統．對連續系統的計算機模擬只能是近似的，只要這種近似到達一定的精度，也就可以滿足要求．例追逐問題:如圖,正方形ABCD的4個頂點各有1人.在某一時刻,4人同時出發以勻速v=1m/s按順時針方向追逐下一人,如果他們始終保持對準目標,那么最終按螺旋狀曲線于中心點O.試求出這種情況下每個人的行進軌跡.OBCDA211.建立平面直角坐標系:A(x1,y1),B(x2,y2),C(x3,y3),D(x4,y4).2.取時間間隔為Δt,計算每一點在各個時刻的坐標.4.對每一個點，連接它在各時刻的位置,即得所求運動軌跡.求解過程:返回22v=1;dt=0.05;x=[001010];y=[010100];fori=1:4plot(x(i),y(i),'.'),holdonendd=20;while(d>0.1)x(5)=x(1);y(5)=y(1);

fori=1:4d=sqrt((x(i+1)-x(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),'.'),holdon

endend計算程序:返回23離散系統模擬實例:排隊問題

排隊論主要研究隨機效勞系統的工作過程．

在排隊系統中，效勞對象的到達時間和效勞時間都是隨機的．排隊論通過對每個個別的隨機效勞現象的統計研究，找出反映這些隨機現象平均特性的規律，從而為設計新的效勞系統和改進現有效勞系統的工作提供依據．

對于排隊效勞系統,顧客常常注意排隊的人是否太多,等候的時間是否長,而效勞員那么關心他空閑的時間是否太短.于是人們常用排隊的長度、等待的時間及效勞利用率等指標來衡量系統的性能.24[1]系統的假設：〔1〕顧客源是無窮的；

〔2〕排隊的長度沒有限制；〔3〕到達系統的顧客按先后順序依次進入效勞，即“先到先效勞〞．單效勞員的排隊模型：在某商店有一個售貨員，顧客陸續來到，售貨員逐個地接待顧客．當到來的顧客較多時，一局部顧客便須排隊等待，被接待后的顧客便離開商店．設：1．顧客到來間隔時間服從參數為0.1的指數分布．２．對顧客的效勞時間服從［４，１５］上的均勻分布．３．排隊按先到先效勞規那么，隊長無限制．

假定一個工作日為8小時，時間以分鐘為單位．[1]模擬一個工作日內完成效勞的個數及顧客平均等待時間t．[2]模擬100個工作日，求出平均每日完成效勞的個數及每日顧客的平均等待時間．25[2]符號說明

w：總等待時間；ci：第i個顧客的到達時刻；

bi：第i個顧客開始效勞時刻；ei：第i個顧客效勞結束時刻．

xi：第i-1個顧客與第i個顧客到達之間的時間間隔

yi：對第i個顧客的效勞時間c1b1c3c4c5c2e1b2e2b3e3b4e4b5ci=ci-1+xiei=bi+yibi=max(ci,ei-1)t26[3]模擬框圖初始化：令i=1，ei-1=0，w=0產生間隔時間隨機數xi服從參數為0.1的指數分布

ci=xi

,bi=xi

產生服務時間隨機數yi服從[4,15]的均勻分布

ei=bi+yi累計等待時間：w=w+bi-ci準備下一次服務：i=i+1產生間隔時間隨機數xi服從參數為0.1的指數分布

ci=ci-1+xi

確定開始服務時間：bi=max(ci,ei-1)bi＞480?YNi=i-1,t=w/i輸出結果：完成服務個數：m=i

平均等待時間：t停止[1]模擬1日[2]模擬100日ToMATLAB(simu2)返回27用蒙特卡羅法解非線性規劃問題28根本假設

試驗點的第j個分量xj服從[aj,bj]內的均勻分布．符號假設求解過程

先產生一個隨機數作為初始試驗點,以后那么將上一個試驗點的第j個分量隨機產生，其他分量不變而產生一新的試驗點．這樣，每產生一個新試驗點只需一個新的隨機數分量．當K>MAXK或P>MAXP時停止迭代．29框圖初始化:給定MAXK,MAXP;k=0,p=0,Q:大整數xj=aj+R(bj-aj)j=1,2,…,nj=0j=j+1,p=p+1P>MAXP?YNxj=aj+R(bj-aj)gi(X)≥0?i=1,2…nYNj<n?f(X)≥Q?YNX*=X,Q=f(X)k=k+1k>MAXK?YN輸出X,Q,停止YN30

在MATLAB軟件包中編程，共需3個Ｍ文件:,,

.主程序為.functionz=mylp(x)

%目標函數z=2*x(1)^2+x(2)^2-x(1)*x(2)-8*x(1)-3*x(2);

%轉化為求最小值問題31

function[sol,r1,r2]=randlp(a,b,n)

%隨機模擬解非線性規劃debug=1;a=0;

%試驗點下界b=10;

%試驗點上界n=1000;

%試驗點個數r1=unifrnd(a,b,n,1);

%n1階的［a,b］均勻分布隨機數矩陣r2=unifrnd(a,b,n,1);sol=[r1(1)r2(1)];z0=inf;fori=1:nx1=r1(i);x2=r2(i);lpc=lpconst([x1x2])；

iflpc==1z=mylp([x1x2])；

ifz<z0z0=z；

sol=[x1x2]；

end

endendToMATLAB(randlp)返回32實驗作業1．編一個福利彩票電腦選號的程序．334.某設備上安裝有4只型號規格完全相同的電子管，電子管壽命服從1000～2000h之間的均勻分布．電子管損壞時有兩種維修方案，一是每次更換損壞的那只；二是當其中1只損壞時4只同時更換．更換時間為換1只時需1h，4只同時換為2h．更換時機器因停止運轉每小時的損失為20元，又每只電子管價格10元，試用模擬方法確定哪一個方案經濟合理？5.

導彈追蹤問題：設位于坐標原點的甲艦向位于x軸上點A(1,0)處的乙艦發射導彈，導彈頭始終對準乙艦.如果乙艦以最大的速度(常數)沿平行于y軸的直線行駛，導彈的速度為5.模擬導彈運行的軌跡.乙艦行駛多遠時，導彈將它擊中？34初始化:i=0,k1=0,k2=0,k3=0i=i+1骰子點數?k1=k1+1k2=k2+1k3