差分方程方法

1、第三章差分方程方法在实际问题中，许多事物所研究的变量都是离散的形式，所建立的数学模型也是离散的，比如政治、经济和社会等领域中的实际问题。很多时候，即使所建立的数学模型是连续形式，例如象常见的微分方程模型、积分方程模型等，但是往往因为不能求解析解，而需要用计算机求数值解。这要求将连续变量在一定的条件下进行离散化，从而将连续模型转化为离散模型，最后归结为求离散形式的差分方程的问题。关于差分方程研究和求解方法在建立数学模型、解决实际问题的过程中起着重要的作用。3.1 差分方程和常系数线性差分方程3.1.1 差分和差分方程的概念定义3.1设函数yf(x),记为yx.当x取遍所有的非负整数时,函数值可以

2、排成一个数列：义1，丫2，,yx，.则差yx1yx称为函数yx的差分，也称为一阶差分，记为yx,即yxyx1yx.容易验证差分具有如下性质：(1) M)cyx；(2) (yxzx)yxzx.又(yx)yx1yx(yx2yx1)(yx1yx)yx22yx1yx：2yx称为函数yx的二价差分，类似地，可以定义三阶、四阶差分.定义2.2含有未知函数差分或表示未知函数几个时期值的符号的方程称为差分方程，如F(x,yx,yx1,yx2,yxn)0；G(x,yx,yx1,yx2,yxn)0；2nH(x,yx,yx,yx,yx)0.差分方程中所含未知函数差分的最高阶数称为此差分方程的阶如果差分方程中关于未知

3、函数及未知函数的各阶差分都是线性函数就称此方程为线性的.如果一个函数代入差分方程后，方程恒成立，则这个函数称为该差分方程的解.3.1.2 常系数齐次线性差分方程考虑常系数k阶线性差分方程YnaiYn1a2Yn2akynk0.(3.1)称代数方程ka1k1akiak0.(3.2)为方程(3.1)的特征方程，特征方程的根称为特征根.常系数线性差分方程的解可根据相应的特征根的情况给出.下面分别由特征根为单根、重根和复根的情况写出差分方程的解的情况。1.单根的情况设特征方程(3.2)有k个相异的特征根1,2,k,则差分方程1 3.1)的通解为nnnYnC11C22Ckk，其中G。，0为任意的常数.2

4、.重根的情况设特征方程(3.2)有1个相异重根1,2，i(1lk),它们的重数1分别为m1,m2,mI,m,k.则差分方程(2.1)的通解为i1m1m2mili1ni1ni1nYnC1,n1GE2LCE1.i1i1i13.单复根的情况设1,2i为特征方程的一对单复根，它对应的差分方程的解nnCicosnC2sinn(22,arctan).特别地，一阶、二阶线性齐次差分方程的通解情况如下：(1)一阶差分方程yniayn0的通解为nynca.(2)对于二阶线fi差分方程yn2ayniby00,设其特征方程2ab0的特征根为1,2.当a24b0,12,1,2R时，通解为ynCi1nC22n；当a24

5、b0,1,2I时，通解为yn(C1c2n)(|)n；当a24b0,1,2i时，通解为ynn(ccosnC2sinn)(其中V22,arctan).3.1.3常系数非齐次线性差分方程对于常系数非齐次线性差分方程ynayn1a2yn2akYnkf(n),(3.3)其通解为ynYnyn,在这里Yn为相应的齐次方程的通解，孔为方程(3.3)的一个特解.下面我们给出几个关于一阶、二阶差分方程解的例子.(1)yn1aynYnyn通解为ynbnA（其中a1,A为任意常数），或-Aan(a1).1aYn1aYncbn(c,b1).当ba时，通解为外*bnAan；当ba时,通解为,n1ncnbAa.mYn1aY

6、ncn可设方程有特解*Ynns(BB1nLBmnm),其中B0,B1,L,Bm待定；当a1时，s0；当a1时，s将这个特解代入方程，比较两端同次项系数来确定Bi（i0,1,L,m).Yn2aYn1bYnc.(i)当1ab0时,特解为Vn(ii)当1ab0,且a2时,c；1ab特解为Vn(6)(iii)当1ab0,且a2时,特解为Y；M;c2一n.2Yn2aYn1bYn(i)当q2aqb(ii)当q2(iii)当aqb2qaqYn2aYn1bYn可设方程有特解cqn(c,q1).0时，特解为0,且2qab0,且2qmcn*Yn*Ynncnqq2aqb0时，特解为*Ynn1cnq.;2qaa0时，

