差分方程方法

1、ii=1第四章差分方程方法在实际中，许多问题所研究的变量都是离散的形式，所建立的数学模型也是离散的，譬如，像政治、经济和社会等领域中的实际问题。有些时候，即使所建立的数学模型是连续形式，例如像常见的微分方程模型、积分方程模型等等，但是，往往都需要用计算机求数值解。这就需要将连续变量在一定条件下进行离散化，从而将连续型模型转化为离散型模型，因此，最后都归结为求解离散形式的差分方程解的问题。关于差分方程理论和求解方法在数学建模和解决实际问题的过程中起着重要作用。下面就不同类型的差分方程进行讨论。所谓的差分方程是指：对于一个数列，把数n列中的前n+1项x(i二0,1,2,n)关联起来所得到的方程。i

2、41常系数线性差分方程常系数线性齐次差分方程常系数线性齐次差分方程的一般形式为x+ax+ax+.+ax=0(4.1)n1n-12n-2kn-k其中k为差分方程的阶数，a(i=1,2,k)为差分方程的系数，且a丰0(kn)。对应的ik代数方程九k+a九k-1+a九k-2+a=0(4.2)12k称为差分方程的(4.1)的特征方程，其特征方程的根称为特征根。常系数线性齐次差分方程的解主要是由相应的特征根的不同情况有不同的形式。下面分别就特征根为单根、重根和复根的情况给出差分方程解的形式。特征根为单根设差分方程(4.1)有k个单特征根九，九,九，,九，则差分方程(4.1)的通解为TOC o 1-5 h

3、 z123kx=C九n+C九n+C九n，n1122kk其中c,c，,c为任意常数，且当给定初始条件12k HYPERLINK l bookmark10 o Current Document x=X(0)(i=1,2,k)(4.3)ii时，可以唯一确定一个特解。特征根为重根设差分方程(4.1)有l个相异的特征根九，九，九，,九Glk)重数分别为123lm1，m2，mi且m=k则差分方程()的通解为x二区cni-1九n+芳cni-1九nF瓦Cni-1九nn1i12i2lili=1i=1i=1同样的，3.由给定的初始条件（4.3）可以唯一确定一个特解。特征根为复根设差分方程（4.1）的特征根为一对共

4、轭复根九，九=aip和相异的k-2个单根12九，九,九，则差分方程的通解为34kx=cpncosn0+cpnsinn0+c心+c心HFc九n,n123344kki其中i其中p=尹2+p2,9=arctanP同样由给定的初始条件（4.3）可以惟一确定一个特解。另外，对于有多个共轭复根和相异实根，或共轭复根和重根的情况，都可以类似地给出差分方程解的形式。412常系数线性非齐次差分方程常系数线性非齐次差分方程的一般形式为(4.4)x+ax+axFFax(4.4)n1n-12n-2kn-k其中k为差分方程的阶数，川二1,2，,为差分方程的系数，ak丰-If（n）为已知函数。在差分方程（4.4）中，令f

5、（）=0,所得方程x+ax+axFFax=0（4.5）n1n-12n-2kn-k称为非齐次差分方程（4.4）对应的齐次差分方程，即与差分方程（4.1）的形式相同。求解非齐次差分方程通解的一般方法为首先求对应的齐次差分方程（4.5）的通解x*，然后求非齐次差分方程（4.4）的一个n特解x（o），贝ynx=x*+x（o）nnn为非齐次差分方程（4.4）的特解。关于求x*的方法同求差分方程（4.1）的方法相同。对于求非齐次方程（4.4）的特解nx（0）的方法，可以用观察法确定，也可以根据f（n）的特性用待定系数法确定，具体方法可n参照常系数线性非齐次微分方程求特解的方法。4.2差分方程的平衡点及其稳

6、定性一般来说，差分方程的求解是困难，实际中往往不需要求出差分方程的一般解，而只需要研究它的平衡点及其稳定性即可。一阶线性常系数差分方程一阶线性常系数差分方程的一般形式为x+ax=b,k=0,1,2,k+1k其中a,b为常数，它的平衡点由代数方程x+ax二b求解得到，不妨记为x*.如果limx=x*，则称平衡点x*是稳定的，否则是不稳定的。kT8*为了便于研究平衡点x*的稳定性问题，一般将其转化为求方程xk+1+叮0的平衡点x*=0的稳定性问题。事实上，由x+ax-0k+1k可以解得x=Cax,k0于是x*=0是稳定的平衡点的充要条件是：|a|Y1.4.2.2一阶线性常系数差分方程组一阶线性常系

