實(shí)驗(yàn)人口遷移的動(dòng)態(tài)分析

1、實(shí)驗(yàn)八 人口遷移的動(dòng)態(tài)分析一實(shí)驗(yàn)?zāi)康尼槍?duì)人口遷移問題，建立以每單位時(shí)間為階段的常系數(shù)線性系統(tǒng)動(dòng)態(tài)變化模型使用Mathematica4.0 作矩陣運(yùn)算，并由模型討論該過程的極限狀態(tài)是否有穩(wěn)定解，用于分析、 預(yù)報(bào)、決策和控制該過程.通過討論狀態(tài)方程解的穩(wěn)定性，加深對(duì)矩陣特征值、特征向量的理解.實(shí)驗(yàn)原理把形如Uk 4 二 Auk的矩陣方程稱作 常差分方程組 或狀態(tài)方程，這里uk是列向量，A是矩陣.形式上，它是容 易解的，因?yàn)槊恳淮蔚加肁去乘.于是得到解為山=Aku°Uo是初始條件，A稱作一步狀態(tài)轉(zhuǎn)移矩陣.問題在于尋求某種快速計(jì)算幕Ak的方法，解決的關(guān)鍵是 A的特征值和特征向量.根據(jù)線性

2、代數(shù)的知識(shí)，n階方陣A與對(duì)角陣上相似（即A可對(duì)角化）的充要條件是 A有n個(gè)線性無關(guān)的特征向量.設(shè)A可以對(duì)角化，則存在可逆陣 S和對(duì)角陣上，使得A =這里S =（xX2，Xn）, -I二diagCi,d，n）, Xi是對(duì)應(yīng)A的特征值 打的特征向量（1 < i w n）.將上述結(jié)果用于uk = Au0，則自然有Uk =Aku0 NSIS/XSIS-1）（Stsju。= SVsUo （1）從而由矩陣乘法得出Uk =（X1, |,Xn）diag（雪,川,丸：）SUo丸壯 +| + Cn丸：Xn （2）由此，可以看出一般解是特解的一個(gè)線性組合，其中組合系數(shù)由初始條件決定：C0X1 Cn，：Xn =

3、U°,或 S（G, ,Cn）T =Uo,或（G, , Cn ）J S 仏從（2）中可以看出它與微分方程某些相近的地方，這將為我們下面要討論的狀態(tài)方程解的 穩(wěn)定性帶來方便.一般地，我們有如下結(jié)論：常系數(shù)線性系統(tǒng)uk q = Auk （ A可對(duì)角化），當(dāng)它的所有特征值|耐| <1時(shí)，它是完全穩(wěn)定的，即UkT o（ k TR），這保證了初始條件的微小變化所造成的影響會(huì)隨著 k的增加而趨于零；當(dāng)所有<1時(shí)，它是中性穩(wěn)定的，即uk有界； 當(dāng)至少有一個(gè)特征值 卩>1時(shí)，它是不穩(wěn)定的，即uk是無界的，也就是說穩(wěn)定性依賴于 A 的特征值這些由公式（2）很容易得到.三學(xué)習(xí)Mathem

4、atica命令1.方陣的幕 MatrixPower求方陣的幕An的命令的形式為MatrixPowerA , n 其中n為整數(shù)，當(dāng)n - _1時(shí)即求逆.例如：輸入aa = 1,0, 0, 1,1,0, 0, 1, 1;MatrixPoweraa, 5輸出為1,0, 0, 5, 1, 0, 10, 5, 1如果輸入MatrixPoweraa, -1則得到逆陣1,0, 0, -1, 1,0, 1, -1, 1還有一個(gè)求逆陣的命令,輸入aa1 = Inv erseaa同樣得到逆陣1,0, 0, -1, 1,0, 1, -1, 1.不過，如果求逆陣的幕，則用前一個(gè)命令較好.只須輸入MatrixPower

5、aa, -5得到輸出1,0, 0, -5, 1,0, 15, -5, 12. Do型循環(huán)結(jié)構(gòu)Do型循環(huán)結(jié)構(gòu)根據(jù)循環(huán)描述先計(jì)算循環(huán)次數(shù)，再作循環(huán)體，常用于有確定循環(huán)次數(shù)的循環(huán)結(jié)構(gòu).Do語句的一般形式為Do循環(huán)體，循環(huán)范圍.它有下列形式：Do表達(dá)式，k（計(jì)算表達(dá)式k次.）Do表達(dá)式，i ， imax（計(jì)算表達(dá)式imax次，其中i的值從1變到imax，每次 步長為1.）Do表達(dá)式，i ， imin ， imax （當(dāng) i 的值從 imin 變 到 imax、 步長為1，每次都計(jì)算表達(dá)式.）Do表達(dá)式，i ， imin ， imax， in creme nt（當(dāng) i 的值從 imin變至 U imax

