第五章 差分方程模型.doc

134第五章 差分方程模型在第四章中，我们利用微分方程方法研究了一些连续变化的变量。如果将变量离散化，即可得到相应的差分方程模型，为了方便不熟悉差分方程的读者，先对本章用到的差分方程的知识作一简略介绍。5.1差分方程简介一、差分方程及其通解以表示时间，规定只取非负整数。表示第一周期初，表示第二周期初等。记为变量在时刻时的取值，则称为的一阶差分，称为yt的二阶差分。类似地，可以定义yt的阶差分。由、及的差分给出的方程称为差分方程，其中含的最高阶差分的阶数称为该差分方程的阶。差分方程也可以写成不显含差分的形式。例如，二阶差分方程也可改写成。满足差分方程的序列称为此差分方程的解。类似于微分方程情况，若解中含有的独立常数的个数等于差分方程的阶数时，则称此解为该差分方程的通解。若解中不含任意常数，则称此解为满足某些初值条件的特解，例如，考察两阶差分方程： 易见与均是它的特解，而则为它的通解，其中，为两个任意常数。类似于微分方程，称差分方程 （1）为阶线性差分方程，当时称其为阶非齐次线性差分方程，而 （2）称为方程（1）对应的齐次线性差分方程。若（1）中所有的均为与无关的常数，则称其为常系数差分方程，即阶常系数线性差分方程可写成 （3）其对应的齐次方程为 （4）容易证明，若序列与均为方程（4）的解，则也是方程（4）的解，其中，为任意常数，这说明，齐次方程的解构成一个线性空间（解空间）。若是方程（4）的解，是方程（3）的解，则也是方程（3）的解。方程（3）可用如下的代数方法求其通解：（步一）先求解对应的特征方程 （5）（步二）根据特征根的不同情况，求齐次方程（4）的通解。 情况1 若特征方程（5）有个互不相同的实根，则齐次方程（4）的通解为 (，为任意常数)情况2 若是特征方程（5）的重根，通解中对应于的项为， （）为任意常数。 情况3 若特征方程（5）有单重复根，通解中对应于的项为，其中，为任意常数，为的模，为的幅角。情况4 若为特征方程（5）的k重复根，则通解对应于的项为， （）为任意常数。（步三） 求非齐次方程（3）的一个特解。若为方程（4）的通解，则非齐次方程（3）的通解为。求非齐次方程（3）的特解一般要用到常数变易法，计算较繁。对特殊形式的b(t)也可使用待定系数法。例如，当，为的次多项式时可以证明：若不是特征根，则非齐次方程（3）有形如的特解， 也是的次多项式；若是重特征根，则方程（3）有形如的特解。进而可利用待定系数法求出，从而得到方程（3）的一个特解。例 求解两阶差分方程。解 对应齐次方程的特征方程为，其特征根为，对应齐次方程的通解为： 原方程有形如的特解。代入原方程求得，故原方程的通解为：二、差分方程的平衡点及其稳定性在应用差分方程研究问题时，一般不需要求出方程的通解，在给定初值后，通常可用计算机迭代求解，但我们常常需要讨论解的稳定性。1、 一阶线性常系数差分方程 （6）在（6）式中令得到的代数方程的根称为差分方程（6）的平衡点。如果时，则称平衡点是稳定的，否则是不稳定的。由（6）式得: （7）即方程（6）的解可表为 （8）其中由初始值确定。显然，由（8）式可知差分方程（6）的平衡点稳定的充要条件是 | a | 1 （9）顺便指出，对于维向量x(k)和常数矩阵A构成的方程组 x( k+1) + Ax( k ) = 0 （10）其平衡点稳定的条件是A的特征根（）均有 （11）即均在复平面上的单位圆内。这个结果可由将化为对角阵（或Jordan）阵得到。2、二阶线性常系数差分方程， （12）考察方程（12）的平衡点（0）的稳定性，设（12）的特征方程的根为，则不难验证，（12）的通解可表为 （13）其中常数，由初始条件，确定。由（13）知，当且仅当 ， （14）时方程（12）的平衡点才是稳定的。3、 非齐次线性差分方程 （15）方程（15）的平衡点的稳定性和方程（12）相同。 二阶方程的上述结果可以推广到阶线性方程，即稳定平衡点的条件是特征根次代数方程的根（）均有。 4、一阶非线性差分方程 (16)方程（16）的平衡点由代数方程解出。为分析的稳定性，将方程（16）的右端在点作Taylor展开，只取一次项，（16）近似为 (17)（17）是（16）的近似线性方程，也是（17）的平衡点。关于线性方程（17）的稳定平衡点的讨论已由（6）（9）给出，而当时方程（16）与（17）平衡点的稳定性相同。于是得到当 （18）时，非线性方程（16）的平衡点是稳定的；当 （19）时，非线性方程（16）的平衡点是不稳定的。5.2市场经济中的蛛网模型问题 在自由贸易市场上常会出现这样的现象：一个时期以来当某种消费品如猪肉的上市量远大于需求时，由于销售不畅致使价格下跌，生产者发现养猪赔钱，于是转而经营其它农副业。过一段时间猪肉上市量就会大减，供不应求将导致价格上涨。生产者看到有利可图，又重操旧业，这样下一个时期会重现供大于求、价格下降的局面。商品数量和价格的这种振荡现象在自由竞争的市场经济中常常是不可避免的。进一步观察可以发现，振荡有两种完全不同的形式，一种是振幅逐渐减小，市场经济趋向平稳，另一种是振幅越来越大，如果没有外界如政府的干预，将导致经济崩溃。试建立数学模型描述这种现象，研究经济趋向平稳的条件，并讨论当经济趋向不稳定时政府可能采取的干预措施。蛛网模型 商品在市场上的数量和价格出现反复的振荡，是由消费者的需求关系和生产者的供应关系决定的。记商品第时段的上市数量为，价格为。这里我们把时间离散化为时段，1个时段相当于商品的1个生产周期，如蔬菜、水果是一个种植周期，肉类是牲畜的饲养周期。同一时段商品的价格取决于数量，设 （1）它反映消费者对这种商品的需求关系，称需求函数。因为商品的数量越多价格越低，所以在图1中用一条下降曲线表示它，称需求函数。 下一时段商品的数量由上一时段价格决定，设 ，或 （2）这里是的反函数。或反映生产者的供应关系，称供应函数。因为价格越高生产量（即下一时段的商品数量）越大，所以在图中供应曲线是一条上升曲线。 图中两条曲线相交于点。是平衡点，其意义是，一旦在某一时段有，则由（1），（2）可知，即以后各时段商品的数量和价格将永远保持在点。但是实际生活中的种种干扰使得数量和价格不可能停止在点，不防设偏离（如图1）。我们分析随着的增加，的变化。yfgyfg0xP0P1P2P3P4图1 需求曲线 f 和供应曲线 g , P0是稳定平衡点图2 P0 是不稳定平衡点xx1x3x0x20y0y1y2P1P0P2P3P4商品数量给定后，价格由曲线上的点决定，下一时段的数量由曲线上的点决定，又由曲线上的点决定，这样得到一系列的点，在图1上这些点将按箭头所示方向趋向，表明是稳定的平衡点，意味着市场经济（商品的数量和价格）将趋向稳定。但是如果需求函数和供应函数由图2的曲线所示，则类似的分析发现，市场经济将按，的规律变化而远离，即不是稳定的平衡点，意味商品数量和价格将出现越来越大的振荡。图1和图2中折线形似蛛网，所以这种用需求曲线和供应曲线分析市场经济稳定性的图示法在经济学中称蛛网模型。实际上，需求曲线和供应曲线的具体形式通常是根据各个时段商品的数量和价格的一系列统计资料得到的。一般地说，取决于消费者对这种商品的需要程度和他们的消费水平，则与生产者的生产能力、经营水平等因素有关。比如当消费者收入增加时，会向上移动；当生产能力提高时，将向右移动。一旦需求曲线和供应曲线被确定下来，商品数量和价格是否趋向稳定，就完全由这两条曲线在平衡点附近的形状决定。只要分析一下图1和图2的不同之处就会发现，在附近，图1的比平缓，而图2的比陡峭。记在点斜率的绝对值（因为它是下降的）为，在点的斜率为，图形的直观告诉我们，当 （3）时是稳定的（图1），当 （4）时是不稳定的（图2）.由此可见，需求曲线越平，供应曲线越陡，越有利于经济稳定。为了进一步分析这种现象，下面给出蛛网模型的另一种表达形式差分方程。 差分方程模型 在点附近可以用直线来近似曲线和，设（1），（2）式分别近似为 ， （5） ， ， （6）从二式中消去可得 （7）（7）是一阶线性常系数差分方程，对递推不难得到 ， （8）容易看出，当，时，即点稳定的条件是 或 （9）当，时，即点不稳定的条件是 或 （10）注意到（5），（6）式中的定义，有，,所以条件（9），（10）与蛛网模型中的直观结果（3），（4）式是一致的。 模型解释 首先考察参数，的含义。由（5）式可知，表示商品供应量减少1个单位时价格的上涨幅度；由（6）式可知，表示价格上涨1个单位时（下一时期）商品供应的增加量。所以的数值反映消费者对商品需求的敏感程度，如果这种商品是生活必需品，消费者处于持币待购状况，商品数量稍缺，人们立即蜂拥抢购，那么会比较大；反之，若这种商品非必需品，消费者购物心理稳定，或者消费水平低下，则较小。的数值反映生产经营者对商品价格的敏感程度，如果他们目光短浅，热衷于追逐一时的高利润，价格稍有上涨就大量增加生产，那么会比较大；反之，若他们素质较高，有长远的计划，则较小。 根据的意义容易对市场经济稳定与否的条件（9），（10）作出解释。当供应函数即固定时，越小，需求曲线越平，表明消费者对商品需求的敏感程度越小，（9）式越容易成立，有利于经济稳定。当需求函数即固定时，越小，供应曲线越陡，表明生产者对价格的敏感程度越小，（9）式也容易成立，有利于经济稳定。反之，当较大，表明消费者对商品的需求和生产者对商品的价格都很敏感，则会导致（10）式成立，经济不稳定。应该指出，和都是有量纲的，它们的大小都应在同一量纲单位下比较。同时，和的量纲互为倒数，所以是无量纲量，就可以与1比较大小了。经济不稳定时的干预办法 基于上述分析可以看出，当市场经济趋向不稳定时政府有两种干预办法。一种办法是使尽量小，不妨考察极端情况，即需求曲线变为水平（图3），这时不论供应曲线如何，即不管多大，（9）式总成立，经济总是稳定的。实际上这种办法相当于政府控制物价，无论商品数量多少，命令价格不得改变。另一种办法是使尽量小，极端情况是，即供应曲线变为竖直（图4），于是不论需求曲线如何，即不管多大，（9）式总成立，经济也总是稳定的。实际上这相当于控制市场上的商品数量，当供应量少于需求时，从外地收购或调拨，投入市场；当供过于求时，收购过剩部分，维持商品上市量不变。显然，这种办法需要政府具有相当强大经济实力。x0yfgy0yfg0xx0图3 第一种干预办法示意图图4 第二种干预办法示意图模型的推广 如果生产者的管理水平和素质更高一些，他们在决定商品生产数量时，不是仅根据前一时期的价格，而是根据前两个时期的价格和。为简单起见不妨设根据二者的平均值（）/2，于是供应函数（2）式表为 （11）相应地，（2）式的线性近似表达式（6）修改为 （12）其中是平均价格上涨1个单位时的增量。又设需求函数仍由（1），（5）式表示。则由（5），（12）式得 (13)(13)是二阶线性常系数差分方程。为寻求时，即点稳定的条件，不必解方程（13），只须利用判断稳定的条件方程特征根均在单位圆内由方程（13）的特征方程容易算出其特征根为 （14）当时显然有从而，在单位圆外。下面设，由（14）式可以算出 （15）要使特征根均在单位圆内，即，必须 （16）这就是点稳定的条件。与原有模型中点稳定的条件（9）式相比，参数，的范围放大了。可以想到，这是因为生产者的管理水品和素质提高，对市场经济的稳定起着有利影响的必然结果。5.3汽车租赁公司的运营图1 3个城市间汽车数量比例的转移ABC0.60.10.10.30.10.70.30.20.6当我们研究的对象是若干变量构成的一个微量的离散动态过程时，用差分方程组描述比较方便。问题 一家汽车租赁公司在3个相邻的城市运营，为方便顾客起见公司承诺，在一个城市租赁的汽车可以在任意一个城市归还。根据经验估计和市场调查，一个租赁期内在A市租赁的汽车在A、B、C市归还的比例分别为0.6，0.3，0.1；在B市租赁的汽车在A、B、C市归还的比例分别为0.2，0.7，0.1；在C市租赁的汽车在A、B、C市归还的比例分别为0.1，0.3，0.6。若公司开业时将600辆汽车平均分配到3个城市，建立运营过程中汽车数量在3个城市间转移的模型，并讨论时间充分长以后的变化趋势。模型及其求解 记第k个租赁期末公司在A，B，C市的汽车数量分别为x1(k), x2(k), x3(k)，3个城市间汽车数量比例的转移可直观地用图1给出，按照这个转移关系容易写出第k+1个租赁期末公司在A，B，C市的汽车数量为（k0，1，2， ） （1）记向量x(k)= x1(k), x2(k), x3(k)T ，矩阵 A （2）则（1）式可表为 ， （3）给定初始值，可以用(3)式计算各个租赁期3个城市汽车数量的变化。如开始时600辆汽车平均分配到3个城市，用MATLAB计算的程序如下：A=0.6,0.2,0.1;03,0.7,0.3;0.1,0.1,0.6;x(:,1)=200,200,200; % 赋初值n=10;for k=1:n x(:,k+1)=A*x(:,k); % 按照（3）式迭代计算endround(x),k=0:10;plot(k,x),grid,gtext(x1(k),gtext(x2(k),gtext(x3(k),得到的结果见表1和图2（a）。表1 （3）式的计算结果（开始时600辆汽车平均分配到3个城市）k012345678910x1(k)200180176176178179179180180180180x2(k)200260284294297299300300300300300x3(k)200160140130125123121121120120120可以看到，时间充分长以后3个城市汽车数量趋于180，300，120.这个结果与初始条件有关吗？如果开始时600辆汽车全分配给A市，可得结果如表2和图2(b). 由上面的计算可以猜想，时间充分长以后3个城市的汽车数量趋向稳定，并且稳定值与600辆汽车的初始分配无关。表2 （3）式的计算结果（开始时600辆汽车全分配到A市）k012345678910x1(k)600360258214195187183181181180180x2(k)0180252281292297299300300300300x3(k)06090105113116118119120120120x1(k)x2(k)x3(k)x2(k)x1(k)x3(k)x1(k)x2(k)x3(k) (a) 表1的图示 (b) 表2的图示图2 （3）式的计算结果结果分析 为了证实上面的猜想，记稳定值为x ，由（3）式x应满足A x = x （4）这表明矩阵A的一个特征根1 ，且x是对应的特征向量。事实上，从矩阵各列之和均为1 立即可知，A有特征根1（为什么？），且特征向量x可由（4）式直接算出（当然应该利用x的各分量之和等于初始值600）.5.4 按年龄分组的种群增长问题 野生或饲养的动物因繁殖而增加，因自然死亡和人为屠杀而减少，不同年龄动物的繁殖率、死亡率有较大差别，因此在研究某一种群数量的变化时，需要考虑按年龄分组的种群增长。将种群按年龄等间隔地分成若干个年龄组，时间也离散化为时段，给定各年龄组种群的繁殖率和死亡率（在稳定环境下不妨假定它们与时段无关），建立按年龄分组的种群增长模型，预测未来各年龄组的种群数量，并讨论时间充分长以后的变化趋势。模型及其求解 将种群按年龄大小等间隔地分成个年龄组，比如每10岁或5岁为1个年龄组。与年龄的离散化相对应，时间也离散化为时段，并且时段的间隔与年龄区间大小相等，即以10年或5年为1个时段。种群是通过雌性个体的繁殖而增长的，所以用雌性个体数量的变化为研究对象比较方便，下面提到的种群数量均指其中的雌性。记时段第年龄组的种群数量为，；，第年龄组的繁殖率为，即第年龄组每个（雌性）个体在1个时段内平均繁殖的数量；第年龄组的死亡率为，即第年龄组1个时段内死亡数与总数之比，称为存活率。和可由统计资料获得，而的变化规律由以下的基本事实得到：时段第1年龄组种群数量是时段各年龄组繁殖数量之和，即 （1）而时段第年龄组的种群数量是时段第年龄组在存活下来的数量，即 （2）记时段种群按年龄组的分布向量为 （3）将和排成如下的矩阵 （4）则（1），（2）可表为 （5）当矩阵和按年龄组的初始分布向量已知时，可以预测任意时段种群各年龄组的分布为 ， （6）为方便分析，将归一化后的向量记做。例 设一种群分成5个年龄组，繁殖率，存活率，,各年龄组现有数量均为100只，由（4）式，（5）式用MATLAB计算，的程序如下：b=0,0.2,1.8,0.8,0.2;s=diag0.5,0.8,0.8,0.1; %对角阵，对角元素为0.5,0.8,0.8,0.1L=b;s,zeros(4,1); %按照（4）式构造矩阵Lx(:,1)=100*ones(5,1); %赋初值n=30;for k=1:n x(:,k+1)=L*x(:,k); %按照（5）式迭代计算endround(x),y=diag(1./sum(x); %为向量x归一化做的计算z=x*y, % z是向量x的归一化k=0:30;subplot(1,2,1),plot(k,x),grid, %在一个图形窗内画两张图subplot(1,2,2),plot(k,x),grid得到的结果见表1，表1和图1表1 的计算结果012342627282930100300220155265393403412423434100501501107719019620120621110080401208814915215716116510080643296117120122126129100108631112121213表2 的计算结果0123272829300.20000.57690.45640.36580.45640.45530.45590.45620.20000.09620.31120.25990.22240.22290.22180.22220.20000.15380.08300.28360.17260.17370.17370.17300.20000.15380.13280.07560.13550.13490.13540.13580.20000.01920.01660.01510.01320.01320.01310.0132 （a）x(k) 的图形（自上而下为x1(k)至x5(k)）（b）的图形（自上而下为至）图1 （5）式的计算结果结果分析 从表2和图1（b）可以看出，时间充分长以后种群按年龄组的分布向量趋向稳定，这种状况完全由矩阵L的性质决定。 根据和的定义，矩阵L中的元素满足 ，（若某个，那么第年龄组应取消） （7），且至少一个 （8）满足（7），（8）的矩阵L称Leslie矩阵（简称L矩阵）。关于L矩阵我们不加证明地叙述两个定理。定理1 L矩阵有唯一的正特征根，且它是单重的，对应特征向量为 （9）L矩阵的其他个特征根都满足(=2，3，) （10） 定理2 若L矩阵第一行有两项顺次的元素，都大于零，则（10）式中仅不等号成立，即0)和N(0)为何值，当时都有.那么对于方程（1）的差分形式（2），是否也有同样的性质，即时呢？下面将会看到，回答这个问题并不简单，而且将引出一个十分有趣的现象。平衡点及稳定性 代替（2）式我们讨论方程（6）的平衡点及其稳定性（非线性差分方程平衡点稳定性判据见5.1节。因为，由（4）知，为了求（6）的平衡点，解代数方程 （7）得（6）的非零平衡点为 （8）利用（4）、（5）式可以验证，x*相当于原方程（2）的非零平衡点y*N。为分析x*的稳定性，计算 （9）根据x*稳定的条件| f (x*) | 3图2 方程（6）的图解法（xx*)在条件（10）下xk收敛于x*的状况可以通过方程（6）的图形解法清楚地表示出来。以x为横坐标作和的图形（图1），曲线和直线交点的横坐标为平衡点x*。对于初值x0由方程（6）求x1 ，x2，的过程表示为图上带箭头的折线。当时的过程基本上是单调的（见图1（1）；而当时的过程则会出现形如蛛网模型图1那样的衰减振荡（图1（2）。图1 方程（6）的图解法（xx*)0xyy=xx0x1x*1x2y=f(x)1/2b/4(1) 1b20xyy=xx0x1x*1x2y=f (x)1/2(2) 2b3当时，虽然方程（6）仍可形式地求解，但不稳定，其图解法如图2所示，出现形如蛛网模型图2那样的发散振荡（）。数值计算 b由小到大取不同的数值，用方程（6）作计算（初值均取x00.2），结果xk（k1，2，100）见表1。表1 不同值下方程（6）计算结果xkkb=1.7b=2.6b=3.3b=3.45b=3.55b=3.5700.2000 0.2000 0.2000 0.2000 0.2000 0.2000 10.2720 0.4160 0.5280 0.5520 0.5680 0.5712 20.3366 0.6317 0.8224 0.8532 0.8711 0.8744 30.3796 0.6049 0.4820 0.4322 0.3987 0.3921 40.4004 0.6214 0.8239 0.8466 0.8510 0.8509 50.4081 0.6117 0.4787 0.4480 0.4500 0.4529 60.4107 0.6176 0.8235 0.8532 0.8786 0.8846 70.4114 0.6141 0.4796 0.4322 0.3785 0.3645 80.4117 0.6162 0.8236 0.8466 0.8351 0.8270 90.4117 0.6149 0.4794 0.4479 0.4888 0.5109 100.4118 0.6157 0.8236 0.8531 0.8871 0.8921 110.4118 0.6152 0.4794 0.4322 0.3557 0.3437 120.4118 0.6155 0.8236 0.8467 0.8135 0.8053 130.4118 0.6153 0.4794 0.4479 0.5385 0.5598 140.4118 0.6154 0.8236 0.8531 0.8822 0.8797 150.4118 0.6154 0.4794 0.4323 0.3688 0.3777 160.4118 0.6154 0.8236 0.8467 0.8264 0.8391 170.4118 0.6154 0.4794 0.4479 0.5093 0.4820 180.4118 0.6154 0.8236 0.8531 0.8872 0.8913 190.4118 0.6154 0.4794 0.4323 0.3553 0.3457 200.4118 0.6154 0.8236 0.8467 0.8132 0.8076 210.4118 0.6154 0.4794 0.4478 0.5394 0.5548 220.4118 0.6154 0.8236 0.8531 0.8820 0.8818 230.4118 0.6154 0.4794 0.4323 0.3695 0.3722 240.4118 0.6154 0.8236 0.8467 0.8270 0.8342 250.4118 0.6154 0.4794 0.4478 0.5079 0.4939 260.4118 0.6154 0.8236 0.8531 0.8873 0.8924 270.4118 0.6154 0.4794 0.4323 0.3550 0.3429 280.4118 0.6154 0.8236 0.8467 0.8129 0.8044 290.4118 0.6154 0.4794 0.4478 0.5399 0.5617 300.4118 0.6154 0.8236 0.8531 0.8818 0.8789 310.4118 0.6154 0.4794 0.4324 0.3699 0.3800 320.4118 0.6154 0.8236 0.8467 0.8274 0.8411 330.4118 0.6154 0.4794 0.4478 0.5070 0.4772 340.4118 0.6154 0.8236 0.8531 0.8873 0.8906 350.4118 0.6154 0.4794 0.4324 0.3549 0.3477 360.4118 0.6154 0.8236 0.8467 0.8128 0.8097 370.4118 0.6154 0.4794 0.4477 0.5402 0.5501 380.4118 0.6154 0.8236 0.8531 0.8818 0.8836 390.4118 0.6154 0.4794 0.4324 0.3701 0.3673 400.4118 0.6154 0.8236 0.8467 0.8276 0.8296 410.4118 0.6154 0.4794 0.4477 0.5065 0.5046 420.4118 0.6154 0.8236 0.8531 0.8874 0.8924 430.4118 0.6154 0.4794 0.4324 0.3549 0.3427 440.4118 0.6154 0.8236 0.8467 0.8127 0.8042 450.4118 0.6154 0.4794 0.4477 0.5403 0.5621 460.4118 0.6154 0.8236 0.8531 0.8817 0.8787 470.4118 0.6154 0.4794 0.4325 0.3702 0.3805 480.4118 0.6154 0.8236 0.8468 0.8277 0.8415 490.4118 0.6154 0.4794 0.4477 0.5062 0.4762 500.4118 0.6154 0.8236 0.8530 0.8874 0.8905 510.4118 0.6154 0.4794 0.4325 0.3548 0.3482 520.4118 0.6154 0.8236 0.8468 0.8127 0.8102 530.4118 0.6154 0.4794 0.4476 0.5404 0.5489 540.4118 0.6154 0.8236 0.8530 0.8817 0.8840 550.4118 0.6154 0.4794 0.4325 0.3