差分方程及matlab求解

1、1 一个半球状雪堆，其体积融化的速率与半一个半球状雪堆，其体积融化的速率与半球面面积球面面积S成正比，比例系数成正比，比例系数k 0。设融化中。设融化中雪堆始终保持半球状，初始半径为雪堆始终保持半球状，初始半径为R且且3小时小时中融化了总体积的中融化了总体积的7/8，问雪堆全部融化还需，问雪堆全部融化还需要多长时间？要多长时间？ 差分方程初步差分方程初步第一节第一节 差分方程的基本概念差分方程的基本概念一、一、 差分的概念差分的概念定义1 设函数yt=f(t)在t=,-2,-1,0,1,2,处有定义,对应的函数值为,y-2,y-1,y0,y1,y2,则函数yt=f(t)在时间t的一阶差分定义为

2、Dyt=yt+1-yt=f(t+1)-f(t)依此定义类推,有Dyt+1=yt+2-yt+1=f(t+2)-f(t+1),Dyt+2=yt+3-yt+2=f(t+3)-f(t+2),一阶差分的性质 (1) 若yt=C(C为常数),则Dyt=0;(2) 对于任意常数k,D(kyt)=kDyt;(3) D(yt+zt)=Dyt+Dzt定义2 函数yt=f(t)在时刻t的二阶差分定义为一阶差分的差分,即 D2yt= D (D yt)= D yt+1- D yt =(yt+2-yt+1)-(yt+1-yt)=yt+2-2yt+1+yt依此定义类推,有D2yt+1= Dyt+2- Dyt+1=yt+3-