6、、步長為in creme nt ，每次都計(jì)算表達(dá)式.）Do表達(dá)式，i , imin , imax , j , jmin , jmax，（當(dāng) i 的值從 imin 變到 imax、步長為1、當(dāng)j的值從jmin 變到j(luò)max、步長為1，每次都計(jì)算表達(dá)式.當(dāng)j完成一次 循環(huán)后，i的值增加1，以此類推.這就是所謂的Do循環(huán)嵌套.）Do表達(dá)式，i , imin , imax, in creme nt , j , jmin , jmax, jn creme nt ,（形 成一個(gè)Do循環(huán)嵌套，這時(shí)步長是指定值.）例如：輸入t = x; Dot = 1/（1 + k*t）, k, 2, 6, 2; t輸出為輸

7、入得到輸出1,12,12,241 2xDoPrinti, j, i, 2, j, i四實(shí)驗(yàn)內(nèi)容例1 對(duì)城鄉(xiāng)人口流動(dòng)作年度調(diào)查，發(fā)現(xiàn)有一個(gè)穩(wěn)定的向城鎮(zhèn)流動(dòng)的趨勢(shì)：每年，農(nóng)村居民的2.5%移居城鎮(zhèn)，而城鎮(zhèn)居民的1%遷出現(xiàn)在總?cè)丝诘?60 %位于城鎮(zhèn)假如城鄉(xiāng)總?cè)丝诒３植蛔?，并且人口流?dòng)的這種趨勢(shì)繼續(xù)下去，那么一年以后住在城鎮(zhèn)人口所占比例是多少？兩年以后呢？十年以后呢？最終呢？解 為了分析這個(gè)問題，開始時(shí)，令鄉(xiāng)村人口為y0，城鎮(zhèn)人口為z0, 一年以后有或?qū)懗删仃囆问絻赡暌院?，有y2<y2丿975125991000 y0 100zz1鄉(xiāng)村人口城鎮(zhèn)人口975100025<1000110099陸

8、丿100丿*975100025<10001 ,100 5 j99乜.丿100*9751 2/ 、/ 、即y21000100y。<Z2 J25990<1000100丿十年以后，有10 <z10 丿975100025110001、0100y990丿100丿事實(shí)上，它給出了一個(gè)差分方程組:9751、uk+ = Auk，這里，A =1000100，uk2599<1000100丿A的特征值和特征向量：輸入yk<Zk丿根據(jù)前面的討論，我們首先計(jì)算結(jié)果是:Eige nsystem975/1000,1/100,25/1000,99/100193/200, 1, -1, 1,

9、 2/5, 1從而可知,A可以對(duì)角化,并且對(duì)角陣為193、200I 1丿，相應(yīng)的25 再輸入1丿In verse-1,2/5,1,1得逆陣s57572757于是利用公式(1)，得到k年之后的分布:yky。-12 官 193 5200 丿=(y。 z。)r 572<193 V7+ ( z°-y°)55(200 丿5<7丿I 7丿7丿這就是我們所要的解容易看出：當(dāng)kis時(shí)，這個(gè)解會(huì)達(dá)到一個(gè)極限狀態(tài)y：Zdc7=(y° Z0)這里，總?cè)丝谌允莥0 - z0，與開始時(shí)一樣.但在此極限中人口的 鄉(xiāng)村易見：極限狀態(tài)與初始分布無關(guān)，且正是A的屬于特征值 上述例子有一

10、些很好的性質(zhì).因?yàn)榫仃嘇的每一列的和等于5/7在城鎮(zhèn)，而2/7 在 1的一個(gè)特征向量.1, 所以人口總數(shù)保持不 變又因?yàn)榫仃?A沒有負(fù)元素，而y0和z0也是非負(fù)的，所以鄉(xiāng)村和城鎮(zhèn)的人口數(shù)決不能為 負(fù).除此之外，每一個(gè)新的狀態(tài) uk1僅依賴于前一個(gè)uk，而與狀態(tài)u0,IH,ukj無關(guān).這個(gè)性質(zhì)稱為無后效性具有無后效性的隨機(jī)過程稱為馬爾可夫過程.由于A - I不可逆,5=1總是狀態(tài)轉(zhuǎn)移矩陣 A= (aij) nxn的一個(gè)特征值當(dāng)A可以 對(duì)角化時(shí)，A沒有某一個(gè)特征值的絕對(duì)值大于1.(可由Gerschgorin 圓盤定理得到.)如果其它的所有特征值的絕對(duì)值嚴(yán)格小于= 1,則公式(2)中的第一項(xiàng)完全處于

11、支配地位；其它的丄：將迅速地趨于0 ,從而ukcn = u：- ( kim).即，在一般情況下，' = 1 的一個(gè)特征向量是穩(wěn)定狀態(tài).本例中的狀態(tài)轉(zhuǎn)移矩陣A是二階方陣，所以我們可以利用 Mathematica4.0作出散點(diǎn)圖動(dòng)畫來觀察狀態(tài)方程的穩(wěn)定性u(píng) = 0.6, 0.4;A = 0.975, 0.01, 0.025, 0.99;pi_ := MatrixPowerA, i.u;fg=DoListPlotTablepi, i, 1, j, 5,PlotRange -> 0.29, 0.59, 0.41, 0.72, j, 1, 160, 5fp = Tablepi, i, 1,