7、数齐次差分方程组的一般形式为x(k+1)+Ax(k)=B,k=0,1,2,.其中x(k)为口维向量，A为nxn阶常数矩阵。它的平衡点x*=0是稳定的充要条件是A的所有特征根都有卜|1(i=l,2,n).i对于一阶线性常系数非齐次差分方程组x(k+1)+Ax(k)=B，k=0,1,2,.的情况同样给出4.2.3二阶线性常系数差分方程二阶线性常系数齐次差分方程的一般形式为x+ax+ax=0,k=0,1,2,k+21k+12k其中a,a为常数，其平衡点x*=0是稳定的充要条件是特征方程九2+a九+a二0,1212的根九,九2满足卜J1，卜J1。对于一般的x+ax+ax二0的平衡点的稳定性问题同样给出

8、。类似地，也可直k+21k+12k接推广到n阶线性差分方程的情况。4.2.4一阶非线性差分方程一阶非线性方程的一般形式为x=f(x)，k=0,1,2,，k+1k其中f为已知函数，其平衡点定义为方程x=f(x)的解x*。-f(x*)(-x-f(x*)(-x*)+f(*),kx则x*也是一阶线性万程xk+1已f(x*)(-x*)+f(x*)的平衡点，故此，平稳衡点x*稳定kk+1的充要条件是|f(x*)1。4.3连续模型的差分方法微分的差分方程已知f(x)在点x处的函数值f(x)k=0,1,.,n+1)，且kka=xx.x=b，试求函数的导数值f(x),(k=1,2,.n)。01n+1k根据导数的