3、2yt+2+yt+1,D2yt+2= Dyt+3- Dyt+2=yt+4-2yt+3+yt+2,类推,计算两个相继的二阶差分之差,便得到三阶差分D3yt= D2yt+1- D2yt=yt+3-3yt+2+3yt+1-yt,D3yt+1= D2yt+2- D2yt+1=yt+4-3yt+3+3yt+2-yt+1, 一般地,k阶差分(k为正整数)定义为 这里 ), 3 , 2 , 1()1( )(01111= =- -= =D D- -D D= =D DD D= =D D = =- -+ +- -+ +- - -kyCyyyykiiktikitktktktk)!( !ikikCik- -= =二、

4、二、 差分方程差分方程定义3 含有未知函数yt=f(t)以及yt的差分Dyt, D2yt,的函数方程,称为常差分方程(简称差分方程);出现在差分方程中的差分的最高阶数,称为差分方程的阶. n阶差分方程的一般形式为F(t,yt, Dyt, Dnyt)=0, 其中F是t,yt, Dyt, Dnyt的已知函数,且Dnyt一定要在方程中出现 定义3 含有两个或两个以上函数值yt,yt+1,的函数方程,称为(常)差分方程,出现在差分方程中未知函数下标的最大差,称为差分方程的阶 n阶差分方程的一般形式为F(t,yt,yt+1,，yt+n)=0, 其中F为t,yt,yt+1,，yt+n的已知函数,且yt和y

5、t+n一定要在差分方程中出现.三、三、 差分方程的解差分方程的解定义4 如果将已知函数yt=j(t)代入方程F(t,yt,yt+1,，yt+n)=0,使其对t=,-2,-1,0,1,2,成为恒等式,则称yt=j(t)为方程的解.含有n个任意(独立)常数C1,C2,Cn的解yt=j(t,C1,C2,Cn)称为n阶差分方程的通解.在通解中给任意常数C1,C2,Cn以确定的值所得的解,称为n阶差分方程的特解. 例如,函数yt=at+C(a为已知常数,C为任意常数)是差分方程yt+1-yt=a的通解.而函数yt=at,yt=at-1,均是这个差分方程的特解. 由差分方程的通解来确定它的特解,需要给出确

6、定特解的定解条件.n阶差分方程F(t,yt,yt+1,，yt+n)=0常见的定解条件为初始条件.y0=a0, y1=a1，,yn-1=an-1，这里a0,a1,a2,，an-1均为已知常数 只要保持差分方程中的时间滞后结构不变,无论对t提前或推后一个相同的等间隔值,所得新方程与原方程是等价的,即二者有相同的解.例如,方程ayt+1-byt=0与方程ayt+2-byt+1=0都是相互等价的 四、四、 线性差分方程及其基本定理线性差分方程及其基本定理 形如yt+n+a1(t)yt+n-1+a2(t)yt+n-2+an-1(t)yt1+an(t)yt=f(t) 的差分方程,称为n阶非齐次线性差分方程

7、.其中a1(t),a2(t),an-1(t),an(t)和f(t)都是t的已知函数,且an(t)0,f(t)0.而形如yt+n+a1(t)yt+n-1+an-1(t)yt+1+an(t)yt=0 的 差 分 方 程 , 称 为 n 阶 齐 次 线 性 差 分 方 程 . 其 中ai(t)(i=1,2,，n)为t的已知函数,且an(t)0. 如果ai(t)=ai(i=1,2,n)均为常数(an0),则有yt+n+a1yt+n-1+a2yt+n-2+an-1yt+1+anyt=f(t)， yt+n+a1yt+n-1+a2yt+n-2+an-1yt+1+anyt=0 分别称为n阶常系数非齐次线性差分

8、方程和n阶常系数齐次线性差分方程. 定理定理1(齐次线性差分方程解的叠加原理齐次线性差分方程解的叠加原理) 若y1(t),y2(t),ym(t)是齐次线性差分方程yt+n+a1yt+n-1 +a2yt+n-2+an-1yt+1+anyt=0的m个特解(m2),则其线性组合y(t)=A1y1(t)+A2y2(t)+Amym(t)也是方程的解,其中A1，A2，Am为任意常数定理定理2 n阶齐次线性差分方程yt+n+a1yt+n-1 +a2yt+n-2 +an-1yt+1+anyt=0一定存在n个线性无关的特解定理定理3(齐次线性差分方程通解结构定理齐次线性差分方程通解结构定理) 如果y1(t),y

9、2(t),yn(t)是齐次线性差分方程yt+n+a1yt+n-1 +a2yt+n-2 +an-1yt+1+anyt=0的n个线性无关的特解,则方程的通解为：yA(t)A1y1(t)+A2y2(t)+Anyn(t)， 其中A1，A2，An为n个任意(独立)常数 定理定理4(非齐次线性差分方程通解结构定理非齐次线性差分方程通解结构定理) 如果 (t)是非齐次线性方程yt+n+a1(t)yt+n-1+a2(t)yt+n-2 +an-1(t)yt1+an(t)yt=f(t)的一个特解,yA(t)是其对应的齐次线性方程yt+n+a1yt+n-1 +a2yt+n-2 +an-1yt+1+anyt=0的通解

10、,那么,非齐次线性差分方程的通解为：y(t)=yA(t)+ (t) 即y(t)=A1y1(t)+A2y2(t)+Anyn(t)+ (t)，这里A1，A2,An为n个任意(独立)常数yyy第二节第二节 一阶常系数线性差分方程一阶常系数线性差分方程 一阶常系数线性差分方程的一般形式为yt+1+ayt=f(t) 和yt+1+ayt=0, 其中f(t)为t的已知函数,a0为常数.分别称为一阶常系数非齐次线性差分方程和其对应的齐次差分方程. 一、一、 齐次差分方程的通解齐次差分方程的通解将方程yt+1+ayt=0改写为:yt+1=-ayt, t=0,1,2,假定在初始时刻(即t=0)时,函数yt取任意值

11、A,那么由上式逐次迭代,算得y1=-ay0=-aA,y2=-ay1=(-a)2A,方程的通解为yt =A(-a)t, t=0,1,2,如果给定初始条件t=0时yt=y0,则A=y0,此时特解为：yt =y0(-a)t 二、二、 非齐次方程的通解与特解非齐次方程的通解与特解1. 迭代法求通解迭代法求通解将方程改写为 yt+1=(-a)yt+f(t), t=0,1,2,逐步迭代,则有y1=(-a)y0+f(0),y2=(-a)2y0+(-a)f(0)+f(1),y3=(-a)3y0+(-a)2f(0)+(-a)f(1)+f(2), 由数学归纳法,可得 yt=(-a)ty0+(-a)t-1f(0)+

12、(-a)t-2f(1)+f(t-1) =(-a)ty0+ , (t=0,1,2,)， ty.)1()( )1()1()()0()(1021为为方方程程的的特特解解其其中中 - -= =- - - - - -= =- -+ +- -+ +- -= =tiitttitfatffafayyA(t)=(-a)ty0为对应的齐次方程的通解. 解解例例.2211的的通通解解求求差差分分方方程程tttyy= =- -+ +ttfa2)(,21= = =)12()21(31411)41(12)41(22)21(22)21(211101101101- -= =- - - = = = = = = =- - - -

13、= =- - -= =- - - -= =- - - tttttiittiiittiitity121231)21()12()21(31)21(+ +- -+ += =- -+ += =ttttttAAy方程的通解方程的通解 .32为为任任意意常常数数- -= = AA2.待定系数法求特解待定系数法求特解情形 f(t)为常数方程变为yt+1+ayt=b, a,b均为非零常数试以 (为待定常数)形式的特解代入方程得 +a =(1+a) =b = =ty当a-1时,可求得特解abyt+ += =1当a=-1时,改设特解 (为待定系数),将其代入方程得 (t+1)+a t=(1+a) t+ =b ty

14、t = =求得特解btyt= =方程的通解为 .1 ,1,1)()( 为为任任意意常常数数其其中中AabtAaabaAytyyttAt - -= =+ +- - + + +- -= =+ += =解解例例.521的的通通解解求求差差分分方方程程= =- -+ +ttyy5, 12-= =- - = =ba., 52为为任任意意常常数数AAytt- - = =情形情形 f(t)为为t的多项式的多项式不妨设f(t)=b0+b1t(t的一次多项式),即 yt+1+ayt=b0+b1t, t=1,2,， 其中a,b0,b1均为常数,且a0,b10试以特解 =a+bt,(a,b为待定系数)代入方程得a+

15、b (t+1)+a(a+bt)=b0+b1t，ty上式对一切t值均成立,其充分必要条件是： = =+ += =+ + +10)1()1(babab bb ba a当1+a0时,即a-1时，ababab+ += =+ +- -+ += =1)1(11210b ba a方程的特解为 tabababy+ + + +- -+ += =1)1(11210当a=-1时,改设特解 =(a+bt)t=at+bt2 ty将其代入方程可求得特解211021)21(tbtbby+ +- -= =方程的通解为 - -= =+ +- -+ + + + + +- -+ + +- -= =. 1,21)21(, 1,1)1