7、特解为*Yn2n1cnq4qans(BBeLBmnm),其中B。,0,L,Bm待定;当1ab0时，s0；当1ab0,且a2时，s1；当1ab0,且a2时，s2.将这个特解代入方程，比较两端同次项系数来确定Bi(i0,1,L,m).3.2差分方程的平衡点及其稳定性3.2.1 一阶常系数线性差分方程考虑一阶常系数线性差分方程yniaynb,(3.4)当a1时，方程(3.4)的平衡点y*可由代数方程xaxb求得，即*by胃如果lim%y*,则称平衡点y*为稳定的，否则就是不稳定的.n方程(3.4)的稳定性问题可通过变换(xnyny*)转化为齐次差分方程XniaXn0的平衡点X*0的稳定性问题.而由差

8、分方程XniaXn0的解Xn(a)%。可知，其平衡点X*0是稳定的充要条件为|a|1.因此，方程(3.4)的平衡点y*是稳定的充要条件为|a|1.3.2.2 一阶常系数线性差分方程组考虑一阶常系数线性齐次差分方程组y(k1)Ay(k)0,k0,1,(3.5)其中y(k)为n维向量，A为nn阶常值方阵.它的平衡点为y*0.方程(3.5)的平衡点y*0是稳定的充分必要条件是A的所有特征根的绝对值小于1,即|J1(i1,2,n).对于一阶常系数线性非齐次差分方程组y(k1)Ay(k)B,k0,1,这里，B为n维向量.其平衡点由线性方程组yAyB给出，平衡点的稳定性同样为A的所有特征根的绝对值小于1.

9、3.2.3 二阶常系数线性差分方程考虑二阶常系数齐次线性差分方程yn2aiyn132Yn0,k0,1,2,.其中ai,a2为常数，其的平衡点x0是稳定的充分必要条件是它的特征方程2ala20的根满足|il1(i1,2).对于二阶常系数非齐次线性差分方程yn2a1yn1a2ynb(k0,1,2,),其平衡点的稳定性有相同的结果.3.2.4一阶非线性差分方程考虑一阶非线性差分方程yn1f(yn).(3.6)它的平衡点y*由方程yf(y)求得.方程(3.6)的近似线性方程为*_*zyn1f(y)(yny)f(y).(3.7)由于y*也是方程(3.7)的平衡点，因此它们的平衡点有相同的稳定性，即方程(

10、3.6)平衡点y*是稳定的充要条件为|f(y*)|1.3.3连续模型的差分方法3.3.1 微分的差分方法假设已知函数f(x)在区间a,b上的分点aX0XiX2Xn1b的函数值，下面我们来用差商代替该函数在区间内分点的函数微商1 .向前差分f(Xk)f(Xk1)f(k),k1,2,n；Xk1Xk2 .向后差分f(Xk)f(Xk)f(k1),k1,2,n；XkXk13 .中心差分f(Xk)f(Xk1)f(k1),k1,2,n.Xk1Xk1如果把区间a,b分为n等分，步长h，分点nXkakh(k0,1,n),我们也常常用二阶差商代替二阶微商：f(Xk)f(Xk1)2f(Xk)f(k1)h1,2,n.

