1 大学数学实验 Experiments in Mathematics 实验2 差分方程和数值微分 清华大学数学科学系 差分方程 ~ 离散时段上描述变化过程的数学模型 ? 一年期存款年利率为 r，存入 M，记第 k年本息为 x k Mxkxrx kk ==+= + 01 ,,2,1,0,)1( L n年后本息为 Mrx n n )1( += ? 污水处理厂每天将污水浓度降低比例 q, 记第 k天 的污水浓度为 c k , L,2,1,0,)1( 1 =?= + kcqc kk 离散动态过程（系统），实际的变化可以是连续的 0 )1( cqc n n ?= )1lg( 2lg q n ? ?= 天后污水浓度降低一半 1.一阶线性常系数差分方程 2. 高阶线性常系数差分方程 3. 线性常系数差分方程组 4. 非线性差分方程 数值微分简介 ? 建立离散动态过程的数学模型; ? 用 MATLAB计算数值解; ? 作理论分析 (平衡点及其稳定性 ). 差分方程 例 1 濒危物种 (Florida 沙丘鹤 )的自然演变 和人工孵化 一阶线性常系数差分方程 ? 在较好自然环境下，年平均增长率为 1.94% ? 在中等自然环境下，年平均增长率为 -3.24% ? 在较差自然环境下，年平均增长率为 -3.82% 如果在某自然保护区内开始有 100只鹤，建立 描述其数量变化规律的模型，并作数值计算 . 生态学 家估计 如果每年人工孵化 5只鹤放入该保护区， 在中等自然环境下鹤的数量将如何变化 ? 模型及其求解 例 1 濒危物种 (Florida 沙丘鹤 )的自然演变 和人工孵化 记第 k年沙丘鹤的数量为 x k ， 自然环境下年平均增长率为 r L,2,1,0,1, 1 =+== + kraaxx kk 0 5 10 15 20 40 60 80 100 120 140 160 r=0.0194 r=-0.0324 r=-0.0382 设每年人工孵化的数量为 b， L,2,1,0, 1 =+= + kbaxx kk 0 5 10 15 20 40 60 80 100 120 140 r=-0.0324,b=5 r=-0.0324,b=0 结果分析 例 1 濒危物种 (Florida 沙丘鹤 )的自然演变 和人工孵化 时间充分长后 (k→∞ )沙丘鹤数量的变化趋势 L,2,1,0, 0 == kxax k k a>1(r>0)时 x k →∞ , a<1(r<0)时 x k → 0 自然环境下 ,1, 1 raaxx kk +== + 在中等及较差的自然环境下沙丘鹤将濒于灭绝。 L,2,1,0, 1 1 0 = ? ? += k a a bxax k k k 人工孵化条件下 baxx kk += +1 a<1(r<0)时 x k → x=b/(1-a) 0 50 100 150 200 100 110 120 130 140 150 160 r=-0.0324,b=5 x=5/0.0324=154.32 2 一阶线性常系数差分方程的平衡点及其稳定性 差分方程的一般形式 L,2,1,0, 1 =+= + kbaxx kk 差分方程的平衡点 ~代数方程 x=ax+b 的根 x=b/(1-a) L,2,1,0),1/( =?+= kabcax k k 差分方程的解 c= x 0 -b/(1-a)由初始值 x 0 和 a、 b确定 若 k→∞时 x k → x, 平衡点 x稳定 , 否则平衡点 x不稳定 平衡点稳定的充要条件是 ?a ?<1 L,2,1,0, 1 1 0 = ? ? += k a a bxax k k k 高阶线性常系数差分方程 例 2 一年生植物的繁殖 一年生植物春季发芽，夏天开花，秋季产种。 没有腐烂、风干、被人为掠去的那些种子可以活过冬 天，其中的一部分能在第二年春季发芽，然后开花、 产种，其中的另一部分虽未能发芽，但如又能活过一 个冬天，则其中一部分可在第三年春季发芽，然后开 花、产种，如此继续。 一年生植物只能活一年，且近似地认为，种子最 多可以活过两个冬天 。 建立数学模型研究植物数量的变化规律， 及它能够一直繁殖下去的条件 . 模型及其求解例 2 一年生植物的繁殖 设一棵植物平均产种数为 c, 种子能够活过冬天 的比例为 b, 活过冬天的那些种子在来年春季发芽的 比例为 a 1 ，未能发芽的那些种子又活过一个冬天的 比例仍为 b, 在下一年春季发芽的比例为 a 2 。 x k ~第 k年的植物数量 设今年种下 (并成活 )的数量为 x 0 011 bcxax = += ?11 kk bcxax 记 p= -a 1 bc, q= -a 2 b(1-a 1 ) bc L,3,2,0,0 2101 ==++=+ ?? kqxpxxpxx kkk L,3,2,0)1( 212 ==? ? kbcxaba k 寻找形如 x k =λ k 的解 设 c=10， a 1 =0.5， a 2 =0.25， b=0.18, 0.19, 0.20， x 0 =100 模型及其求解例 2 一年生植物的繁殖 0 2 =++ qpλλ 0 21 =++ ?? kkk qxpxx 特征方程 特征根 2 4 2 2,1 qpp ?±? =λ L,2,1,0, 2211 =+= kccx kk k λλ 差分方程 常数 c 1 , c 2 由 x 0 , x 1 确定 差分方程的解 b 2 305 2,1 ± =λ kk k xb )043.0(36.4)943.0(64.95,18.0 ?+== kk k xb )0477.0(35.4)0477.1(65.95,2.0 ?+== ?λ 1, 2 ?<1时 x k →0(k→∞) ?λ 1, 2 ? >1时 x k →∞ (k→∞) 植物能够一直繁殖下去的条件为 b>0.191 例 3 蛛网模型的推广 )( 00 xxyy kk ??=? α )( 001 yyxx kk ?=? + β ) 2 ( 0 1 02 y yy xx kk k ? + =? + + β 生产者的管理水平和素质提高 蛛网模型 产量由前两时段平均价格决定 L,2,1,)1(22 012 =+=++ ++ kxxxx kkk αβαβαβ 02 2 =++ αβαβλλ 4 8)( 2 2,1 αβαβαβ λ ?±? = 平衡点 x=x 0 x k → x 0 （ k→∞）的条件是 ?λ 1, 2 ?<1 2/ 2,1 αβλ = 2<αβ经济趋向稳定的条件 原蛛网模型稳定的条件 1<αβ 高阶线性常系数差分方程的平衡点及其稳定性 LL ,2,1, 11110 ==++++ +??++ kbxaxaxaxa knknnknk 0 1 1 10 =++++ ? ? nn nn aaaa λλλ L λ 1 , λ 2 , …λ n LL ,2,1, 2211 =++++= kxcccx k nn kk k λλλ c 1 , …,c n 由初始值 x 1 , …,x n 确定 )/( 110 nn aaaabx ++++= ? L 特征方程 特征根 平衡点 差分方程的解 平衡点稳定的条件 : 所有特征根的模小于 1 3 线性常系数差分方程组 例 4 汽车租赁公司的运营 汽车租赁公司在 3个相邻的城市运营，在一个 城市租赁的汽车可以在任意一个城市归还 . 在 A市租赁在 A, B, C市归还的比例分别为 0.6, 0.3, 0.1 在 B市租赁在 A, B, C市归还的比例分别为 0.2, 0.7, 0.1 在 C市租赁在 A, B, C市归还的比例分别为 0.1, 0.3, 0.6 公司开业时将 600辆汽车平均分配到 3个城市， 建立运营中汽车数量在 3个城市间转移的模型， 讨论时间充分长以后的变化趋势 . 例 4 汽车租赁公司的运营 模型及其求解 0.1 0.2 0.1 0.1 0.3 0.3 0.7 0.60.6 A B C ? ? ? ? ? ++=+ ++=+ ++=+ )(6.0)(1.0)(1.0)1( )(3.0)(7.0)(3.0)1( )(1.0)(2.0)(6.0)1( 3213 3212 3211 kxkxkxkx kxkxkxkx kxkxkxkx x 1 (k), x 2 (k), x 3 (k)~第 k 个租赁期末公司在 A, B, C市的汽车数量 x(k)=[x 1 (k), x 2 (k), x 3 (k)] T ? ? ? ? ? ? ? ? ? ? = 6.01.01.0 3.07.03.0 1.02.06.0 A L,2,1,0),()1( ==+ kkAxkx 120120120121121123125130140160200x 3 (k) 300300300300300299297294284260200x 2 (k) 180180180180179179178176176180200x 1 (k) 109876543210k 开始时 600辆汽车平均分配到 3个城市 开始时 600辆汽车全分配给 A市 12012012011911811611310590600x 3 (k) 3003003003002992972922812521800x 2 (k) 180180181181183187195214258360600x 1 (k) 109876543210k 例 4 汽车租赁公司的运营 模型及其求解 时间充分长后 3个城市的汽车数量趋向稳定，稳定值与初始分配无关 例 4 汽车租赁公司的运营 模型及其求解 设稳定值为 x L,2,1,0),()1( ==+ kkAxkx xAx = ? ? ? ? ? ? ? ? ? ? = 6.01.01.0 3.07.03.0 1.02.06.0 A 矩阵 A的一个特征根 λ=1， 且 x是对应的特征向量。 A各列之和均为 1 A有特征根 λ=1 x可由 Ax=x, x 1 +x 2 +x 3 =600 算出 : x 1 =180, x 2 =300, x 3 =120 例 5 按年龄分组的种群增长 ? 动物因自然或人工繁殖而增加，因自然死亡 和人为屠杀而减少 ; ? 不同年龄动物的繁殖率、死亡率有较大差别 ; ? 将种群按年龄等间隔地分成若干个年龄组， 时间离散化为时段 ; ? 在稳定环境下假定各年龄组种群的繁殖率和 死亡率与时段无关 ; ? 建立按年龄分组的种群增长模型 ; ? 预测未来各年龄组的种群数量 ; ? 讨论时间充分长以后的变化趋势 . ? 种群按年龄大小等分为 n个年龄组，记 i=1,2,…n ? 时间离散为时段，长度与年龄组区间相等，记 k=1,2,… ? 第 i 年龄组 1雌性个体在 1时段内的繁殖率为 b i ? 第 i 年龄组在 1时段内的死亡率为 d i , 存活率为 s i =1-d i 1,,2,1),()1( 1 ?==+ + nikxskx iii L x i (k)~时段 k第 i 年龄组的种群数量 )()1( 1 1 kxbkx i n i i ∑ = =+ (设至少 1个 b i >0) 例 5 按年龄分组的种群增长 模型及其求解 4 例 5 按年龄分组的种群增长 )()1( kLxkx =+ )0()( xLkx k = T n kxkxkxkx )](),(),([)( 21 L= ~按年龄组的分布向量 ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? = ? ? 00 0 000 1 2 1 121 n nn s s s bbbb L O M O L ~Leslie矩阵 (L矩阵 ) 模型及其求解 1,,2,1 ),()1( 1 ?= =+ + ni kxskx iii L )()1( 1 1 kxbkx i n i i ∑ = =+ ∑ = = n i i kxkxkx 1 )(/)()( ~ 归一化 : 例 5 按年龄分组的种群增长 模型及其求解 设一种群分成 n=5个年龄组 , 繁殖率 b 1 ~b 5 = 0, 0.2, 1.8, 0.8, 0.2, 存活率 s 1 ~s 4 = 0.5, 0.8, 0.8, 0.1, 各年 龄组现有数量均为 100 ． 0 5 10 15 20 25 30 0 0.2 0.4 0.6 0.8 0 5 10 15 20 25 30 0 100 200 300 400 500 x 1 (k)~x 5 (k) 的图形 （自上而下 ) 的图形 （自上而下 ) )( ~ ~)( ~ 51 kxkx )0()( xLkx k = 结果分析 例 5 按年龄分组的种群增长 nk k L,3,2, 1 =≤ λλ? L矩阵存在正单特征根 λ 1 ， ? 若 L矩阵存在 b i , b i+1 >0, 则 nk k ,,3,2, 1 L=< λλ T n n ssssss x ? ? ? ? ? ? = ? ? 1 1 121 2 1 21 1 1* ,,,,1 λλλ L L特征向量 * 1 )( lim cx kx k k = ∞→ λ , c是由 b i , s i , x(0)决定的常数且 时间充分长后 x(k), 的稳定性)( ~ kx )0()( xLkx k = Leslie矩阵的性质 结果分析例 5 按年龄分组的种群增长 时间充分长后 x(k), 的稳定性)( ~ kx 按年龄组的分布向量趋向稳定分布 (与 x(0)无关 ) * 1 )( lim cx kx k k = ∞→ λ )()1()2 kxkx λ≈+ ~各年龄组种群数量按同一倍 数 λ (固有增长率 ) 增减 )()1( kLxkx =+ 与基本模型 比较 3） λ=1时 * )()1( cxkxkx ≈≈+ ~ 各年龄组种群 数量不变 * x xkx ~ )( ~ )1 ≈ (归一化的特征向量 ) * x ~ 1个个体在整个存活 期内的繁殖数量为 1 1 121121 =+++ ?nn sssbsbb LL 3） λ=1时 ** xLx = [] T n ssssssx 121211 * ,,,1 ? = LL ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? = ? ? 00 0 000 1 2 1 121 n nn s s s bbbb L O M O L 结果分析例 5 按年龄分组的种群增长 T n ssssx ],,,,1[ 1211 * ? = L,)()4 * xckx k λ≈ ~存活率 s i 是同一时段的 x i+1 与 x i 之比 （与 s i 的定义 比较） )()1( 1 kxskx iii =+ + 1,2,1),()( 1 ?=≈ + nikxskx iii L 非线性差分方程 例 6 离散形式的阻滞增长模型 (Logistic模型 ) )1()( N x rxtx ?=& t→∞, x→N (与 r大小无关 ) 离散形式 x k ~某种群第 k代的数量 (k=0,1,2,… ) 指数 增长 模型 连续形式 阻滞 增长 模型 rxtx =)(& kkk rxxx =? +1 k k kk x N x rxx )1( 1 ?=? + 对于不同的 r，研究 k →∞ 时 x k 的趋势 线性方程 非线性方程 5 0 10 20 30 40 0 0.2 0.4 0.6 0.8 1 1.2 r=0.3 0 10 20 30 40 0 0.2 0.4 0.6 0.8 1 1.2 1.4 r=1.8 0 10 20 30 40 0 0.2 0.4 0.6 0.8 1 1.2 1.4 r=2.5 设 N=1，取 r=0.3, 1.8, 2.5, 初值 x 0 =0.1， 例 6 离散形式的阻滞增长模型 (Logistic模型 ) k k kk x N x rxx )1( 1 ?=? +x k ~某种群第 k代的数量 x k 单调地趋向 N x k 振荡地趋向 Nx k 不收敛 结果分析例 6 离散形式的阻滞增长模型 L,2,1,0,)1( 1 =?+= + kx N x rxx k k kk k k kk x N x rxx )1( 1 ?=? + x N x rxx )1( ?+= 差分方程的平衡点 x=N, x=0 的根 : 若 x k → x (k→∞) , 则平衡点 x稳定 ; 否则， x不稳定 . ? 平衡点 x=0 不稳定 ; ? 平衡点 x=N 的稳定性与 r 有关 . 研究非线性差分方程平衡点稳定的条件 L,2,1,0),( 1 == + kyfy kk 非线性差分方程平衡点稳定的条件 平衡点 ~代数方程 y=f(y) 的根 y * 非线性差分方程 : L,2,1,0,))(( *** 1 =+?′= + kyyyyfy kk f 在 y * 点作 Taylor展开，保留线性项 近似线性方程 : 1)( * <′ yf y * 对近似线性方程是稳定平衡点 , 对非线性方程也是稳定平衡点 . 1)( * >′ yf y * 对近似线性方程是不稳定平衡点 , 对非线性方程也是不稳定平衡点 1, )1( += + = rbx Nr r y kk 例 6 离散形式的阻滞增长模型 k k kk x N x rxx )1( 1 ?+= + 变量和 参数代换 )1( 1 kkk ybyy ?= + 平衡点 x=N, x=0 平衡点 y=1-1/b, y=0 )1()( ybyyf ?= )21()( ybyf ?=′ bbf ?=?′ 2)/11( 1)0( >=′ bf y=0不稳定 平衡点 y * =1-1/b稳定的条件 : )2(31 <<< rb 若 b>3 (即 r>2), 平衡点 y * 不稳定 1)( * <′ yf 例 7 寄主 —寄生模型 ? 寄主靠自然资源为生，假定其数量用离散形式的 阻滞增长模型描述 ; ? 寄生物的存在会减少其增长，寄生物数量越多， 寄主的增长率减少得越多，假设寄主的减少率与 寄生物数量成正比 ; ? 寄生物完全靠寄主为生，假定其相邻两代数量之 比与寄主数量成正比 . 建立寄主 —寄生模型研究二者数量变化的规 律，讨论时间充分长以后的趋势 . 例 7 寄主 —寄生模型 模型及其求解 x k ~第 k代寄主的种群数量 (最大容量 N，固有增长率 r); y k ~第 k代寄生物的种群数量（ k=0, 1, 2, …） . ??=? + k k kk x N x rxx )1( 1 kkk ybxy = +1 a~寄生物由寄主处攫取营养，阻滞寄主增长的能力 b~寄主供养寄生物，使寄生物增长 (或减少 )的能力 非线性差分方程组 kk yax 6 例 7 寄主 —寄生模型 模型及其求解 kkk k kk yaxx N x rxx ??=? + )1( 1 kkk ybxy = +1 设 N=100， r=1.5 寄生物在最好的条件下每代的数量可以翻一番， 即 x k =N时 y k+1 =2y k ， bN=2， b=2/N=0.02; 设 a=0.025。 0 20 40 60 80 100 0 20 40 60 80 100 x k y k → 50 → 30 例 7 寄主 —寄生模型 模型及其求解 3个平衡点： )) 1 1(, 1 ( bNa r b ? kkk k kk yaxx N x rxx ??=? + )1( 1 kkk ybxy = +1 0 20 40 60 80 100 0 20 40 60 80 100 x k y k N=100， r=1.5, b=0.02, a=0.025 =(50,30) 寄主 —寄生物共存的平衡点 (x, y) = (0, 0), (N, 0)， 数值微分 用离散方法近似计算函数 y=f(x)在某点 x=a的导数值 h afhaf af )()( )( ?+ ?′ h hafaf af )()( )( ?? ?′ 前差公式 后差公式 h hafhaf af 2 )()( )( ??+ ?′ 中点公式 )()( 2 )()()( 3 2 hOaf h afhafhaf ±′′+′±=± 误差为 O(h) 误差为 O(h) 误差为 O(h 2 ) 泰勒展开 : 函数 y=f(x)在等间距 h的分点 x 0 <x 1 <… <x n 上用离散数值表示为 y 0 , y 1 , … , y n 数值微分 1,2,1, 2 )( 11 ?= ? ?′ ?+ nk h yy xf kk k L h yyy xf h yyy xf nnn n 2 34 )( 2 43 )( 12210 0 +? =′ ?+? =′ ?? ， 三点公式，误差为 O(h 2 ) 在中间点 x 1 , … , x n-1 在两端点 x 0 , x n 实验练习 实验目的 1.练习用差分方程建立离散动态过程的数学模型， 并用 MATLAB计算其数值解； 2.了解差分方程平衡点及其稳定性的基本知识； 3. 练习数值微分的计算。 实验内容 ? 老师课上布置 ? 同学通过 “清华网络学堂 ”提交作业 ? 网址： http://learn.tsinghua.edu.cn