4微分方程的解及解的稳定性

1、第四讲微分方程解的稳定性k(t)k(t)上一讲，我们利用最大值原理讨论了新古典经济增长模型，得到了两个方 程，一个是状态变量的转移方程，另一个是欧拉方程。这两个方程构成了包含 状态变量和控制变量的二元一次方程组。c(t)-ok(t)3 =1.k(t)；,一 ：一、.c(t)这个方程组是一个非线性微分方程组，一般情况下，非线性方程组不存在 解析解，即方程组的解不能用初等函数来表示。因此，他们的性质需要借助其 他方法来了解。微分方程：变量为导数的方程叫做微分方程。常微分方程：只有一个自变量的微分方程叫做常微分方程。偏微分方程：有两个或两个以上自变量的方程叫做偏微分方程。微分方程的阶：微分方程中变量

2、的导数最高阶叫做方程的阶。线性方程：方程的形式是线性的。例如,方程ai y(t) + a2 y(t)十ag y(t)+x(t) = 0是一个二阶线性常微分方程。又如，索洛-斯旺模型的基本方程是一个非线性方程:k(t) =sk(t) -、.k(t)再如，拉姆齐模型的动态是下列微分方程组的解：处=k(t)L.皿一 k(t)k(t)也=二 k(t) P-、.c(t)一阶微分方程一阶微分方程可以用下面的方程表示半=f(x,y)(1.1)dx其中，函数f:RMRT R是连续可微函数。最简单的微分方程是(1.2)(1.3)f(x) dx它的解可表示为不定积分:y = f (x)dx c其中，F(x)=f(

3、x)dx表示任意一个被被积函数，c为任意常数。当然,我们也 . 、- 一一x可以确定任意一个被积函数，例如，令 F(x)= Jf(x)dx=L f(t)dt,则(2.2)的不定 积分可表示为 xy= j f (t)dt +c这时，不定积分仍然代表无穷多条曲线，如果给出初始条件 y(0) = y。，则，上面 微分方程的解就是 xy= 0 f (t)dt +y0(1.4)二、常见的一阶微分方程解法一阶线性微分方程一阶线性微分方程的一般形式为dy p(x)y = g(x)(2.1)dx边界条件(即初始条件)y(0) = y。 x为求解线性微分万程，在万程的两边同乘以exp j0 p(t)dt ,则方

4、程的左边为dy-exp 0 p(t)dtxp(x)exp。p(t)dt yd y exp 0 p(t)dtdx所以x TOC o 1-5 h z d y exp 0 p(t)dtx =g(x)exp p(t)dt(2.2) dx0方程(2.2)的解为xxxy =exp - 0 p(t)dt 刊.0 g(x)exp ,0 p(t)dt c (2.3).可分离变量的微分方程一个方程是可分离变量的，如果它可以写成下列形式f (x)dx = g(y)dy这类方程的解只需在方程两边同时积分即可。(2.4)f (x)dx = g(y)dy.可化为可分离变量或线性方程的贝努利方程方程dy p(x)y =yn