12、200, 50.7A0.65A «0.6 *0.55-#:0.350.40.450.50.550.45H4-*圖 22.1這里fp是包含U1, U6, U11,，U200的數(shù)值表，從中可以看出uk收斂于(0.286 , 0.714 ),即(2/7，5/7). fg是以年距為5 一張張給出的逐點(diǎn)增加的散點(diǎn)圖.我們可以選中作出的圖形，即用鼠標(biāo)選擇圖像單元括號(hào)，點(diǎn)擊菜單Cell | An imate Selected Graphcis(若是 Mathematica2.2則單擊菜單 Graph | Ani mate Selected Graphcis)就可以觀看動(dòng)畫了 .從中可以看到散點(diǎn)最終

13、的聚集狀態(tài).圖1是動(dòng)畫序列中的最后一張散點(diǎn)圖.例2設(shè)狀態(tài)方程由矩陣 A= c 1定義因?yàn)?A是上三角陣，所以它的特征值是0 _<2丿主對(duì)角元素0和1/2.因?yàn)樗鼈兌夹∮?1,由公式(2)可知：從任何一個(gè)初始向量 u0出 發(fā)，狀態(tài)方程Uk1.=AUk的解趨于零也可以象例 1那樣作出散點(diǎn)圖動(dòng)畫，請(qǐng)讀者完成 .廣 1-10、例3 由矩陣A= -12-1定義的狀態(tài)方程.先求出它的特征值和特征向量，<0-11輸入Eige nsystem1,-1,0,-1,2,-1,0,-1,1結(jié)果是：0, 1,3, 1, 1, 1, -1,0, 1, 1, -2, 1.A有3個(gè)線性無關(guān)的特征向量，從而A可以

14、對(duì)角化.又因?yàn)樗幸粋€(gè)特征值3>1，由公式(2)可知：從任何一個(gè)非零初始向量U0出發(fā)，狀態(tài)方程的解都將是無界的.例4 討論由矩陣函數(shù)A(k)sin k(-1)kek01-cosk0e k2u = 0.6, 0.4;A = (-1)TExpSini/i, 0, 0, Exp(1 - Cosi)/(iA2);定義的Uk1=A(k)Uk的穩(wěn)定性，這里初值U0= ( 0.6 , 0.4 ).這個(gè)方程已經(jīng)超出了前面講畫出它的散點(diǎn)動(dòng)畫圖:過的常系數(shù)線性系統(tǒng)，但我們?nèi)岳肕athematica4.0pi := MatrixPowerA, i.u;DoListPlotTablepi, i, 1, j,P

15、lotRange -> -2, 2, 0.4, 0.47, j, 1,80僅從已繪出的uk的散點(diǎn)圖，我們可以發(fā)現(xiàn) uk是有界的，但是它卻是發(fā)散的 我們可以畫出 更多的點(diǎn)，如圖 22.2 :輸入ListPlotTableMatrixPowerA, i.0.6, 0.4, i, 1, 8000-1.5-1-0.50.511.5圖 22.2事實(shí)上，例4是一種混沌現(xiàn)象，利用它的偽隨機(jī)性可用于作 “一次性”的隨機(jī)密碼.若 每組密碼只用一次，那末就不容易被破譯出來五實(shí)驗(yàn)作業(yè)要求按如下步驟完成實(shí)驗(yàn)作業(yè).1）由實(shí)際問題寫出轉(zhuǎn)移矩陣，列出狀態(tài)方程2） 利用Mathematica4.0求出轉(zhuǎn)移矩陣 A的特征值、特征向量并判斷能否對(duì)角化？ 推導(dǎo)出極限狀態(tài)的表達(dá)式.注意要用精確值計(jì)算.3）如果轉(zhuǎn)移矩陣是二階的，做出迭代的散點(diǎn)圖，利用圖形觀察解的性態(tài)4） 討論最后結(jié)果的實(shí)際含義.1假定有一種昆蟲傳染病，在每個(gè)月內(nèi)，健康昆蟲的一半會(huì)染上病，而染病昆蟲的1/4會(huì)死亡，那么最終是否會(huì)有昆蟲健在？2. 假定某一物質(zhì)能以液態(tài)與氣態(tài)存在，又設(shè)在一段很短的時(shí)間內(nèi)，液體的1/10蒸發(fā),氣體的2/10凝結(jié)，那么該物質(zhì)是否存在一個(gè)平衡狀態(tài)，能保證有物質(zhì)的60%是氣態(tài)的？3. 假設(shè)有三個(gè)大型貨運(yùn)卡車中心.每個(gè)月中，在北京和在上海的卡車的一半開往廣州， 而其余一半留在原地廣州的卡車分成相等