9、定义，用差商代替微商，则有下面的差分公式。向前差：f(、f(x)f(x)(k12)f(x)沁屮k(k=1,2,.,n),kx-xk+1k向后差：f(xk)-f向后差：f(xk)-f(x)-f(x)kx-xkk-1(k=1,2,.,n),中心差：f(xk)-f(x)-f(x)(k=1,2,.,n),kk-1定积分的差分方法已知函数f(X)在点x处的函数值f(x),(k=1,2,.,n),且在a,b上可积，试求函kk数在a,b上的积分值ibf(x)dx。a根据定积分的定义，则有一般的求积公式ibf(x)dxibf(x)dxakkk=0其中A为求积系数，它与x的选取方法有关。取不同的求积系数，可以得

10、不同的求积公式。kkba对于等距节点x=a+kh(k=0,1,.,n)，其中步长h=为很小的数，则有如kn下的求积公式。复化矩阵公式；2)复化梯形矩阵；2)复化梯形矩阵；Jbf(xJbf(x)dxqh艺f(x)+f(x)=a2kk+1k=0复化辛普森(Simpson)矩阵公式；ibf(x)dxQa+2刃f(x)+f(b)kk=1呢f(xk)+4f(xk)+f(xk+1)k=0k+2=f(a)+逻f(xk+1盪f(f(b)k=0k2k=1其中xk+其中xk+12(4)=2(xk+x)为子区间x,x(k=0,1,.,n1)的中点。k+1kk+1复化柯特斯(Cotes)公式;(1)a+k+hI2丿J

11、ibf(x)dxqh艺fak=032艺f(xk+3)+艺f(7f(b)k=0k4k=1丄.32艺f(xk+3)+艺f(7f(b)k=0k4k=1丄.f(x丄)+k+k+4k=02k=0其中x,x,x为子空间x,x(k二0,1,.,n-1)中的四等分点。,Jkk+1k+k+k+2444.3.3常微分方程的差分方法一阶常微分方程的差分方法设一阶常微分方程的定解问题为y=f(x,y),y(x)=y,00其中函数f(x,y)关于y满足李普希兹条件，即保证问题解的存在唯一性。现在的问题是求方程在一系列节点xx.x.处的近似数值解12ny,y,y,不妨假设步长为h二x-x为常数。在此，我们根据微分的差分方

12、法，12nn+1n即用差商来近似代替微商，再利用“步进式”方法，可以给出求解问题(4-6)的差分方法。单步欧拉(Euler)公式用差商竺*丄近似代替y(x)二f(x,y(x)中的导数，则可以得差分公式hnnny二y+hf(x,y),n二0,1,2,.n+1nnn其精度为O(h2)阶的。两步欧拉公式用差商兰化+1二近似代替y(x)二f(x,y(x)中的导数，则可得差分公式2hnnny二y+2hf(x,y),n二0,1,2,.n+1n-1nn两步法需要用到前两步的方信息，一般不能自行起步，需先用单步方法求出y，其精度是O(h2)阶的。梯形公式对于方程y二f(x,y)的两边在x,x上求积分得nn+1

13、y(x)二y(x)二y(x)+ixn+if(x,y(x)dx.n+1nnxn利用积分的差分方法中梯形公式求解积分,y(x)n+1h,y(x)n+1x+1f(x,y(x)dx沁-f(x,y(x)+f(x2nnn+1hy(x)-y(x)+-f(x,y(x)+f(x,y(x)n+1n2nnn+1n+1离散化即可得到微分方程的梯形差分公式-y=y+-f(x,y)+f(x,y),N=0,1,2,n+1N2nNn+1N+1这是一个隐式格式，计算量大，一般不单独使用。其镜的也是O(-3)阶的。改进的欧拉公式由于单步欧拉公式色精度低，但计算量小；矩形公式精度高，但是计算量大，为此我们综合运用这两种方法舅老爷得

14、到改进的欧拉公式，其精度为O3)阶的。预报：yn+yn+1,y),n=0,1,nn校正：yn+1,ynn,yyn+1,ynn,yxin+1-n+1丿,n=0,1,或写成平均化形式；cnncnn+1py=三(y+y),n=0,1,2,.n+12pc龙格库塔法龙格库塔方法的基本思想：y(x)-y(x)对于微分方程的定解问题(4.6),考虑差商八n+1：丿n，根据阿格朗日微分中值-定理可得gx，n+1y(x)=y(x)+-Y*。现在n+1ny(x)gx，n+1y(x)=y(x)+-Y*。现在n+1nn+1nn记Y*二f(g,y(g)的平均变化率,nn+1的问题只要找到寻找一种好的计算Y*的近似方法。

15、如果取Y*沁f(x,y)=Y，则就是欧拉公式。nn1如果取Y*沁f(x,y)+fC,y力=Y，则相应的就是改进的欧拉公式。2nnn+1n+12现在，我们取m个点(x+ah,y+Ph)gEx,x(=1,2,m)，用f在这m个nininn+1点的函数值的加权平均作为Y*的近似值，即33艺艺f(x+ah,y+i=10h)ni其中为权系数。则有in+1+h艺f(xii=1+ah,y+0n+1+h艺f(xii=1+ah,y+0h)ini4.7)其中a，i为待定系数。i实际上，适当选择a，i0，使得公式有更咼的精度，这是龙格库塔方法的思i想。二阶龙格-库塔公式在lx,x内取中点xnn+1=xn+211+-

16、h，则可取巴=0,笃=1，a=卩=-代人(4.7)22式得到二阶龙格-库塔公式，其精度为O(h3阶y=y+hY,n+1n2Y=f(x,y),1nnhhjn2n21丿y2=fX七,y+：Y1,n=0,1,2,三阶龙格-库塔公式：在lx,x内任取二点x=x+ph，x=x+qh，(0pq1)类似的方法nn+1n+pnn+qn可得到三阶的龙格-库塔公式其精度是O(4)阶的常用的是三阶的情况。y=y+-(Y+4Y+Y)n+1n6123=f(x,y),nn/hh)=fIx+,y+Y(n2n21丿=f(x+h,y+hCY+2Y)nn12,n=0,1,2,四阶龙格库塔公式：类似的方法可以得到四阶龙格-库塔公式

17、，其精度是O(h5y=y+-(Y+2Y+2Y+Y)TOC o 1-5 h zn+1n61234Y=f(x,y),nn“/hh“)小-宀Y=fx+,y+Y,n=0,1,2,丿In2n2i丿+h+hY22Y=f(x,y+hY)4n+1n3一阶常微方程组的差分方法将前面的单个方程中的变量和函数视为向量，相应的差分方法即可用于由多个方程组的一阶方程组的情形。对于二个方程的方程组00(4.8)y=f(x.y,zz=g(X,y,z00(4.8)设以y,z表示函数在节点x=x+nh,n=1,2,上的近似解，则有改进的欧拉公式:nnn0预报:校正:y=yn+1nn+1hy,预报:校正:y=yn+1nn+1hy

18、,=y+2n+1n2,y,znnn,y,znnnf(x,y,zn+1,yn+1,zn+1g(xg(x,y,znnn(gX,yn+1n+1z_n+1丿四阶龙格-库塔公式y=y+-(Y+2Y+2Y+Y)n+1n61234C1c,n=0,1,2,z=z+-(Z+2Z+2Z+Z)n+1n61234其中Y=f(xZ=gY=f(xZ=g(x,y,z;nnn,y,znnn);hhh、 HYPERLINK l bookmark125 o Current Document n2n21n21丿hhh、n2n21n21丿hhh、n2n22n22丿hhh、n2n22n22/Y=fa,y+hY,z+hZ)4n+1n3n

19、3Z=gvc,y+hY,z+hZ丿4n+1n3n3其他的公式也都可以类似得到，即相当于同时求解多个一阶方程，从方法上没有本质的差别。3.高阶常微分方程的差分方法对于某些高阶方法的定解问题，原则上可以转化为一阶方程组来求解。譬如，对于如下的二阶微分方程的定解问题Y=fx+,y+Y,z+Z2Z=gx+,y+Y,z+Z2ccfVY=fx+,y+Y,z+-Z3Z=gx+,y+-Y,z+-Z3丄)y=fW,y,y丿y(x)=y,y(x)=y0000若令z二y，则可以化为一阶方程组的定解问题z=fC,y,z)(4.9)Vy=z,y(x)=y,y(x)=y0000(4.9)实际上,(4.9)式可以视为(4.

20、8)式的特例，类似地可以得到相应的求解差分公式4.4最优捕鱼问题4.3.1问题的提出假设鯷鱼可分为4个年龄组：称1、2、3、4龄鱼。各年龄组每条鱼的平均重量分别为5.07，11.55，17.86，22.99(g)；各年龄组鱼的自然死亡率均为0.8(1/年)；这种鱼为季节性集中产卵繁殖，产卵孵化期为每年的最后4个月，平均每条4龄鱼的产卵量为1.109*105(个)，3龄鱼的产卵量为这个数的一半，2龄和1龄鱼不产卵。卵孵化并成活为1.22x10111龄鱼，成活率(1龄鱼条数与产卵量n之比)为.22x1011+n渔业部门规定，每年只允许在产卵孵化期前的8个月内进行捕捞作业。如果每年投入的捕捞能力固定

21、不变，即固定努力量捕捞，这时单位时间捕捞量将与各年龄组鱼群条数成正比比例系数称为捕捞强度系数。通常使用13mm网眼的拉网，这种网只能捕捞3、4龄鱼，其两个捕捞系数之比为0.42：1.要解决的问题是：建立数学模型，分析如何实现可持续性捕捞(即每年开始捕捞时渔场中各年龄组鱼群条数不变)，并且在此前提下得到最高的年收获量(总质量)。模型的假设与符号说明模型的假设只考虑鱼的繁殖和捕捞的变化，不考虑鱼群迁入与迁出；各龄鱼在一年的任何时间都会发生自然死亡；所有鱼都在每年最后四个月内(后1/3年)完成产卵孵化的过程，成活的幼鱼在下一年初成为1龄鱼；产卵发生于后四个月之初，产卵鱼的自然死亡发生与产卵之后；相邻

22、两个年龄组的鱼群在相邻两年之间变化是连续的，即第k年底i龄鱼的条数等于第k+1年初i+1龄鱼的条数；4龄以上的鱼全部死亡；采用固定努力量捕捞意味着捕捞的速率正比于捕捞时各龄鱼群的条数，比例系数为捕捞强度系数。符号的说明用x(t)表示t时刻(年)i龄鱼的条数；r表示鱼的平均自然死亡率，即r=0.8；fiiA表示龄i鱼的产卵数，即(f,f,f,f)=(0,0,A)，A=1.109x105；w表示龄i鱼群TOC o 1-5 h z12342i的捕捞强度系数，即(w,w,w,w)二(5.07,11.55,17.86,22.99)；q表示i龄鱼群的捕捞1234i-2强度系数，即(q,q,q,q)二(0,

23、0,0.42e,e)，e为捕捞努力量；t=表示产卵开始的12343Y月份；Y表示i龄鱼的捕捞量；C表示i龄鱼的捕捞率，即C。iiixi模型的建立与求解无捕捞时鱼群的自然增长模型由假设(1)和(2)得dx(/dx(/)idt=-rxC),ii=1,2,3,4；ktk+1，k=0,1,2,又由假设(3)和(4)得x(kx(k+1)=11.22*10111.22*10ii+xk+-x(k+tA=x(k+t+Ax(k+10(丿23(丿4(丿由假设(5)和(6)得x(0)=x,i1limx(t)tk+11x(Kx(0)=x,i1limx(t)tk+11I+1i+(I=1,2,3;k=0,1,)2。固定努

24、力量捕捞鱼群的增长和捕捞模型由假设知，捕捞期为ktk+1，则有则有则有dxC)-i=-rxdtI(t)-q(E)xC),dxC)-i=-rxdtI(t)-q(E)xC),kIItk+t-(4.10)dx(t)i-dt-rx(t)k+ttk+1(4.11)x(0)=xI1x(k+1)=x(k+1)，I=1,2,3I+1I+(4.12)x(k+1)=11.22*1011(1.22*1011+xk+1(4.13)(A(A(-)k+t=xk+1+Ax1k+1V丿23V丿4V丿x0k（1）鱼群的增长规律求解方程（1）鱼群的增长规律求解方程（4.10）和（4.11），并利用连续条件（4.12）式可得x(k

25、+1)=sl(E)X(k),I=1,2,3,I+1II(4.15)x(k+1)=1Bx(k+1)=1Bb+x(k)0(k),(4.16)x(k+1)=0-s-I(E)xAs-tl(E)x44(4.17)=1.22x1011。其中s=e-r=0.4993,I(E)=I(E)=1,/(E)=e-0.421-,A=1.109=1.22x1011。123（2）捕捞量单位时间第i龄鱼的捕捞量（条数）为yC）=q（E）xC）,ktk+1,iii第k年全年（8个月）第i龄鱼的捕捞量（条数）为()(1-sH(Ex(k)iY(k)=Jty(t)dt=f()(1-sH(Ex(k)ii0i0ii于是，第k年总捕捞量

26、（质量）为W（k）=wY（k）+wY（k）3344（3）可持续性捕捞模型可持续捕捞，即意味着由于自然死亡和捕捞使鱼群减少，而通过产卵繁殖补充，使得鱼群能够在每年初开始捕捞时保持平衡不变，这样的捕捞策略就可以年复一年地一直持续下去。因此，可持续捕捞的鱼群数应是（4.15）、（4.16）、（4.17）式的平衡解，即模型不依赖于时间t的解x*（i=0,1,2,3,4）。求解（4.15）、（4.16）、（4.17）得ix*i+1=si(E)x*，i=1,2,3x*i+1iix*2xx*2x*0将（4.18）式代入 HYPERLINK l bookmark179 o Current Document =

27、six*,x*=sx*=s2x*,x*=si(E)x*=s3i(E)x*132143331bx*=0,b+bx*0=0.5As-ti(E)x*+Ast-i(E)x*33444.20)式得代入（4.19）式有求解可得8x代入（4.19）式有求解可得8x*=0.5+si(E)As31(E)x0*31x*18b0.5+si(E)As31(E)x*=43，8b+0.5+si(E)As3i(E)x*431x*x*1(1)=巾-B(E)代入（4.18）式得到x*x*2(1)=sb1-Be)J，x*3x*3x*=s3l（Eb431、B（E）其中B（E）=0.5+si（E）As31（E）.43当B（E）1时，x*0即意味着捕捞过度，致使鱼群灭绝。当B（E）=1时,0条数）为E=31.4称之为过度捕捞力量，因此，可以在EE的范围内找最优捕捞策略。0条数）为在可持续性捕捞的条件下，第i龄鱼的年捕捞量qiY=qiY=i（E（1s-l（E）；+q（E）ix*ii=3,4。整个鱼群的年捕捞量（重量）为Y（E）=wy+wy33440.42