11、4 .3.2定积分的差分方法在这里，讨论定积分的近似计算问题.设函数f(X)在区间a,b上连续，把区间a,b分为n等分，分点Xkakh(k0,1,n),步长hU.在小区间Xk,Xk1上任取点k(k1,2,n),根据定积分的n定义可得，nf(X)dXk).balimf(nnk1因此，我们有如下的求积公式(1)矩形公式baf(x)dx(2)复化矩形公式baf(x)dxnf(Xk).k1bf(x)dx一xk1xkf(Ukjk).2(3)梯形公式f(x)dxba2nnf(Xk1)f(Xk)h2f(a)n1f(b)f(Xk)k1(4)抛物线公式af(X)dxba6n1f(Xk)4f(x1)f(Xk1)k

12、-2f(x1)k一2n12f(Xk),k1其中x1k-21,、一国Xk1),2k0,1,L,n.nbaf(Xk1)；nk13.3.1常微分方程的差分方法1.一阶常微分方程考虑一阶常微分方程的初值问题(3.8)yf(x,y),y(x0)v。,其中函数f(x,y)连续，且关于y满足李普希茨条件，即保证问题(3.8)的解是存在唯一的.一阶常微分方程的初值问题的数值解通常是将区间X0,+)按照一定的步长h进行划分，分点为x0,x1,x2,L,其中xnx0nh(n0,1,2,L),记y(xn)的近似值为yn,建立关于yn的递推公式，再由问题(3.8)的初始条件求出yi,y2,y3,L.(1)欧拉(Eul

13、er)法用一阶向前差商近似地代替导数y(xn),y(xnl)y(xn)y(*ni)丫函)y(xn)；,xn1xnh用yn代替y(xn),于是问题(3.8)就化为差分方程的初值问题yn1ynhf(xn,yn),yoy(x).这就是解一阶常微分方程的初值问题的欧拉法，其误差为2、y(xn1)yn1O(h).(2)隐式欧拉法和二步欧拉法用一阶向后差商近似地代替导数y(xn),可得到隐式欧拉法公式Yn1Ynhf(xn1,Yn1),VoY(X0).在这个公式中，没有直接给出yn1的计算公式，而是关于yn1的方程，因此称为隐式欧拉法.这种方法的误差也为O(h2).用一阶中心差商近似地代替导数y(Xn),可