5、g(x)(2. dz ,、,、P(x)z =g(x) n -1 dx 这样，贝努利方程就转化为线性方程。 .恰当方程考虑非线性方程M(x,y) N(x,y)dy =0(2.7)dx或者M (x, y)dx N(x, y)dy = 0 如果存在函数(x, y) 满足M (x, y)dx N(x, y)dy = d (x, y)则称方程(2.7)是恰当方程，其解6(x,y)=c。dx叫做贝努利方程。其中，n为正整数。方程(2.5)两边同除以yn/ p(x)y1g(x)1 -n1 dy /、j ,、 一P(x)y = g(x) n1 dx1 .n(2.6)z 二 y例1,方程xdx + ydy =0

6、的解是xy=c.方程 dy +(ln y)dx =0 的解是 xln y = c . y三、一阶常微分方程的图解法对于线性常微分方程而言，目前已经有完整的理论，方程的解也可以用明确 的解析表达式来表示。但是，对于非线性方程而言，除了个别特殊的形式之外， 一般是没有办法获得解析表达式的，甚至根本不存在解析表达式。我们希望在没 有明确的解析表达式的情况下，仍然了解方程的解的性质。例2,索罗-斯旺模型的基本方程k(t) =sk(t)二 一、.k(t)(3.1)k表示资本存量，6表示资本折旧率，a表示资本的收入份额。该方程表示资本 存量的净增加等于总储蓄与总折旧之间的差额。先求稳定点。令k(t) =0

7、,得s(k(t)r5 -k(t) = 0可以求得两个解，k(t) = 0,(k(t) ) =9/s严力由于k -0,再判断稳定点稳定性。* 0,k kk(t) =s(k(t)尸6 k(t)=0,k=k*k，J根据 处为=口 S(k(t)产-6= 0 ,可得(k(t)* =(每/sa产)，k(t)在 dk(k(t) * =g/ss -3)有最大值，在(k(t)=付6支)1(d)的左边大于0,是k的增函数；在(k(t)* =(6/sa3的右边小于0,是k的减函数。即：L -*. 0, k k=a s(k(t)产-6 = * = 0, k = k*dk*k四、一元高阶线性微分方程与多元微分方程组以二

8、阶线性微分方程为例：ay(t) a2 y(t) a3 y(t) x(t) =0令z(t) = y(t),则，z(t) = y(t),于是该二阶线性微分方程就可以用一个一阶线性微分方程组来表示：az(t) a?z a3 y(t) x(t) =0y(t) = z(t)或者 TOC o 1-5 h z a2a31z(t)= - z(t)-y(t) x(t) aiaiaiJ(t) =z(t)由此看来，一个二阶微分方程就可以用一个一阶线性微分方程组来表示。同样道理，任何一个更高阶的微分方程，可以化成一个一阶微分方程组。因此， 要了解高阶微分方程的性质，只要研究一阶微分方程组的性质即可。1.最简单的线性线

9、性方程组：对角矩阵系统。yi(t) =a1iyi(t)y2(t) =a22 y2(t)写成矩阵形式就是：5y；(t) A 0 、,yi(t)、M(t)八 0a22 Xy2(t),系数矩阵有两个特征根，分别是aii和a?。方程的解yi(t) =ea11t - ciy2(t) =ea22t c2情形 1, aii 0,且 a22 A0： yi(t0) 0,y2(t0) 0, yi, y2 都随着时间的推移而 增加。状态不稳定。情形 2, a11c0,且 a22 0, y2(t0) 0, yi, y2都随着时间的推移而 下降。状态稳定。情形3,a11A 0,且a220,y2(t0)0,yi都随着时间

10、的推移而增加y2随着时间的推移而下降。状态为鞍点稳定。2. 一般非对角线性系统:y；(t)、由2 (t) /aii22ai2丫 yi甲 xi(t)a22 人y2(t)JlX2 该方程组的矩阵形式为Y(t) -AY(t) X(t)根据矩阵理论，对于矩阵 A,存在矩阵V,使V-AV = D为一个对角矩阵。其中，对角线上的元素是矩阵 A的特征根。令Z(t)=V,Y(t),则Z(t)=V 二AVZ(t) V 二X(t)= DZ(t) V JX(t)这个方程组定义了两个独立的一阶常系数线性微分方程：1Zi(t)= iZ(t) Vi X(t)其中，叫是矩阵A的第i个特征根，V是V-1第i行。乙。)=e J

11、eV1X(t)dt +bie酬再通过变换Z(t) =V Y(t)求Yo二维系统稳定性的一般讨论：对角例子的稳定性性质依赖于对角元的符号。所以，依此类推，非对角系统的稳定性性质依赖于其特征值的符号。于是会产生以下几种可能性：1)两个特征值不同且都是正实数，在这种情况下系统是不稳定的。2)两个特征值不同且都是负实数，在这种情况下系统是稳定的。3)两个特征值是实数但符号相反，在这种情况下系统是鞍点路径稳定的。止匕外，当系统是鞍点路径稳定时，稳定臂对应于与负特征值相关的特征向量。 同理，不稳定臂对应于正特征值相关的特征向量。这里的直观想法仍然是与对角 矩阵相关的轴就是特征向量。正如我们前面的例子中看到

12、的，当系统是对角的时与对角矩阵的负分量相关的轴是稳定臂，与正分量相关的是不稳定臂。4)两个特征值都是负实部的复数，在这种情况下系统以一种振荡方式收敛到 稳态。5)两个特征值都是有正实部的复数，系统是不稳定且振动的。6)两个特征值是有零实部的复数，所示其轨迹是环绕着稳态运动的椭圆。7)两个特征值相等。在这种情况下特征向量矩阵不可逆，所以前面概括的解析解法不适用，此时的解的形式为yi (t) = (bi bi2t)e：t其中bi和源是积分常数和矩阵A中的系数函数。汽是唯一的特征值。若a0,解是不稳定的。更高维系统的稳定性也有类似的性质。如果所有的特征值都为正，这系统是不稳定的。如果所有的特征值都为