16、(1)(21101210atbtbbAatabababaAytt解解例例.231的的通通解解求求差差分分方方程程tyytt+ += =- -+ +2, 3, 110= = =- -= =bba.,22为为任任意意常常数数AttAyt+ + += =情形情形 f(t)为指数函数为指数函数 不妨设f(t)=bdt, b,d均为非零常数,方程变为 yt+1+ayt=bdt, t=0,1,2, 求得特解ttddaby + += =当a+d0时,设方程有特解 =dt, 为待定系数.将其代入方程得 dt+1+adt=bdt, ty当a+d=0时,改设方程的特解 =tdt,为待定系数,将其代入方程可求得特解

17、= =btdt tytyty当a+d=0时,改设方程的特解 =tdt,为待定系数,将其代入方程可求得特解= =btdt ty求得特解ttddaby + += =当a+d=0时,改设方程的特解 =tdt,为待定系数,将其代入方程可求得特解= =btdt ty方程的通解为 = =+ + +- - + + + + +- -= =+ += =. 0,)(, 0,)(dabtdaAdaddabaAyyytttttAt解解例例.21的的通通解解求求差差分分方方程程tttyy= =- -+ +01, 2, 1, 1 = =+ += = =- -= =dadba.,2为为任任意意常常数数AAytt+ += =

18、情形 f(t)为正弦、余弦型三角函数 设f(t)=b1cost+b2sint,其中b1,b2,均为常数,且 0,b1与b2不同时为零.于是非齐次方程变为yt+1+ayt=b1cost+b2sint,a0, t=0,1,2, 设方程有特解 =acost+bsint,a,b均为待定系数. ty将其代入方程得acos(t+1)+bsin(t+1)+aacost+absint =b1cost+b2sint, (acos+bsin +aa)cost+(-asin +bcos +ab)sint=b1cost+b2sint (acos+bsin +aa)cost+(-asin +bcos +ab)sint=

19、b1cost+b2sint 上式对t=0,1,2,恒成立的充分必要条件是 = =+ + + - -= = + + +.)cos(sin,sin)cos(21babab b a a b b a a 其系数行列式 22sin)cos(cossinsincos+ + += =+ +- -+ += =aaaD当D0时,则可求得其解 + + += =- -+ += =;sin)cos(1,sin)cos(11221 b b a ababDbabD当D=(a+cos)2sin2=0时,则有 )(. 1,12. 1,2为整数为整数或或kakak = =+ += = - -= = = 改设特解 .,),sin

