实验人口迁移的动态分析

1、实验八 人口迁移的动态分析一实验目的针对人口迁移问题，建立以每单位时间为阶段的常系数线性系统动态变化模型使用Mathematica4.0 作矩阵运算，并由模型讨论该过程的极限状态是否有稳定解，用于分析、 预报、决策和控制该过程.通过讨论状态方程解的稳定性，加深对矩阵特征值、特征向量的理解.实验原理把形如Uk 4 二 Auk的矩阵方程称作 常差分方程组 或状态方程，这里uk是列向量，A是矩阵.形式上，它是容 易解的，因为每一次迭代都用A去乘.于是得到解为山=Aku°Uo是初始条件，A称作一步状态转移矩阵.问题在于寻求某种快速计算幕Ak的方法，解决的关键是 A的特征值和特征向量.根据线性

2、代数的知识，n阶方阵A与对角阵上相似（即A可对角化）的充要条件是 A有n个线性无关的特征向量.设A可以对角化，则存在可逆阵 S和对角阵上，使得A =这里S =（xX2，Xn）, -I二diagCi,d，n）, Xi是对应A的特征值 打的特征向量（1 < i w n）.将上述结果用于uk = Au0，则自然有Uk =Aku0 NSIS/XSIS-1）（Stsju。= SVsUo （1）从而由矩阵乘法得出Uk =（X1, |,Xn）diag（雪,川,丸：）SUo丸壮 +| + Cn丸：Xn （2）由此，可以看出一般解是特解的一个线性组合，其中组合系数由初始条件决定：C0X1 Cn，：Xn =

3、U°,或 S（G, ,Cn）T =Uo,或（G, , Cn ）J S 仏从（2）中可以看出它与微分方程某些相近的地方，这将为我们下面要讨论的状态方程解的 稳定性带来方便.一般地，我们有如下结论：常系数线性系统uk q = Auk （ A可对角化），当它的所有特征值|耐| <1时，它是完全稳定的，即UkT o（ k TR），这保证了初始条件的微小变化所造成的影响会随着 k的增加而趋于零；当所有<1时，它是中性稳定的，即uk有界； 当至少有一个特征值 卩>1时，它是不稳定的，即uk是无界的，也就是说稳定性依赖于 A 的特征值这些由公式（2）很容易得到.三学习Mathem

4、atica命令1.方阵的幕 MatrixPower求方阵的幕An的命令的形式为MatrixPowerA , n 其中n为整数，当n - _1时即求逆.例如：输入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还有一个求逆阵的命令,输入aa1 = Inv erseaa同样得到逆阵1,0, 0, -1, 1,0, 1, -1, 1.不过，如果求逆阵的幕，则用前一个命令较好.只须输入MatrixPower

5、aa, -5得到输出1,0, 0, -5, 1,0, 15, -5, 12. Do型循环结构Do型循环结构根据循环描述先计算循环次数，再作循环体，常用于有确定循环次数的循环结构.Do语句的一般形式为Do循环体，循环范围.它有下列形式：Do表达式，k（计算表达式k次.）Do表达式，i ， imax（计算表达式imax次，其中i的值从1变到imax，每次 步长为1.）Do表达式，i ， imin ， imax （当 i 的值从 imin 变 到 imax、 步长为1，每次都计算表达式.）Do表达式，i ， imin ， imax， in creme nt（当 i 的值从 imin变至 U imax

6、、步长为in creme nt ，每次都计算表达式.）Do表达式，i , imin , imax , j , jmin , jmax，（当 i 的值从 imin 变到 imax、步长为1、当j的值从jmin 变到jmax、步长为1，每次都计算表达式.当j完成一次 循环后，i的值增加1，以此类推.这就是所谓的Do循环嵌套.）Do表达式，i , imin , imax, in creme nt , j , jmin , jmax, jn creme nt ,（形 成一个Do循环嵌套，这时步长是指定值.）例如：输入t = x; Dot = 1/（1 + k*t）, k, 2, 6, 2; t输出为输

7、入得到输出1,12,12,241 2xDoPrinti, j, i, 2, j, i四实验内容例1 对城乡人口流动作年度调查，发现有一个稳定的向城镇流动的趋势：每年，农村居民的2.5%移居城镇，而城镇居民的1%迁出现在总人口的 60 %位于城镇假如城乡总人口保持不变，并且人口流动的这种趋势继续下去，那么一年以后住在城镇人口所占比例是多少？两年以后呢？十年以后呢？最终呢？解 为了分析这个问题，开始时，令乡村人口为y0，城镇人口为z0, 一年以后有或写成矩阵形式两年以后，有y2<y2丿975125991000 y0 100zz1乡村人口城镇人口975100025<1000110099陆