13、负，则系统稳定的。如果特征值异号，则系统是鞍点路径稳定的。由于像前面所说的一样，稳定臂对应于与负特征值相关的特征向量，那么稳 定臂的维数就是负特征值的个数。例如在有一个负特征值的3*3系统中，稳定臂(有时被称为稳定流形(stable manifold)是一条通过稳态且对应于这一负特征向 量的直线。如果有两个负特征值，则稳定流形是一个通过稳态的平面。这一平面 由这两个负特征值生成。在一个 n*n系统中，稳定臂是由相关特征向量生成的一 个超平面，其维数等于负特征值的个数。例3,考虑一个非对角系统。yi(t) =0.06yi(t)-y2(t) 1.4y2 (t) = -0.004yl (t) 0.0

14、4以矩阵符号表示，这一系统可被写成jy；(t)M(t),0.060.004一1,i0 大y2(t)J 0.04；它有初始条件为 y1(0)=1, lim.e_0.06t y1 (t) = 0在这个例子中x是一个常数向量X(t)=1.4904对角特征值矩阵D和特征向量矩阵100 0.04,-0.040.10.1/0.140.04/0.14-1/0.141/0.14定义4、=VyBarro的经济增长中译本有误，其100/14 一项的符号为正是错误的。 ,则系统可改写成 ）y2 ）乙=0.1乙 1 01 4Z2 =一0.0钮 9.61 4这样一个由两个微分方程组成的系统，而我们已经知道如何求解。其解

15、为 01tZ1（t） = 100/14 13.1Z2（t） =240/14 bze。04t其中“和b2是要通过边界条件加以确定的积分常数。通过在Z之前乘以V我们可以把对和而言的解转换成对 Y的解y1（t） =10 b4.1t b2e，.4ty2（t） =2 0.04b1e.1t 0.母-04t我们现在需要确定常数b1和b2的值。初始条件y1（0） = 1暗示，b1+b2 = -9。利用条件 lim_e06ty1（t）1=0当t趋于无穷大时中间那个表达式中的第一和第三项都趋于0,但是除非等于0,这意味着“ =0 ,所以，b2 = -9。.因此这个ODE系统的精确解为y1（t） =10 9e04t

16、y2=2 0.9e4t稳定臂为y2（t） =0.1y1（t） 1注意在t=0时，%（0） =1 ,随后持续增加并趋于其稳态值10,变量y2在t=0时等于1.1,随后持续增加并趋于其稳态值 2。换言之，边界条件选择了能使经济 终结于其稳态的的初始值。所选的值使系统位于稳定臂上。在初始点，朝向稳态 的向量为，9 L或者把第一个元素标准化为1单位后，这一向量为负特征向量0.9I1 L因此，稳定臂通过稳态且对应于与负特征值的特征向量。10l撇开解析解，也可以讨论方程组解的性质yi=0.06，(t) 72 1.4 = 0y2(t) - -0.004yi(t) 0.04 = 0的解，稳态解为*yi (t)

17、 =10* 一一 _y2 (t) = 2y； =0的轨迹是向上倾斜的直线y2 =1.4 + 0.06%。如果开始时，y1这条直线 在左边，即y1比均衡值减少一点，y； 0, %是递减的。y =0的轨迹是垂直与 必轴的直线y =10。如果开始时，y1这条直线在左边，即y1比均衡值减少一点，丫20, y2是递增的。反之，如果开始时，y1在这条直线在右边，即y1比均衡值增加一点，y20, y2是递减的。另外，一阶导数等于0产生了两条轨迹，将整个平面分成了四个区域。 可以 综合起来画出在四个区域的箭头走向。为了评价系统的稳定性，我们可以用相位图研究在那个区域系统朝稳定状态 移动，答案是有两个区域系统可

18、能趋于稳定状态，所以，该系统是鞍点稳定的， 在这两个区域，有一条路径，系统沿着它从非稳态开始，在稳态点结束。该路径 叫做稳定臂。在任何非稳定臂上的一点开始，系统都将偏离稳态。在另外两个区 域，有一条路径是通过稳态点的路径叫做非稳定臂。对微分方程系统的图形解基本上也一样.稳定臂和不稳定臂对应于两个特征 向量。如果我们把这两条臂理解为一组新的轴 -也就是说如果我们把旧的轴和线 抹去-那么老的矩阵A就可被表示为对角特征值矩阵.非对角情形下的相位图相应 的看上去就向对角情形下那一个扭曲的版本。例4,新古典经济增长模型的稳定性。考虑拉姆齐模型:处二k(t)尸-凶一k(t)k(t)c(t) = k(t)冷

19、-： c(t)稳定状态所满足的条件为:c(t) = k(t)F -、. k(t)fP+N3)dc . - j=: k 一：dk(二严 k* ,意味着稳定点在消费所以，在k=k* =但i 时，c有最大值。注意k5最大值点的左侧。在k二士I CLI (二.)左侧，消费增长率大于0,在其右侧，消费增长率小于00在c(t) = k(t)R-6,k(t)上侧，资本增长率小于0,在其下侧，资本增长率小 于0。据此，我们可以画出，相位图。五、非线性系统的线性化如同新古典拉姆齐模型那样，经济学中许多最优化问题的解都与非线性常微 分分成(ODE)系统有关。我们可以利用前面讨论过的相位图技巧直接研究稳态 路径的性

20、质，也可以利用泰勒级数展开来对非线性方程进行线性化近似。考虑如下的线性ODE系统：yi(t) = f1yi(t),y2(t),yn(t)y，2(t) = f 2yi(t), y2(t), ,yn(t)yn(t) = fnyi(t), y2(t), ,yn(t)其中，函数1()=12一一，是非线性的。我们可以利用泰勒级数展开来研究系 统在稳态领域中的稳态。一阶展开可被写成yi(t) = f1()fy：()(yi-yi)f；n()(ynTn)Ryn=fn() fy：() (yi -y；)fy：() (yn 7；) Rn其中，是函数f i( ),i =1,2,，n ,在稳态的值，f；j (),i =

21、 1,2,n; j = 1,2,，n是函数在稳态对yi的偏导数。Ri项是泰勒余项。如果系统接近于稳态，那么这些余项很小而且可被忽略。在稳态附近线性化的方便之处在于根据稳态的定义每个方程的第一项1()=1,2,，n为0;也就是说，对所有的i, yi稳态值为0。去掉余项后的线性系统可用矩阵符号表示为*Y = A (Y -丫)其中A是一个对应于稳态的函数的一阶偏导数f；j (),i =1,2，n； j =1,2，n ,构成的nxn常数矩阵。这一线性系统与前面分析过的系统类似。例如5,通过非线性系统线性化讨论下列非线性系统的稳定性。0.3=k-c0.3 k = k-c_0 7一c -c (0.3k .

22、 -0.06)其边界条件为它有初始条件为k(0) = 1 , limb061k(t)=0。稳态值为 t J二二_ _-0 7 一一c = c (0.3k . -0.06)其边界条件为它有初始条件为k(0)=1, limeq6tk(t)=011穆利甘和萨拉伊马丁(1991)提出了一种有效的数值技巧，称之为时间消去方 法(time-elimination method) 0这种方法的关键是从方程中消去时间，就像我们在 构造相位图时所作的那样。稳定臂把 c表示为k的一个函数。在动态规划中，这 一函数有时被称为策略函数。暂时假设我们有这一策略函数的一个闭式解 c=c(k)0时间消去方法为算出策略函数 c=c(k)提供了一种数值技巧。窍门是注意 到这一函数的斜率由对的比率给出：0 7dc _ g _ c _ c(k)0.3 k -0.06= c(k) 0-3dkk k c(k)注意到这是一个关于c的微分方程。其中导数dc/dk是对k而非对t的。为用标 准数值方法解这个方程，我们需要一个边界条件；也就是说我们要知道位于稳定 臂上的一点(c,k)。尽管我们并不了解初始的一对c(0), k(0),但我们清楚策略函