20、cos(为待定系数为待定系数b ba a b b a atttyt+ += =代入方程并整理可得 - -= =- -= = = = =.,2121bbbbb ba ab ba a或或方程的通解为 = =+ += =+ + + +- - - -= = =+ + + + + +- -= =. 1,)12(,)12sin()12cos()1(, 1,2),2sin2cos(, 0,sincos)(2121aktkbtkbtAaktkbtkbtADttaAyttt b b a a例例 求差分方程求差分方程yt+1- -2yt= =cost的通解的通解解解 对应齐次方程的通解为对应齐次方程的通解为 yA

21、(t)= =A2t设非齐次方程的特解为设非齐次方程的特解为 = =a acost+b bsint,其中其中a a, b b为待定系数为待定系数 ty将其代入原方程将其代入原方程,并利用三角函数的和角公式并利用三角函数的和角公式,得得 = =- -+ +- -= =+ +- -. 0)21(cos1sin, 11sin)21(cosb ba ab ba a1cos451sin,1cos4521cos- -= =- - -= =b ba a所给方程的通解为所给方程的通解为 ttAyttsin1cos451sincos1cos451cos22- -+ +- - - - = =第三节第三节 二阶常系数

22、线性差分方程二阶常系数线性差分方程 二阶常系数线性差分方程的一般形式为yt+2+a1yt+1+a2yt=f(t)，t=0,1,2,， 其中f(t)为t的已知函数,a1,a2为已知常数,且a20,称为二阶常系数非齐次线性差分方程 特别地,当f(t)0时,方程变为yt+2+a1yt+1+a2yt=0 称为对应的齐次差分方程一、 齐次差分方程的通解 称2a1+a2=0为二阶常系数非齐次线性差分方程或其对应的齐次差分方程的特征方程特征方程它的解(或根)称为方程的特征根特征根(值值) 特征方程的两个根为 )4(),4(2122122112, 1aaaaa- -= =D D- - - -= = (1) 特

23、征根为相异的两实根特征根为相异的两实根当D0时,1, 2为两相异的实根. y1(t)= 1t与y2(t)=2t是齐次差分方程的两个线性无关的特解. 齐次差分方程的通解 ttAAAty2211)( + + = =1,2由特征方程确定,A1,A2为两任意(独立)常数 例例 求差分方程yt+2-7yt+1+12yt=0的通解解解 特征方程为2-7+12=( -3)( -4)=0,有两相异实特征根 1=3, 2=4 原方程的通解为 .,43)(2121为为任任意意常常数数AAAAtyttA + + = =(2) 特征根为两相等的实根特征根为两相等的实根当当D D=0时时, = 1= 2= 为两相等的实

24、根为两相等的实根. 21a- -方程的一个特解：方程的一个特解：yt(t)= t 方程的另一个特解为方程的另一个特解为y(t)=t t，且与且与 t线性无关线性无关. 方程的通解为方程的通解为 .,2)()()()(2112121为为任任意意常常数数其其中中或或AAatAAtytAAtytAtA - -+ += =+ += = 例例 求差分方程求差分方程yt+ +2- -4yt+ +1+ +4yt=0的通解的通解.解解 特征方程为特征方程为 2- -4 + +4=( - -2)2=0，方程有重特征根方程有重特征根 = 1= 2=2 原方程的通解为原方程的通解为yA(t)=(A1+ +A2t)2

25、t, A1,A2为任意常数为任意常数(3) 特征根为一对共轭复根特征根为一对共轭复根当当D D0时时, 1, 2为一对共轭复根为一对共轭复根. 1,2=a aib b=r(cos isin ) .,20 ,tan,21sin,2cos2221为为复复特特征征根根的的辐辐角角为为复复特特征征根根的的模模 a ab b b ba a b b a ararrar = = =+ += =D D= = =- -= = =y1(t)=rtcos t, y2(t)=rtsin t是方程的两个线性无关特解是方程的两个线性无关特解. 方程的通解为方程的通解为yA(t)=rt(A1cos t+ +A2sin t)

26、 其中其中 A1，A2为任意常数为任意常数.例例 求差分方程求差分方程yt+ +2- -2yt+ +1+ +2yt=0的通解的通解解解 特征方程特征方程 2- -2 + +2=( - -1)21=0 特征根为一对共轭复根特征根为一对共轭复根 1,2=1i 方程的通解为方程的通解为 4, 1tan,2 = = = =r.,)4sin4cos(2)( 21212为为任任意意常常数数其其中中AAtAtAtytA + += =二、二、 非齐次方程的特解与通解非齐次方程的特解与通解例例 求差分方程求差分方程yt+ +2- -7yt+ +1+ +12yt=6的通解的通解解解 对应的齐次方程的通解为对应的齐

27、次方程的通解为 yA(t)=A13t+ +A24t，原方程的通解为原方程的通解为yt=yA(t)+ +=A13t+ +A24t+ +1，这里这里A1,A2为任意常数为任意常数 由于由于1+ +a1+ +a2=1- -7+ +120,设特解设特解 =B,B为待定常为待定常数数,将其代入原方程将其代入原方程,求得求得B=1.ty例例 求差分方程求差分方程yt+ +2- -3yt+ +1+ +2yt=4的通解的通解解解 特征方程为特征方程为 2- -3 + +2=( - -1)( - -2)=0,特征根特征根 1=1, 2=2. 对应齐次方程的通解为对应齐次方程的通解为 yA(t)=A1+ +A22

28、t因因1+ +a1+ +a2=1- -3+ +2=0,故应设非齐次方程的特解为故应设非齐次方程的特解为 =Bt,B为待定系数为待定系数,将其代入原方程将其代入原方程,求得求得B=- -4 ty原方程的通解为原方程的通解为yt=yA(t)+ + =A1+ +A22t- -4t，这里这里A1,A2为任意常数为任意常数ty例例 求差分方程求差分方程yt+ +2- -4yt+ +1+ +4yt=3+ +2t的通解的通解. 解解 对应齐次方程的通解为对应齐次方程的通解为yA(t)=(A1+ +A2t)2t此式对此式对t=0,1,2,恒成立的充要条件是恒成立的充要条件是B0- -2B1=3, B1=2.由

29、此解得：由此解得：B0=7，B1=2 设非齐次方程有特解设非齐次方程有特解 =B0+ +B1t,B0,B1为待定系数为待定系数.将其代入原方程中将其代入原方程中,得得(B0- -2B1)+ +B1t=3+ +2t, ty所求非齐次方程的特解为所求非齐次方程的特解为 tyt27+ += =原方程的通解为原方程的通解为 ttAAytytyttA272)()()(21+ + + + += =+ += =A1,A2为任意常数为任意常数 例例 求差分方程求差分方程yt+ +2- -4yt+ +1+ +4yt=5t的通解的通解解解 对应齐次方程的通解为对应齐次方程的通解为yA(t)=(A1+ +A2t)2

30、t设所给非齐次方程的特特为设所给非齐次方程的特特为 =B5t,B为待定系数为待定系数. ty将其代入所给方程将其代入所给方程,可得可得 B5t+ +2- -4B5t+ +1+ +4B5t=5t91= =B非齐次方程的特解为非齐次方程的特解为 tty591 = =所给方程的通解为所给方程的通解为 其中其中A1,A2为任意常数为任意常数tttAtAAytyty5912)()()(21 + + + += =+ += =差分方程在经济学学模型差分方程在经济学学模型一、 存款模型 设St为t期存款总额,i为存款利率,则St与i有如下关系式：St+1=St+iSt=(1+i)St, t=0,1,2,，其中

31、S0为初始存款总额 二、 动态供需均衡模型(蛛网定理) 设Dt表示t期的需求量,St表示t期的供给量,Pt表示商品t期价格,则传统的动态供需均衡模型为： = =+ += =+ += =- -)3(,)2()1(, 111ttttttSDPbaSbPaD其中a,b,a1,b1均为已知常数. (1)式表示t期(现期)需求依赖于同期价格； (2)式表示t期(现期)供给依赖于(t-1)期(前期)价格 (3)式为供需均衡条件 若在供需平衡的条件下,而且价格保持不变,即 Pt=Pt- -1=Pe，静态均衡价格 bbaaPe- - -= =11需求曲线与供给曲线的交点(Pe,Qe)即为该种商品的静态均衡点动

32、态供需均衡模型的等价差分方程 baaPbbPtt- -= =- - -111方程的一个特解 etPbbaaP= =- - -= =11方程的通解为 ettPbbAP+ + = =1若初始价格P0已知时,将其代入通解,可求得任意常数A=P0-Pe,此时,通解改写为 etetPbbPPP+ + - -= =10)(如果初始价格P0=Pe,那么Pt=Pe,这表明没有外部干扰发生,价格将固定在常数值Pe上,即静态均衡如果初始价格P0Pe,那么价格Pt将随t的变化而变化. ,11时时 bbeetetttPPbbPPP= =+ + - -= =+ + + + )(limlim10动态价格Pt随着t的无限增

33、大逐渐地振荡趋近于静态均衡价格Pe.普通商品的价格与供需关系图三、 凯恩斯(Keynes.J.M)乘数动力学模型 设Yt表示t期国民收入,Ct为t期消费,It为t期投资,I0为自发(固定)投资,DI为周期固定投资增量.凯恩斯国民经济收支动态均衡模型为： D D+ += =+ += =+ += =- -)3(,)2()1(,0, 1IIIbYaCICYtttttt(1)式为均衡条件,即国民收入等于同期消费与同期投资之和;(2)式为消费函数,即现期消费水平依赖于前期国民收入(消费滞后于收入一个周期),a(0)为基本消费水平,b为边际消费倾向(0b1);(3)式为投资函数,这里仅考虑为固定投资 在(

34、1)(2)(3)式中消去Ct和It,得到一阶常系数非齐次线性差分方程：Yt-bYt-1=a+I0+DI 方程的一个特解 bIIaYt- -D D+ + += =10方程的通解为 bIIabAYtt- -D D+ + + + = =10其中A为任意常数. 称系数 为凯恩斯乘数. b- -11四、 哈罗德(Harrod.R.H)经济增长模型 设St为t期储蓄,Yt为t期国民收入,It为t期投资,s称为边际储蓄倾向(即平均储蓄倾向),0s1,k为加速系数.哈罗德宏观经济增长模型为：其中其中s,k为已知常数为已知常数11,01 (1)()0(2) (3)tttttttSsYsIk YYkSI-=-=，

35、(1)式表示t期储蓄依赖于前期的国民收入;(2)式表示t期投资为前两期国民收入差的加速,且预期资本加速系数k为常数;(3)式为均衡条件.经整理后得齐次差分方程01= =+ +- - -ttYkskY其通解为ttksAY)1( + += =其中A为任意常数, ,哈罗德称之为“保证增长率” 0 ks其经济意义就是：如果国民收入Yt按保证增长率 增长,那么就能保证t期储蓄与t期投资达到动态均衡,即It=St, t=0,1,2, ks假定t-1期收入Yt-1满足于通解,而t期收入Yt由于某种外部干扰满足), 0()1(称称为为外外部部干干扰扰 + + += =BBksAYtt设B0,那么有 kBSkB

36、sYkBkssABksAkskYYkIttttttt+ += =+ += =+ + += =+ + += =- -= =- - - - -1111)1( )1()(因kB0,故ItSt. 表示:总投资将大于总供给(由储蓄提供),从而对收入产生一个向上的压力,迫使收入较以前增加得更多. 充分地说明了,“保证增长率”保证了国民收入的增长. 五、 萨缪尔森(Samuelson P.A)乘数加速数模型 设Yt为t期国民收入,Ct为t期消费,It为t期投资,G为政府支出(各期均相同).萨缪尔森将乘数和加速数两个参数同时引进而得到国民经济收支均衡模型(也称为乘数-加速数模型)： - -= = = =+ +

37、 += =- - -)3(, 0),()2(, 10,)1(,11kCCkIbbYCGICYtttttttt其中G0为常数,b称为边际消费倾向(常数),k为加速数. 将(2)(3)两式代入(1)并经整理后得：Yt-b(1+k)Yt-1+bkYt-2G 其特解 bGYt- -= =1其经济意义为：国民收入的均衡值等于凯恩斯乘数与政府支出自发投资G的乘积.b- -11对应的齐次方程为 Yt-b(1+k)Yt-1+bkYt-2=0， 其特征方程为 2-b(1+k)+bk=0， 特征方程的判别式 D=b2(1k)2-4bk=bb(1+k)2-4k 当D0时,特征方程有两相异实根 )1(21)1(212

38、1D D+ + += =D D- -+ += =kbkb 齐次方程的通解为：YA(t)=A11t+A22t (A1,A2为任意常数)当D=0时,特征方程有一对相等实特征根 )1(21kb+ += = 齐次方程的通解为： (A1,A2为任意常数为任意常数) tAtAAtY + += =)()(21当D y1=sqh(20,0.0194); y2=sqh(20,-0.0324); y3=sqh(20,-0.0382); round(k,y1,y2,y3)利用plot 绘图观察数量变化趋势n可以用不同线型和颜色绘图nr g b c m y k w 分别表示 红绿兰兰绿洋红黄黑白色： + o * .

39、X s d 表示不同的线型 n plot(k,y1,k,y2,k,y3) 在同一坐标系下画图 plot(k,y2,:) plot(k,y2,-) plot(k,y2,r) plot(k,y2,y) plot(k,y2,y,k,y1,:) plot(k,y2,k,y1,:) plot(k,y2,oy,k,y1,:)用gtext(r=0.0194),gtext(r=-0.0324),gtext(r=-0.0382)在图上做标记。n人工孵化是挽救濒危物种的措施之一，如果每年孵化5只鹤放入保护区，观察在中等自然条件下沙丘鹤的数量如何变化Xk+1=aXk +5 ,a=1+r如果我们想考察每年孵化多少只比

40、较合适，可以令Xk+1=aXk +b ,a=1+rnfunction x=fhsqh(n,r,b)na=1+r;nX=100;nFor k=1:nnX(k+1)=a*x(k)+b;nendnk=(0:20) ; %一个行向量ny1=(20,-0.0324,5); 也是一个行向量nround( k , y 1 ) 对k,y1四舍五入，但 是 不改变变量的值 plot( k , y1) k y1 是行向量列向量都可以也可以观察200年的发展趋势，以及在较差条件下的发展趋势，也可以考察每年孵化数量变化的影响。一阶线性常系数差分方程的解、平衡点及其稳定性一阶线性常系数差分方程的解、平衡点及其稳定性n自

41、然环境下,b=0n人工孵化条件下n令xk=xk+1=x得 差分方程的平衡点nk时，xkx,称平衡点是稳定的0kkxax=10(1)kkkxa xbaa-=+011kkaaxba-=+-1bxa=-1kkxaxb+=+高阶线性常系数差分方程n n 如果第k+1时段变量Xk+1不仅取决于第k时段变量Xk，而且与以前时段变量有关，就要用高阶差分方程来描述一年生植物的繁殖n一年生植物春季发芽，夏天开花，秋季产种，没有腐烂，风干，被人为掠取的那些种子可以活过冬天，其中一部分能在第2年春季发芽，然后开花，产种，其中的另一部分虽未能发芽，但如又能活过一个冬天，则其中一部分可在第三年春季发芽，然后开花，产种，

42、如此继续，一年生植物只能活1年，而近似的认为，种子最多可以活过两个冬天，试建立数学模型研究这种植物数量变化的规律，及它能一直繁殖下去的条件。模型及其求解n记一棵植物春季产种的平均数为c,种子能活过一个冬天的(1岁种子)比例为b,活过一个冬天没有发芽又活过一个冬天的（2岁种子）比例仍为b,1岁种子发芽率a1，2岁种子发芽率a2。n设c,a1,a2固定，b是变量，考察能一直繁殖的条件n记第k年植物数量为Xk，显然Xk与Xk-1,Xk-2有关，由 Xk-1决定的部分是 a1bcXk-1,由Xk-2决定的部分是 a2b(1-a1)bcXk-2 Xk= a1bcXk-1 + a2b(1-a1)bcXk-

43、2 Xk= a1bcXk-1 + a2b(1-a1)bcXk-2n实际上，就是Xk= pXk-1 + qXk-2 我们需要知道x0,a1,a2,c, 考察b不同时，种子繁殖的情况。在这里假设nX0=100,a1=0.5,a2=0.25,c=10,b=0.180.20n这样可以用matlab计算了Xk= a1bcXk-1 + a2b(1-a1)bcXk-2nfunction x=zwfz(x0,n,b)nC=10;a1=0.5;a2=0.25;np=a1*b*c;q=a2*b*(1-a1)*b*c;nx(1)=x0;nx(2)=p*x(1);nfor k=3:nnx(k)=p*x(k-1)+q*

44、x(k-2);nendnk=(0:20);ny1=zwfz(100,21,0.18);ny2=zwfz(100,21,0.19);ny3=zwfz(100,21,0,20);nRound(k,y1,y2,y3)nPlot(k,y1,k,y2,:,k,y3,o),nGtext(b=0.18),gtext(b=0.19),gtext(b=0.20)结果分析：Xk= pXk-1 + qXk-2 (1) x1+px0=0 (2) n对高阶差分方程可以寻求形如的解。代入(1)式得称为差分方程的特征方程。差分方程的特征根：方程(1)的解可以表为C1,c2 由初始条件x0,x1确定。kkx=20pq+=21

45、 , 242ppq-=1122kkkxcc=+1,21,0()kxk n本例中，用待定系数的方法可以求出b=0.18时,c1=95.64, c2=4.36 ， 这样实际上，植物能一直繁殖下去的条件是b0.19112( ,)(0.9430, 0.0430) =-95.64(0.9430)4.36( 0.0430)kkkx =+-1,21,()kxk 1 , 251 02b=线性常系数差分方程组n汽车租赁公司的运营n一家汽车租赁公司在3个相邻的城市运营，为方便顾客起见公司承诺，在一个城市租赁的汽车可以在任意一个城市归还。根据经验估计和市场调查，一个租赁期内在A市租赁的汽车在A,B,C市归还的比例分

46、别为0.6,0.3,0.1;在B市租赁的汽车归还比例0.2,0.7,0.1;C市租赁的归还比例分别为0.1,0.3,0.6。若公司开业时将600辆汽车平均分配到3个城市，建立运营过程中汽车数量在3个城市间转移的模型，并讨论时间充分长以后的变化趋势。0.60.3nA B CnA B CnA B C假设在每个租赁期开始能把汽车都租出去，并都在租赁期末归还0.10.70.20.10.60.30.1模型及其求解n记第k个租赁期末公司在ABC市的汽车数量分别为x1(k),x2(k),x3(k)（也是第k+1个租赁期开始各个城市租出去的汽车数量），很容易写出第k+1个租赁期末公司在ABC市的汽车数量为(k

47、=0,1,2,3)112321233123(1)0.6 ( )0.2( )0.1 ( )(1)0.3 ( )0.7( )0.3( )(1)0.1 ( )0.1( )0.6( )x kx kx kx kx kx kx kx kx kx kx kx k+=+=+=+n用矩阵表示用matlab编程，计算x(k),观察n年以后的3个城市的汽车数量变化情况112233(1)0.60.20.1( )(1)0.30.70.3( )(1)0.10.10.6( )x kx kx kx kx kx k+=+nfunction x=czqc(n)nA=0.6,0.2,0.1;0.3,0.7,0.3;0.1,0.1,0.6;nx(:,1)=200,200,200;nfor k=1:nn x(:,k+1)=A*x(:,k);nend如果直接看10年或者20年发展趋势，可以直接在命令窗口（commond