第3章 差分方程方法

第三章 差分方程方法在实际问题中，许多事物所研究的变量都是离散的形式，所建立的数学模型也是离散的，比如政治、经济和社会等领域中的实际问题。很多时候，即使所建立的数学模型是连续形式，例如象常见的微分方程模型、积分方程模型等，但是往往因为不能求解析解，而需要用计算机求数值解。这要求将连续变量在一定的条件下进行离散化，从而将连续模型转化为离散模型，最后归结为求离散形式的差分方程的问题。关于差分方程研究和求解方法在建立数学模型、解决实际问题的过程中起着重要的作用。3.1 差分方程和常系数线性差分方程3.1.1 差分和差分方程的概念定义 3.1 设函数 , 记为 . 当 取遍所有的非负整数时, 函数)(xfyxy值可以排成一个数列: . 则差 称为函数 的差 ,210 xy1xy分, 也称为一阶差分, 记为 , 即 .xyxy容易验证差分具有如下性质:(1) ;xxyc)(2) .xzz又 xxx xyyyy212 11: )()()( 称为函数 的二价差分, 类似地, 可以定义三阶、四阶差分.x定义 2.2 含有未知函数差分或表示未知函数几个时期值的符号的方程称为差分方程, 如 .0),( ;,21xnxxxxxyyHGF差分方程中所含未知函数差分的最高阶数称为此差分方程的阶. 如果差分方程中关于未知函数及未知函数的各阶差分都是线性函数, 就称此方程为线性的. 如果一个函数代入差分方程后, 方程恒成立, 则这个函数称为该差分方程的解.3.1.2 常系数齐次线性差分方程考虑常系数 阶线性差分方程k. (3.1)021 knnnyayay称代数方程. (3.2)11kkkk为方程(3.1)的特征方程, 特征方程的根称为特征根.常系数线性差分方程的解可根据相应的特征根的情况给出. 下面分别由特征根为单根、重根和复根的情况写出差分方程的解的情况。1.单根的情况设特征方程(3.2)有 个相异的特征根 , 则差分方程kk,21(3.1)的通解为,nknnccy21其中 为任意的常数.kc,212. 重根的情况设特征方程(3.2)有 个相异重根 , 它们的重数l )1(,21kll分别为 , .则差分方程(2.1)的通解为lm,21 liik1.112112lmmminininn li i iyccc3. 单复根的情况设 为特征方程的一对单复根, 它对应的差分方程的解i2,1为( ).ncncsios21arctn,2特别地, 一阶、二阶线性齐次差分方程的通解情况如下:(1)一阶差分方程 的通解为01nay.nca(2) 对于二阶线性差分方程 , 设其特征方程012nnbyy的特征根为 . 21,当 , , 时, 通解为 ;042ba2R21 nncy21当 , 时, 通解为 ;,1anac)(21当 , 时, 通解为 2i2, )sino21cyn(其中 ).arctn3.1.3 常系数非齐次线性差分方程对于常系数非齐次线性差分方程, (3.3)(21 nfyayayknnn 其通解为 , 在这里 为相应的齐次方程的通解, 为方程*nYY*ny(3.3)的一个特解.下面我们给出几个关于一阶、二阶差分方程解的例子.(1) .1nyab通解为 (其中 , 为任意常数), 或 (A1aA1nnbyAa).(2) .1(,)nnyacb当 时, 通解为 ;当 时, 通解为1nncybAab.1nnycA(3) .1mnayc可设方程有特解 , 其中 待定; *01()s mnyBn 01,mB当 时, ; 当 时, . 将这个特解代入方程, 比较两端0sa同次项系数来确定 .()i(4) .21nnyabyc(i)当 时, 特解为 ;0*1ncyab(ii)当 , 且 时, 特解为 ;2a*2ncya(iii)当 , 且 时, 特解为 .1ab(5) .2(,1)nnyycq(i)当 时, 特解为 ;0qab*2nncqyab(ii)当 , 且 时, 特解为 ;220q1*2ncqya(iii)当 , 且 时, 特解为 .20qaba1*4nn(6) 21mnnyyc可设方程有特解 , 其中 待定; *01()s mBn 01,mB当 时, ; 当 , 且 时, ;当 , 10absab2as0ab且 时, .将这个特解代入方程, 比较两端同次项系数来确2as定 .(01)iBm3.2 差分方程的平衡点及其稳定性3.2.1 一阶常系数线性差分方程考虑一阶常系数线性差分方程, (3.4)1nyab当 时, 方程(3.4)的平衡点 可由代数方程 求得, 即1a*yxab.*by如果 , 则称平衡点 为稳定的, 否则就是不稳定的.*limny*y方程(3.4)的稳定性问题可通过变换( )转化为齐次差*nxy分方程 的平衡点 的稳定性问题. 而由差分方程10nxa*0x的解 可知, 其平衡点 是稳定的充要条件为1n()nn*0x. 因此, 方程(3.4)的平衡点 是稳定的充要条件为 .|a*y|1a3.2.2 一阶常系数线性差分方程组考虑一阶常系数线性齐次差分方程组, (3.5)01k ,)(1(Ayk其中 为 维向量, 为 阶常值方阵. 它的平衡点为 .)(ynn0y*方程(3.5)的平衡点 是稳定的充分必要条件是 的所有特* A征根的绝对值小于 ,即 ( ).1|ini,21对于一阶常系数线性非齐次差分方程组 0,k )(BAyk这里, 为 维向量. 其平衡点由线性方程组 给出, 平衡Bn BAy点的稳定性同样为 的所有特征根的绝对值小于 .A13.2.3 二阶常系数线性差分方程考虑二阶常系数齐次线性差分方程.,210 ,212kyaynnn其中 为常数, 其的平衡点 是稳定的充分必要条件是它的特21a*x征方程 0212a的根满足 ( ).1|i,对于二阶常系数非齐次线性差分方程, 其平衡点的稳定性有相同的结果.)2,0( 212 kbyaynnn3.2.4 一阶非线性差分方程考虑一阶非线性差分方程. (3.6)(1nnyf它的平衡点 由方程 求得. 方程(3.6)的近似线性方程为*y. (3.7)()(*1 yfyfnn由于 也是方程(3.7)的平衡点, 因此它们的平衡点有相同的*y稳定性, 即方程(3.6)平衡点 是稳定的充要条件为 .*y*|()|1fy3.3 连续模型的差分方法3.3.1 微分的差分方法假设已知函数 在区间 上的分点)(xf ,ba的函数值, 下面我们来用差商代替该函数bxxan1210在区间内分点的函数微商.1.向前差分;,nkxffxfkk 21 ,)()12.向后差分;,nkxffxfkk 21 ,)()(13.中心差分.,nkxffxfkk 21 ,)()1如果把区间 分为 等分, 步长 , 分点 ,banabh, 我们也常常用二阶差商代替二阶微商:),10( khaxk.,nkhfxffxf kkk 21 ,)(2)(1 3.3.2 定积分的差分方法在这里, 讨论定积分的近似计算问题. 设函数 在区间)(xf上连续, 把区间 分为 等分, 分点 , ba ,ban,10 nkhaxk步长 . 在小区间 上任取点 ( ), 根据定积nh1kxn,21分的定义可得,.nknba fabxf1)(limd)(因此, 我们有如下的求积公式.(1)矩形公式;nkbaxfaxf1)(d)(.nkbaff1)()(2)复化矩形公式.nkkba xfxf1)2(d)(3)梯形公式 )()(2 d)( 11nknkkkba xfbfahxf(4)抛物线公式 1110201()d()4()()6 ,nb kka kknnkkafxfxffxhfbff 其中 , .112()kkxx0,3.3.1 常微分方程的差分方法1.一阶常微分方程考虑一阶常微分方程的初值问题(3.8)0(,)yfx其中函数 连续, 且关于 满足李普希茨条件, 即保证问题(3.8)的()fxy解是存在唯一的.一阶常微分方程的初值问题的数值解通常是将区间 按照0, +)x一定的步长 进行划分, 分点为 , 其中 (h012,x 0nxh), 记 的近似值为 , 建立关于 的递推公式, 再由0,12n ()nyxnyny问题(3.8)的初始条件求出 . 123,(1)欧拉(Euler)法用一阶向前差商近似地代替导数 ,()nyx,11()()n nnyxyh用 代替 , 于是问题(3.8)就化为差分方程的初值问题ny()nx10(,).nnyfxy这就是解一阶常微分方程的初值问题的欧拉法, 其误差为.21()O()nyxh(2)隐式欧拉法和二步欧拉法用一阶向后差商近似地代替导数 , 可得到隐式欧拉法公式()nyx110().nyhf在这个公式中, 没有直接给出 的计算公式, 而是关于 的方程, 1n 1ny因此称为隐式欧拉法. 这种方法的误差也为 .2O()h用一阶中心差商近似地代替导数 , 可得到二步欧拉法公式()nyx102().nyf利用这个公式在计算 时, 需要调用前面两步的信息 和 , 因1n ny1此称之为二步欧拉法. 这种方法的误差为 .3O()h(3)梯形公式和改进的欧拉法将方程 的两端从 到 积分, 得),(yxf nx1 d)()(1nxn xyf利用数值积分中的梯形公式, 并用 和 分别代替 和 , 1n)(ny)(1nx于是可导出如下公式.),(),(211 nnn yxfyxfhy此公式称为梯形公式. 它是欧拉法和隐形欧拉法的算术平均值. 这种方法的误差为 .3O()h由于欧拉法是一个显式算法, 计算量小, 但精度很低 . 梯形公式虽然提高了精度, 但是一种隐式算法, 计算量大. 综合使用这两种方法, 先用欧拉法求得一个初步的近似值, 记作 , 称为预报值; 1ny预报值的精度不高, 但用它替代梯形公式中右端的 , 直接计算, 得到校正值 . 这样建立的预报-校正系统为1ny1 1 (,)(,)2nn nhfxyfxy 预 报校 正这就是改进的欧拉法, 或称为预报-校正法. 这种方法也可以写成下面的平均化形式 ).(21,)(1cpnpncnpyyxhf这种方法的误差为 .3O()h(4)龙格-库塔(Runge-Kutta)法龙格-库塔法的基本思想对于微分方程的初值问题(3.8), 为了求其在 的精确01,nx 值 , 即建立关于精确值的递推关系, 根据拉格朗01(),()nyxyx 日微分中值定理可得.1 1()()(,)nnnnnyxhyxhfyx 记 , 称为区间 上的平均变化率, 则*,Yf1. 问题归结为寻找一个计算 的方法. 由于在实际*1()(nnyxhY *Y的计算中 无法确定, 从而精确的 也就无法确定, 因此总是取*的一个近似值. 例如, 取 , 就得到了欧拉法公式; 取*Y*()nYfxy, 就得到了改进的欧拉法公式.11(,)()2nnfxyfy现在, 我们取区间 上的 个点 ( ), 相应1,nxmnixh1,2m地取 ( ), 用 在这 个点的函数值的加权平均作niyh,2 f为 的近似值 , 即*Y,*1(,)mininiYwfxhy其中 为权系数, 为待定系数. 于是有递推公式iw,i. (3.9)11(,)mnininiyhfxhy实际计算中, 适当选择 , 可使公式(3.9)有较高的精度. 这就iiw是龙格-库塔方法的思想.二阶龙格-库塔公式取区间 的中点 , 又令 , 代1,nx12nxh120,w12入公式(3.9), 得到二阶龙格-库塔公式1221,()0,2.,)nnyhkkfxny这种方法的误差为 .3O()h三阶龙格-库塔公式在区间 中取两个点, 可导出三阶龙格-库塔公式1nx1231213 2(4),6,)0,.(),nnnykkhfxnyfk这种方法的误差为 .4O()h标准四阶龙格-库塔公式(或称为经典公式) 11234213243(),6,)(),().nnnykkkhfxykfxkh这种方法的误差为 .5O()吉尔(Gill)公式11234213 2423()(),6(,),(, ), .nnnykkkkhfxykfxkhy吉尔公式也是一个常用的四阶龙格-库塔公式, 它的误差为 .5O()h2.一阶常微分方程组将前面关于常微分方程中的变量和函数看作向量函数, 相应的差分方法即可用于一阶常微分方程组的情形.考虑二个方程的方程组 0(,) (,.yfxzyzgz令 , 以 表示节点 上的近似值.0, 12nxh ,nnx(1)欧拉法的计算公式 1 0(,), ()nnyhfxyzyxzg(2)隐式欧拉法的计算公式 (0)1(0)111,.nnnnnyhfxyzzgfxyz(3)改进的欧拉法计算公式(0)1 (0)1 11()(,),(,),2,.nnnnnnyhfxyzzgffxyzhzxyzg (4)标准四阶龙格-库塔法计算公式 11234(),6,nhykkzll其中 111 112 22 223 341334133(,), (,), (,),(,), ,().nnn nnn nn nkfxyzlgxyzhkl hklgxyzkfxyzlhklxyhkzl 3.高阶常微分方程对于高阶方程的初值问题, 一般可通过引进新变量, 化为一阶常微分方程组的情形. 譬如, 二阶常微分方程的初值问题 00(,) ,yfxy若令 , 则可化为一阶常微分方程组的初值问题zy00(,), .zfxyz它的龙格-库塔公式为 211234(),6,nnhyzLz这里 1 12 223 241 3(,),),(,4,).nnn nnLfxyzhLzhfxyLhzzL对于一般的高阶常微分方程组, 也有类似的结果.3.4 最优捕鱼问题3.4.1 问题的提出假设鳀鱼分为四个年龄组: 称 1,2,3,4 龄鱼. 各年龄组每条鱼的平均重量分别为 5.07,11.55,17.86,22.99g; 各年龄组鱼的死亡率均为 0.8 条/年; 这种鱼为季节性集中产卵繁殖, 产卵孵化期为每年的后四个月, 平均每条 4 龄与的产卵量为 个, 3 龄鱼的产51.09卵量为 4 龄鱼的一半, 2 龄鱼和 1 龄鱼不产卵. 卵孵化并成活即为1 龄鱼, 孵化成活率(1 龄鱼的条数与总产卵量 之比)为 .n1.20n渔业部门规定, 每年只允许在产卵孵化期前的 8 个月进行捕捞作业, 如果每年投入的捕捞能力固定不变, 即固定努力量捕捞, 这时单位时间捕捞量将与各年龄组鱼群条数成正比, 比例系数称为捕捞强度系数, 通常使用 13mm 网眼的拉网，这种网只能捕捞 3,4 龄鱼, 其两个捕捞系数之比为 .042:1需要解决的问题是:建立数学模型, 分析如何实现可持续性捕捞（即每年开始捕捞时渔场中各年龄组鱼条数不变）, 并且在此前提下得到最高的年收获量（总质量）.3.4.2 模型的假设与符号说明1. 模型的假设(1)只考虑鱼的繁殖和捕捞的变化，不考虑鱼群迁入与迁出;(2)各龄鱼在一年的任何时间都会发生自然死亡;(3)所有鱼都在每年最后四个月内(后 1/3 年)完成产卵孵化的过程, 成活的幼鱼在下一年初成为 1 龄鱼;(4)产卵发生于后四个月之初，产卵鱼的自然死亡发生于产卵之后;(5)相邻两个年龄组的鱼群在相邻两年之间的变化是连续的, 即第年底 龄鱼的条数等于第 年初 龄鱼的条数;ki 1ki(6)4 龄以下的鱼全部死亡;(7)采用固定努力量捕捞意味着捕捞的速率正比于捕捞时各龄鱼群的条数, 因此比例系数即为捕捞强度系数.2. 符号的说明用 表示 时刻（年） 龄鱼的条数; 表示鱼的平均自然死亡()ixttir率, 即 ; 表示 龄鱼的产卵数, 即 , 0.8rif 1234(,)(0,)2Aff; 表示 龄鱼的平均质量, 即51.9Aiw; 表示对 龄鱼的捕捞强度系数, 234(,)(.7,15.86,29)wiqi即 , 其中 为努力捕捞量; 表示产卵开10qE 23始的时刻; 表示对 龄鱼的捕捞量; 表示对 龄鱼的捕捞系数, 即iYi ici.iicx3.4.3 模型的建立与求解1.无捕捞时鱼群的自然增长模型由假设(1)(2)得 (), 12,34; 1, 0,2.iixrtktk 由假设(3)(4)得 11 00.() (),2()xkxk034.Axk由假设(5)(6)得.1 1(0), ()()lim(), 2,3; 0,1iii itkxkxxk2.固定努力量捕捞下鱼群的增长和捕捞模型我们知道, 捕捞期为 , 因此kt(3.10)()(),iiiixrtqExtk(3.11), 1iik(3.12)1(0)()(),2,3iiixxi(3.13)11 00.2(),()kk(3.14)034()(), 1,2.Axxx(1)鱼群的增长规律解方程(3.10)和(3.11), 并注意到连续条件(3.12), 可得(3.15)1()(),12,3iiixkslExk(3.16)1 01.(),2(3.17)034(1)()(),2AxkslExkAslExk其中 , , , .e.49rs210.23e4()eEl(2)捕捞量单位时间对 龄鱼的捕捞量(条数)为i()(), .iytqExitk因此, 第 年全年(前 8 个月)对 龄鱼的捕捞量(条数)为k00()(dd1().iiiii iiqYytqxtslExkr 于是, 第 年全年 总捕捞量(质量)为k34()()().WkwYk(3)可持续性捕捞模型可持续性捕捞, 即意味着由于自然死亡和捕捞使鱼群的数量减少, 但通过产卵繁殖补充, 使得鱼群数量能够在每年开始捕捞时保持不变, 这样的捕捞策略就可以年复一年地一直持续下去. 因此, 可持续捕捞的鱼群数量应是(3.15)(3.16)(3.17)的平衡解, 即模型不依赖于时间 的解, 记平衡解为 . 于是由(3.15)t * (0,1234)ix(3.16)(3.17), 得 *1(), ,iiixslE即 *2*3*2131431034, , ()(),(),