8、丿100丿*975100025<10001 ,100 5 j99乜.丿100*9751 2/ 、/ 、即y21000100y。<Z2 J25990<1000100丿十年以后，有10 <z10 丿975100025110001、0100y990丿100丿事实上，它给出了一个差分方程组:9751、uk+ = Auk，这里，A =1000100，uk2599<1000100丿A的特征值和特征向量：输入yk<Zk丿根据前面的讨论，我们首先计算结果是:Eige nsystem975/1000,1/100,25/1000,99/100193/200, 1, -1, 1,

9、 2/5, 1从而可知,A可以对角化,并且对角阵为193、200I 1丿，相应的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丿这就是我们所要的解容易看出：当kis时，这个解会达到一个极限状态y：Zdc7=(y° Z0)这里，总人口仍是y0 - z0，与开始时一样.但在此极限中人口的 乡村易见：极限状态与初始分布无关，且正是A的属于特征值 上述例子有一

10、些很好的性质.因为矩阵A的每一列的和等于5/7在城镇，而2/7 在 1的一个特征向量.1, 所以人口总数保持不 变又因为矩阵 A没有负元素，而y0和z0也是非负的，所以乡村和城镇的人口数决不能为 负.除此之外，每一个新的状态 uk1仅依赖于前一个uk，而与状态u0,IH,ukj无关.这个性质称为无后效性具有无后效性的随机过程称为马尔可夫过程.由于A - I不可逆,5=1总是状态转移矩阵 A= (aij) nxn的一个特征值当A可以 对角化时，A没有某一个特征值的绝对值大于1.(可由Gerschgorin 圆盘定理得到.)如果其它的所有特征值的绝对值严格小于= 1,则公式(2)中的第一项完全处于

11、支配地位；其它的丄：将迅速地趋于0 ,从而ukcn = u：- ( kim).即，在一般情况下，' = 1 的一个特征向量是稳定状态.本例中的状态转移矩阵A是二阶方阵，所以我们可以利用 Mathematica4.0作出散点图动画来观察状态方程的稳定性u = 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的数值表，从中可以看出uk收敛于(0.286 , 0.714 ),即(2/7，5/7). fg是以年距为5 一张张给出的逐点增加的散点图.我们可以选中作出的图形，即用鼠标选择图像单元括号，点击菜单Cell | An imate Selected Graphcis(若是 Mathematica2.2则单击菜单 Graph | Ani mate Selected Graphcis)就可以观看动画了 .从中可以看到散点最终

13、的聚集状态.图1是动画序列中的最后一张散点图.例2设状态方程由矩阵 A= c 1定义因为 A是上三角阵，所以它的特征值是0 _<2丿主对角元素0和1/2.因为它们都小于 1,由公式(2)可知：从任何一个初始向量 u0出 发，状态方程Uk1.=AUk的解趋于零也可以象例 1那样作出散点图动画，请读者完成 .广 1-10、例3 由矩阵A= -12-1定义的状态方程.先求出它的特征值和特征向量，<0-11输入Eige nsystem1,-1,0,-1,2,-1,0,-1,1结果是：0, 1,3, 1, 1, 1, -1,0, 1, 1, -2, 1.A有3个线性无关的特征向量，从而A可以

14、对角化.又因为它有一个特征值3>1，由公式(2)可知：从任何一个非零初始向量U0出发，状态方程的解都将是无界的.例4 讨论由矩阵函数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的稳定性，这里初值U0= ( 0.6 , 0.4 ).这个方程已经超出了前面讲画出它的散点动画图:过的常系数线性系统，但我们仍利用Mathematica4.0pi := MatrixPowerA, i.u;DoListPlotTablepi, i, 1, j,P

15、lotRange -> -2, 2, 0.4, 0.47, j, 1,80仅从已绘出的uk的散点图，我们可以发现 uk是有界的，但是它却是发散的 我们可以画出 更多的点，如图 22.2 :输入ListPlotTableMatrixPowerA, i.0.6, 0.4, i, 1, 8000-1.5-1-0.50.511.5图 22.2事实上，例4是一种混沌现象，利用它的伪随机性可用于作 “一次性”的随机密码.若 每组密码只用一次，那末就不容易被破译出来五实验作业要求按如下步骤完成实验作业.1）由实际问题写出转移矩阵，列出状态方程2） 利用Mathematica4.0求出转移矩阵 A的特征值、特征向量并判断能否对角化？ 推导出极限状态的表达式.注意要用精确值计算.3）如果转移矩阵是二阶的，做出迭代的散点图，利用图形观察解的性态4） 讨论最后结果的实际含义.1假定有一种昆虫传染病，在每个月内，健康昆虫的一半会染上病，而染病昆虫的1/4会死亡，那么最终是否会有昆虫健在？2. 假定某一物质能以液态与气态存在，又设在一段很短的时间内，液体的1/10蒸发,气体的2/10凝结，那么该物质是否存在一个平衡状态，能保证有物质的60%是气态的？3. 假设有三个大型货运卡车中心.每个月中，在北京和在上海的卡车的一半开往广州， 而其余一半留在原地广州的卡车分成相等