14、得到二步欧拉法公式yn1yn12hf(Xn,yn),V。Y(X0).利用这个公式在计算Yn1时，需要调用前面两步的信息Yn和Yn1，因此称之为二步欧拉法.这种方法的误差为O(h3).(3)梯形公式和改进的欧拉法将方程yf(x,y)的两端从Xn到Xni积分，得Xn1y(Xn1)y(Xn)f(x,y(x)dxXn利用数值积分中的梯形公式，并用yn和yn1分别代替y(Xn)和y(x0i),于是可导出如下公式hYn1Yn2【f(Xn,yn)f(XniYn1).此公式称为梯形公式.它是欧拉法和隐形欧拉法的算术平均值.这种方法的误差为O(h3).由于欧拉法是一个显式算法，计算量小，但精度很低.梯形公式虽然

15、提高了精度，但是一种隐式算法，计算量大.综合使用这两种方法，先用欧拉法求得一个初步的近似值，记作亍n1,称为预报值；预报值的精度不高，但用它替代梯形公式中右端的ym,直接计算，得到校正值yn1.这样建立的预报-校正系统为预报yn1ynhf(X,yn)h校正yn1yn-f(Xn，yn)f(X1,)这就是改进的欧拉法，或称为预报-校正法.这种方法也可以写成下面的平均化形式ypynhf(Xn,yn),ycynhf(Xn1,yp),1yn1-(ypyJ这种方法的误差为O(h3).龙格-库塔(Runge-Kutta)法龙格-库塔法的基本思想对于微分方程的初值问题(3.8),为了求其在X0,X1,L,Xn

16、,L的精确值y(Xo),y(Xi),L,y(Xn),L,即建立关于精确值的递推关系，根据拉格朗日微分中值定理可得y(Xn1)y(Xn)hy()y(Xn)hf(,y(),XnXn1.记Y*f(,y(),称为区间Xn,Xn1上的平均变化率，则y(Xni)y(Xn)hY*.问题归结为寻找一个计算Y*的方法.由于在实际的计算中无法确定，从而精确的Y*也就无法确定，因此总是取Y*的一个近似值.例如，取Y*f(Xn,yn),就得到了欧拉法公式；取,1Y-f(Xn，Yn)f(Xn1，Yn1),就得到了改进的欧拉法公式.现在，我们取区间Xn,Xn1上的m个点Xnih(i1,2,L,m),相应地取ynih(i1

17、,2,L,m),用f在这m个点的函数值的加权平均作为Y*的近似值，即m一*一,.YWif(Xnih,ynih),i1其中Wi为权系数，i,i为待定系数.于是有递推公式myn1ynhWif(Xn白由).(3.9)i1实际计算中，适当选择i，i，Wi,可使公式(3.9)有较高的精度.这就是龙格-库塔方法的思想.二阶龙格-库塔公式取区间Xn,Xn1的中点X1Xn-h,又令w10,w21,1代n-99入公式(3.9),得到二阶龙格-库塔公式yn1ynhk2,klf(Xn,yn),hk2f(Xn,Vn2n0,1,2,L.2kl)，这种方法的误差为O(h3).三阶龙格-库塔公式在区间Xn,Xn1中取两个点

18、，可导出三阶龙格-库塔公式1k1k2k3Vn1Vn-(k14k2k3),60,1,2,Lhf(Xn,Vn),n一,h1,、hf(Xn-，Vnk)22hf(Xh,Vnk12k2),这种方法的误差为O(h4).标准四阶龙格-库塔公式(或称为经典公式)1八Vn1Vn-(k12k22k3k4),6k1hf(Xn,Vn),h1k2hf(Xn.Vn”1),h1k3hf(Xn产”2),k4hf(Xnh,Vnk3).这种方法的误差为O(h5).吉尔(Gill)公式yn1yn1ki(2.2)k2(2.2)k3kJ62ko,2122-k1-k2),22窘2六人2.kihf(Xn,yn),hk2hf(Xn-,ynh

19、k3hf(xn-，yn2k4hf(Xnh,yn吉尔公式也是一个常用的四阶龙格-库塔公式，它的误差为O(h5).2.一阶常微分方程组将前面关于常微分方程中的变量和函数看作向量函数，相应的差分方法即可用于一阶常微分方程组的情形.考虑二个方程的方程组yf(X,y,z),y(Xo)y0,zg(x,y,z),z(Xo)4.令XnX0nh,n1,2,L,以yn,zn表示节点Xn上的近似值.(1)欧拉法的计算公式Yn1Ynhf(Xn,yn,zn),y0y(X0)zn1znhg(Xn,Yd,4),Zo2(X0)(2)隐式欧拉法的计算公式yn0)1Ynhf(Xn,Yn,Zn),率Znhg(Xn,Yn,Zn),y

20、n1ynhf(Xn1,yn1,zn1),(0)(0)、Zn1Znhg(Xn1,Yn1,Zn1).(3)改进的欧拉法计算公式(0)Yn1Ynhf(xn,yn,zn),Z*Znhg(Xn,Yn,Zn),h(0)(0)1Yn1Yn-f(Xn,Yn,Zn)f(Xn1,Yn1,Zn1),2h(0)(0)、iZn1Zn-g(Xn,Yn,Zn)g(Xn1,Yn1,Zn1).2(4)标准四阶龙格-库塔法计算公式h-Yn1Yn(k12k22k3k4),6h一一一Zn1Zn。1222lD，6其中k1f(Xn，Yn,Zn),I19西$2),k2hhk1f(Xnc，Ync，Zn22hl1Jg(Xn2h2,Ynhk12

21、，ZnhLJk3hhk2f(Xn-,Yn,Zn22hl2)2，l3g(Xnh2，Ynhk22，Zn当2，k4f(Xn1,Ynhk3,ZnM3),l4g(Xn1,Ynhk3,Znhk).3.高阶常微分方程对于高阶方程的初值问题，一般可通过引进新变量，化为一阶常微分方程组的情形.譬如，二阶常微分方程的初值问题Yf(X，Y，Y),Y(%)Y0,Y(X)Y0,若令zy,则可化为一阶常微分方程组的初值问题zf(x,y,z),yz,y(X0)Y0,z(x0)Y0.它的龙格-库塔公式为h2Yn1YnhZn(L1L2L3),6h_一一Zn1Zn(L12L22L3L4),6这里L2f(Xn,ynZ),hL3f(

22、Xn二，y2hf(Xn-,y2hZn2hZnhL1,Zn-2h2L1），f(Xn1,ynhZn2h12,ZnhL3).对于一般的高阶常微分方程组，也有类似的结果.3.4最优捕鱼问题3.4.1 问题的提出假设鳏鱼分为四个年龄组：称1,2,3,4龄鱼.各年龄组每条鱼的平均重量分别为5.07,11.55,17.86,22.99g;各年龄组鱼的死亡率均为0.8条/年；这种鱼为季节性集中产卵繁殖，产卵孵化期为每年的后四个月，平均每条4龄与的产卵量为1.109105个,3龄鱼的产卵量为4龄鱼的一半，2龄鱼和1龄鱼不产卵.卵孵化并成活即为1龄11鱼，孵化成活率（1龄鱼的条数与总产卵量n之比）为1.220.1

23、.221011n渔业部门规定，每年只允许在产卵孵化期前的8个月进行捕捞作业，如果每年投入的捕捞能力固定不变，即固定努力量捕捞，这时单位时间捕捞量将与各年龄组鱼群条数成正比，比例系数称为捕捞强度系数，通常使用13mm网眼的拉网，这种网只能捕捞3,4龄鱼，其两个捕捞系数之比为0.42:1.需要解决的问题是：建立数学模型，分析如何实现可持续性捕捞（即每年开始捕捞时渔场中各年龄组鱼条数不变)，并且在此前提下得到最高的年收获量(总质量).3.4.2模型的假设与符号说明1 .模型的假设(1)只考虑鱼的繁殖和捕捞的变化，不考虑鱼群迁入与迁出；(2)各龄鱼在一年的任何时间都会发生自然死亡；(3)所有鱼都在每年

24、最后四个月内(后1/3年)完成产卵孵化的过程成活的幼鱼在下一年初成为1龄鱼；(4)产卵发生于后四个月之初，产卵鱼的自然死亡发生于产卵之后；(5)相邻两个年龄组的鱼群在相邻两年之间的变化是连续的，即第k年底i龄鱼的条数等于第k1年初i1龄鱼的条数；(6)4龄以下的鱼全部死亡；(7)采用固定努力量捕捞意味着捕捞的速率正比于捕捞时各龄鱼群的条数，因此比例系数即为捕捞强度系数.2 .符号的说明用为表示t时亥(年)i龄鱼的条数；r表示鱼的平均自然死亡率，即r0.8；%表示i龄鱼的产卵数，即(工/)(0,0,-A,A),A1.109105；Wi表示i龄鱼的平均质量，即(W1,W2,W3,W4)(5.07,

25、11.55,17.86,22.99)；qi表示对i龄鱼的捕捞强度系数，即(qqQQ)(0,0,0.42E,E),其中E为努力捕捞量；表示产卵开3始的时刻；Yi表示对i龄鱼的捕捞量；G表示对i龄鱼的捕捞系数，即CiY3.4.3模型的建立与求解1.无捕捞时鱼群的自然增长模型由假设(1)(2)得&rxi(t),i1,2,3,4;ktk1,k0,1,2,L由假设(3)(4)得X1(k1)X)(k由假设(5)(6)得1.2210111.221011x0(k)X0()A)1x3(k)Ax4(k).2.Xi(0)X,xMk1)xi(k1)tlirri1xi(t),i1,2,3;k0,1,2,L.固定努力量捕

26、捞下鱼群的增长和捕捞模型我们知道，捕捞期为ktk，因此*rxi(t)qi(E)xi(t),ktk,(3.10)&rxi(t),ktk1,(3.11)x(0)xi,xi1(k1)xi(k1),i1,2,3,1221011x1(k1)12-x0(k1.2210x0(k),A%(k)-x3(k)Ax4(k),k0,1,2,L.(3.12)(3.13)(3.14)(1)鱼群的增长规律解方程(3.10)和(3.11),并注意到连续条件(3.12),可得x(k1)sli(E)xi(k),i1,2,3,(3.15)1221011x1(k1)1.22ii10%(k),(3.16)1.221011x0(k)(3

27、.17)Axo(k1)2s13(E)x3(k)AsL(E)x4(k),0.42EEe,14(E)e其中ser0.4493,11(E)E)1,h(E)(2)捕捞量单位时间对i龄鱼的捕捞量(条数)为Yi(t)qi(E)xi(t),kt因此，第k年全年(前8个月)对i龄鱼的捕捞量(条数)为丫(k)yi(t)dtq(E)xi(t)dtqi(?、(1sli(E)x(k).00rqi(E)于是，第k年全年总捕捞量(质量)为W(k)w3K(k)w4K(k).(3)可持续性捕捞模型可持续性捕捞，即意味着由于自然死亡和捕捞使鱼群的数量减少，但通过产卵繁殖补充，使得鱼群数量能够在每年开始捕捞时保持不变，这样的捕捞策略就可以年复一年地一直持续下去.因此，可持续捕捞的鱼群数量应是(3.15)(3.16)(3.17)的平衡解，即模型不依赖于时间t的解，记平衡解为x*(i0,1,2,3,4).于是由(3.15)(3.16)(3.17),得*x1即*x2s*,xsx2*bx0*XiT,XobX