第四章种群数量的状态转移-微分方程 人口数量,物种数量,网络系统元件的运行状态 状态转移:变化率(即:某个物理量导数)。 微分方程模型:导数与各个变量之间的某种平衡关系 本章内容:微分方程(组)的模型建立、数值解 图形解;用 Matlab几何直观地展示各种求解方法的求 解结果 §4-1人口问题 马尔萨斯( Malthus),英1766-1834关于人口 数量的指数增长模型: 基本假设:人口增长率(即:单位时间内增长的人口 数量与总人口数量之比)为常数。 模型建立:t=0时刻的人口数量为已知,记为P, 人口增长率常数记为k,人口数量p是时间t的函数, p=p(),在短时间t→t+M内,人口增长数量为 p(t+△)-p(1)=k·p(t) 得 d()=k-p(,p(0)=P dt 解得 p(t)=poe (这就是马尔萨斯人口模型 模型求解:为了确定模型中的参数,不妨借用美 国从1790年至1840年的人口数据: 年份: 179018001810182018301840 人口/103:39295308724096381286617069 p(1) h(p(t)/p0) Po I-Io 代入。=1790,1=39291=1840.,p(1)=17069得 k=0.0294 p()=3929c01. 模型检验:用此结果来预测美国从1850年之后的 人口数量,再与实际人口数量比较。 年份:185018601870188018901900 预测值:229003080041300554007430099700 实际值:231923143338558501566294875995 相对误差:1.13%2.12%706%104%18%31% 年份:19101920193019401950 预测值:133800179500240900323200433700 实际值:91972105711122755131669150697 相对误差:45%70%96%145%188% 年份:19601970198019902000
第四章 种群数量的状态转移----微分方程 人口数量，物种数量，网络系统元件的运行状态。 状态转移：变化率（即：某个物理量导数）。 微分方程模型：导数与各个变量之间的某种平衡关系。 本章内容：微分方程（组）的模型建立、数值解、 图形解；用 Matlab 几何直观地展示各种求解方法的求 解结果。 §4—1 人口问题 马尔萨斯（Malthus），英 1766--1834 关于人口 数量的指数增长模型： 基本假设：人口增长率（即：单位时间内增长的人口 数量与总人口数量之比）为常数。 模型建立： 0 t = t 时刻的人口数量为已知，记为 0 p ， 人口增长率常数记为 k ，人口数量 p 是时间 t 的函数， p = p(t) ，在短时间 t →t + t 内，人口增长数量为 p(t + t) − p(t) = k p(t)t ， 得 0 0 ( ), ( ) ( ) k p t p t p dt dp t = = ， 解得 ( ) . ( ) 0 0 k t t p t p e − = （这就是马尔萨斯人口模型） 模型求解：为了确定模型中的参数，不妨借用美 国从 1790 年至 1840 年的人口数据： 年份： 1790 1800 1810 1820 1830 1840 人口/ 3 10 ：3929 5308 7240 9638 12866 17069 . ln( ( ) ) , ( ) 0 ( ) 0 0 0 t t p t p e k p p t k t t − = = − 代入 t 0 =1790, p0 = 3929,t =1840, p(t) =17069 得 k = 0.0294. ( ) 3929 . 0.0294( −1790) = t p t e 模型检验：用此结果来预测美国从 1850 年之后的 人口数量，再与实际人口数量比较。 --------------------------------------------------------------------- 年份： 1850 1860 1870 1880 1890 1900 预测值：22900 30800 41300 55400 74300 99700 实际值：23192 31433 38558 50156 62948 75995 相对误差：1.13% 2.12% 7.06% 10.4% 18% 31% --------------------------------------------------------------------- 年份： 1910 1920 1930 1940 1950 预测值：133800 179500 240900 323200 433700 实际值： 91972 105711 122755 131669 150697 相对误差： 45% 70% 96% 145% 188% --------------------------------------------------------------------- 年份： 1960 1970 1980 1990 2000
预测值:582000780900104770014058001886300 实际值:179323203212226505248710281416 相对误差:225%284%363%465%570% 该模型作长期预测不合理。 §4-2微分方程的数值解法 略 §43微分方程图解法 1.斜率场表示微分方程的解 给定=f(x,y),则过平面直角坐标系上任意 点(x0,%),都可以按照斜率f(x,yo)做一条短直线 这样的短直线布满整个坐标平面,形成的图形称为 f(x,y)的斜率场 斜率场的构造方法 (1)明确矩形区域 (2)沿着横向(x方向)、纵向(y方向)分别以小 步长h,h2分割该区域 (3)过每个结点分别做短直线,如过结点(x,y) 做直线y-y=f(x,yx-x):若(x,y)≤2/则 代入x=x+1解得y=y+hf(x,y),连接两点 (x,y)(x+1,y)即可:否则代入y=ym解得 x=x.+ f(x,’连接两点(x,y)(xy)即可。 例1:画出 dr sn xsin y F 的斜率场。 解:考虑在区域-9≤x,y≤9上,以步长h=来 分割,x=-9+i,y=-9+九f(x,y)= sin xsin J 程序清单 clear for =1: 40 y=9+(-1)h; f=sin(x)*sin(y); if abs(f yly +h*f
预测值：582000 780900 1047700 1405800 1886300 实际值：179323 203212 226505 248710 281416 相对误差：225% 284% 363% 465% 570% --------------------------------------------------------------------- 该模型作长期预测不合理。 §4—2 微分方程的数值解法 略。 §4—3 微分方程图解法 1．斜率场表示微分方程的解 给定 f (x, y) dx dy = ，则过平面直角坐标系上任意一 点 ( , ) 0 0 x y ，都可以按照斜率 ( , ) 0 0 f x y 做一条短直线。 这样的短直线布满整个坐标平面，形成的图形称为 f (x, y) 的斜率场。 斜率场的构造方法： （1）明确矩形区域； （2）沿着横向（x 方向）、纵向（y 方向）分别以小 步长 1 2 h ,h 分割该区域； （3）过每个结点分别做短直线，如过结点 ( , ) i j x y 做直线 ( , )( ) j i j i y − y = f x y x − x ：若 2 1 f (xi , y j ) h h 则 代 入 = i+1 x x 解 得 ( , ) j 1 i j y = y + h f x y ，连接两点 ( , ),( , ) 1 x y x y i j i+ 即可；否则代入 = j +1 y y 解 得 ( , ) 2 i j i f x y h x = x + ，连接两点 ( , ),( , ) i j j+1 x y x y 即可。 例 1：画出 x y dx dy = sin sin 的斜率场。 解：考虑在区域− 9 x, y 9 上，以步长 40 18 h = 来 分割， x ih y jh f x y x y i j = −9 + , = −9 + , ( , ) = sin sin . 程序清单： clear hold on h=18/40; for i=1:40 x=-9+(i-1)*h; for j=1:40 y=-9+(j-1)*h; f=sin(x)*sin(y); if abs(f)<=1 y1=y+h*f;
plot([x, x+h][y,yl]) else plot([x, x1]y y+hD end end hold off 读书P45中下段 上机练习:中=y-2x,05x的斜率场 2.相平面轨迹表示微分方程的解 关于两个函数x(),y(1)的微分方程组 x()=8(xy) ()=g2(x,y) x(0)=x0,y(t0)=y 如何直观地表示x与y的关系呢? 法一:相平面轨迹图 求出数值解,得到一系列的离散解 再在X0Y面上连接这些点。如书P48图43 法二:利用斜率场作图 ￠82(x,y) dx8(x,y了(x,y),做斜率场即可。如书P48 图44 §4-4 Matlab软件求解 1.解析解 (1)一阶微分方程 求分=1+y2的通解: dsolve(Dy=1+y2x) d x 求2=1+x2-y的通解: dsolve(Dy=1+x^2-y =1+y的特解: dsolve(Dy=1+y^2y(0=1,x) 求{dx y(0)=1 (2)高阶微分方程 求解{xy++(-)y=0 y(x/2)=2,y(x/2)=-2/z,其中,n= 命令为:
plot([x,x+h],[y,y1]) else x1=x+h/f; plot([x,x1],[y,y+h]) end end end hold off 读书 P45 中下段。 上机练习： 1 2 0 1 , 2 = − y x y x y dx dy 的斜率场。 2．相平面轨迹表示微分方程的解 关于两个函数 x(t), y(t) 的微分方程组 = = = = ( ) , ( ) . ( ) ( , ), ( ) ( , ), 0 0 0 0 2 1 x t x y t y y t g x y x t g x y 如何直观地表示 x 与 y 的关系呢？ 法一：相平面轨迹图 求出数值解，得到一系列的离散解 , , , ...... 2 2 1 1 0 0 y x y x y x 再在 X0Y 面上连接这些点。如书 P48 图 4.3 法二：利用斜率场作图 ( , ) ( , ) ( , ) 1 2 f x y g x y g x y dx dy 记 = == ，做斜率场即可。如书 P48 图 4.4 §4—4 Matlab 软件求解 1．解析解 （1）一阶微分方程 求 2 1 y dx dy = + 的通解：dsolve('Dy=1+y^2','x') 求 x y dx dy = + − 2 1 的通解：dsolve('Dy=1+x^2-y','x') 求 = = + (0) 1 1 2 y y dx dy 的特解：dsolve('Dy=1+y^2',’y(0)=1’,'x') （2）高阶微分方程 求解 = = − + + − = ( 2) 2, ( 2) 2 . ( ) 0, 2 2 2 y y x y xy x n y 其中， 2 1 n = ，命令为：
dsolve('x2D2y+x *Dy+(x2-052)*y=0, y(pi/2)=2, y( pi2)=2/pi,'x') 求y"-2y+y-4x=0的通解,命令为 dsolve( D3y-2 Dy+y-4*x=0, x') 输出为 ans=8+4*x+Cl°exp(x)+C2*exp(-12*(5^(12)+1)*x)+C3 *exp(12*(5^(12-1)*x) (3)一阶微分方程组 求 f(x)=3f(x)+4g(x) g(x)=-4f(x)+3g(以的通解 If, gl=dsolve(Df=3*f+4*g, Dg=-4*f+3*g, x') 输出为 exp(3* )(cos(4*)* C1+sin(4* x)C g=-exp(3*x)*(sin(4*x)*Cl-cos(4 x)*C2) 若再加上初始条件f(0)=0,g(0)=1,求特解: If, gldsolve(Df=3*f+4*g, Dg=-4*f+3*g, f(0F0, g(0Fl 输出为:f=exp(3*x)*sin(4*x) g=exp(3*x)°coS(4*x) 2.数值解 (1)一阶微分方程 ,0≤x≤1, 求 的数值解,现以步长h=0.1 用“4阶龙格—库塔公式”求解: 先建立“函数M一文件”: function f=eqsl(xy) 再命令:格式为: 自变量因变量]0de45(函数文件名’,节点数组,初始值) 命令为:[x,y=ode45(eqs,0:0.1:1,1) 输出为 X=00.10000.20000.30000.40000.5000 0.60000.70000.80000.900010000 y=1109541.18321.26491.341614142 1.48321.549216125167331.7321 要画图,继续命令:plot(xy) (2)一阶微分方程组 y=cosx+2y-y2,0≤x≤1, y2=sin x-y +2y2, y(O)=0.2,y2(0)=0.3 只须向量化,即可用前面方法。 function f=eqs 2(x, y) f=[cos(x+2*y(l)y(2)sin(x)y(1)+2*y(2)]
dsolve('x^2*D2y+x*Dy+(x^2-0.5^2)*y=0','y(pi/2)=2,Dy( pi/2)=-2/pi','x') 求 y − 2y + y − 4x = 0 的通解，命令为： dsolve('D3y-2*Dy+y-4*x=0','x') 输出为： ans=8+4*x+C1*exp(x)+C2*exp(-1/2*(5^(1/2)+1)*x)+C3 *exp(1/2*(5^(1/2)-1)*x) （3）一阶微分方程组 求 = − + = + ( ) 4 ( ) 3 ( ). ( ) 3 ( ) 4 ( ), g x f x g x f x f x g x 的通解： [f,g]=dsolve('Df=3*f+4*g','Dg=-4*f+3*g','x') 输出为： f =exp(3*x)*(cos(4*x)*C1+sin(4*x)*C2) g =-exp(3*x)*(sin(4*x)*C1-cos(4*x)*C2) 若再加上初始条件 f (0) = 0, g(0) =1 ，求特解： [f,g]=dsolve('Df=3*f+4*g','Dg=-4*f+3*g','f(0)=0,g(0)=1' ,'x') 输出为： f =exp(3*x)*sin(4*x) g =exp(3*x)*cos(4*x) 2．数值解 （1）一阶微分方程 求 = = − (0) 1. , 0 1, 2 y x y x y dx dy 的数值解，现以步长 h=0.1 用“4 阶龙格—库塔公式”求解： 先建立“函数 M—文件”：function f=eqs1(x,y) f=y-2*x/y; 再命令： 格式为： [自变量,因变量]=ode45(‘函数文件名’,节点数组,初始值) 命令为： [x,y]=ode45('eqs1',0:0.1:1,1) 输出为： x =0 0.1000 0.2000 0.3000 0.4000 0.5000 0.6000 0.7000 0.8000 0.9000 1.0000 y =1 1.0954 1.1832 1.2649 1.3416 1.4142 1.4832 1.5492 1.6125 1.6733 1.7321 要画图，继续命令： plot(x,y) （2）一阶微分方程组 = = = − + = + − (0) 0.2, (0) 0.3. sin 2 , cos 2 , 0 1, 1 2 2 1 2 1 1 2 y y y x y y y x y y x 只须向量化，即可用前面方法。 function f=eqs2(x,y) f=[cos(x)+2*y(1)-y(2);sin(x)-y(1)+2*y(2)];
将此函数文件,以文件名eqs2保存后,再下命令: y]ode45(eqs20.0.1:10.2,0.3]) 输出为 X=00.10000.20000.30000.40000.5000 0.60000.70000.80000.90001.0000 0.20000.3000 0.31930.3434 0.45890.3932 0.621304482 0.81010.5066 1.03000.5655 0.6210 158860.6672 194580.6958 2.37240.6954 2.88710.6500 要画图,继续命令: hold on, plot(x, y( 1), plot(x, y( 2), hold off (3)高阶微分方程 先化成一阶微分方程组,再用前面方法 上机练习:{=2)+)-4=05x51 准备:令y1=y,y2=y,y3=y”",化成 y(0) y2(0) y+2y2+4x八(y3(0)(0 用机器:函数文件eqs3内容? 命令? 画图? §4-4微分方程的应用 逻辑斯谛人口模型 基本假设:地球资源限制,人口数有最大值M,人口增 长率k不再是常数,而是随人口数量p的增大而线性递减, 且当p达到最大值M时k=0,记k=r(M-p) a=(M-P)P,P()=P, 得p() P0+(M-P0)e M=,其中F,M为参数 二.弱肉强食的双种群模型 个孤立的海岛上:青草;兔子;狐狸 青草足够多,兔子吃草大量繁殖,兔子多则狐狸容易捕 食,狐狸数量增加,导致兔子减少,狐狸也减少。这样,兔 子、狐狸数量交替增减,无休止的循环,形成生态平衡。 数学家 Volterra(意)建立了这样的“双种群生态模型”: 基本假设:(1)兔子的繁殖率为常数:(2)任何一只兔子与
将此函数文件，以文件名 eqs2 保存后，再下命令： [x,y]=ode45('eqs2',0:0.1:1,[0.2;0.3]) 输出为： x =0 0.1000 0.2000 0.3000 0.4000 0.5000 0.6000 0.7000 0.8000 0.9000 1.0000 y = 0.2000 0.3000 0.3193 0.3434 0.4589 0.3932 0.6213 0.4482 0.8101 0.5066 1.0300 0.5655 1.2868 0.6210 1.5886 0.6672 1.9458 0.6958 2.3724 0.6954 2.8871 0.6500 要画图，继续命令： hold on,plot(x,y(:,1)),plot(x,y(:,2)),hold off （3）高阶微分方程 先化成一阶微分方程组，再用前面方法。 上机练习： = = = − + − = (0) 1, (0) 0.2, (0) 0. 2 4 0, 0 1, y y y y y y x x 准备：令 y = y y = y y = y 1 2 3 , , ，化成 . 0 0.2 1 (0) (0) (0) , 2 4 3 2 1 1 2 3 2 3 2 1 = − + + = y y y y y x y y y y y 用机器： 函数文件 eqs3 内容？ 命令？ 画图？ §4—4 微分方程的应用 一．逻辑斯谛人口模型 基本假设：地球资源限制，人口数有最大值 M，人口增 长率 k 不再是常数，而是随人口数量 p 的增大而线性递减， 且当 p 达到最大值 M 时 k=0，记 k=r(M-p) . 0 0 ( ) , ( ) ( ) r M p p p t p dt dp t = − = ， 得 ( ) 0 0 0 0 ( ) ( ) rM t t p M p e Mp p t − − + − = ，其中 r,M 为参数。 二．弱肉强食的双种群模型 一个孤立的海岛上：青草；兔子；狐狸。 青草足够多，兔子吃草大量繁殖，兔子多则狐狸容易捕 食，狐狸数量增加，导致兔子减少，狐狸也减少。这样，兔 子、狐狸数量交替增减，无休止的循环，形成生态平衡。 数学家 Volterra（意）建立了这样的“双种群生态模型”： 基本假设：（1）兔子的繁殖率为常数；（2）任何一只兔子与
任何一只狐狸相遇的概率是常数,相遇时兔子被杀的概率是 常数:(3)狐狸的死亡率为常数:(4)狐狸的繁殖率与“兔 子数乘狐狸数”成正比 建立模型:t时刻兔子数x(),狐狸数y(),据基本假设得 dx(o) d=x(1)-a·x()y(t), dy(o) =-s:x(1)+b·x()y() dt 用 Volterra模型做数值实验: 1.利用数值解作相平面轨迹图 参数赋值:r=2;s=0.8;a=0.02,b=0.0002; 给初值:兔子x(0)=4000 狐狸数y(O)分别取10、20、30、…、80 先求出数值解,再画(x,y)的相平面轨迹图。 函数文件 function f-fishp53(t, x) r=2s=0.8:a=0.02;b=0.0002 f=[r*x(1)a*x(1)*x(2)-s*x(2)+b*x(1)x(2) 以文件名fshp53保存后,再执行下面程序即可得“图46” clear hold for k1: 8 [t,x]ode45( fish5300.028,x0) plot(x(:,1)x(:2) hold off 2.利用斜率场作图 模型中,两个方程相除得 sy t bxy ()式 代入数据r=2;s=0.8;a=0.02;b=0.0002;在区域 ∫0<x<1200 0<y<200 内作斜率场图,可得“图47”(程序略) 3.利用等值线作图 分离变量法,可得()式的通解为 In hnx-bx=C,其中C为任意常数 初始条件x(0),y(0)决定C的值。每给定一个初值就可在XOY 面上画出一条等值线。现用前面的参数值作等值线图: [xy]= meshgrid(100:10:12000,10:1:200
任何一只狐狸相遇的概率是常数，相遇时兔子被杀的概率是 常数；（3）狐狸的死亡率为常数；（4）狐狸的繁殖率与“兔 子数乘狐狸数”成正比。 建立模型：t 时刻兔子数 x(t) ，狐狸数 y(t) ，据基本假设得 = − + = − ( ) ( ) ( ). ( ) ( ) ( ) ( ), ( ) s x t b x t y t dt dy t r x t a x t y t dt dx t 用 Volterra 模型做数值实验： 1．利用数值解作相平面轨迹图 参数赋值： r = 2;s = 0.8;a = 0.02;b = 0.0002; 给初值：兔子 x(0) = 4000 ； 狐狸数 y(0) 分别取 10、20、30、… 、80 先求出数值解，再画（x，y）的相平面轨迹图。 函数文件 function f=fishp53(t,x) r=2;s=0.8;a=0.02;b=0.0002; f=[r*x(1)-a*x(1)*x(2);-s*x(2)+b*x(1)*x(2)]; 以文件名 fishp53 保存后，再执行下面程序即可得“图 4.6” clear hold on for k=1:8 x0=[4000,10*k]; [t,x]=ode45('fishp53',0:0.02:8,x0); plot(x(:,1),x(:,2)) end hold off 2．利用斜率场作图 模型中，两个方程相除得 rx axy sy bxy dx dy − − + = ， (*)式 代 入 数 据 r = 2;s = 0.8;a = 0.02;b = 0.0002; 在区域 0 200 0 12000 y x 内作斜率场图，可得“图 4.7”（程序 略） 3．利用等值线作图 用分离变量法，可得(*)式的通解为 r ln y − ay + sln x − bx = C ，其中 C 为任意常数。 初始条件 x(0), y(0) 决定C的值。每给定一个初值就可在XOY 面上画出一条等值线。现用前面的参数值作等值线图： clear [x,y]=meshgrid(100:10:12000,10:1:200);
r2s=0.8:a=0.02b=0.0002 logy )-a 说明:20是要画的等值线的条数。 此图中,区域2000<x≤800050y<150的等值线不明,故 缩小范围,再画: [xy}= meshgrid(2000:10:8500,50:1:150) r=2;s=0.8:a=0.02b=0.0002; Zr"log(y)-a*y+s*log(x)-bX contour(x, y, Z, 10) 30004000500060D0 三.计算机网络可靠性分析问题 (看书P56--157) 一个网络系统的三种状态 (1)无故障工作;(2)带故障工作;(3)不工作。 假设“状态转移参数”都与时间t无关,分别为: 1/2(从(1)到(2) y(从(2)到(1) +4/2(从(1)到(3):A2+41(从(2)到(3) 在任意t时刻,系统分别以概率P(1),P(m),P(m)处 于这三种状态。显然有P()+P()+P(1)=1 初始条件:[P(O),P(0)P(O=[100
r=2;s=0.8;a=0.02;b=0.0002; z=r*log(y)-a*y+s*log(x)-b*x; contour(x,y,z,20) 说明：20 是要画的等值线的条数。 此图中，区域 2000<x<8000,50<y<150 的等值线不明，故 缩小范围，再画： clear [x,y]=meshgrid(2000:10:8500,50:1:150); r=2;s=0.8;a=0.02;b=0.0002; z=r*log(y)-a*y+s*log(x)-b*x; contour(x,y,z,10) 三．计算机网络可靠性分析问题 (看书 P56----P57) 一个网络系统的三种状态： （1）无故障工作； （2）带故障工作； （3）不工作。 假设“状态转移参数”都与时间 t 无关，分别为： t 2 （从（1）到（2））； （从（2）到（1））； p + t 2 （从（1）到（3））； p + t （从（2）到（3））。 在任意 t 时刻，系统分别以概率 ( ), ( ), ( ) 1 2 3 P t P t P t 处 于这三种状态。显然有 ( ) ( ) ( ) 1. P1 t + P2 t + P3 t 初始条件：[ (0), (0), (0)] [1,0,0]. P1 P2 P3 =
根据“马尔可夫状态转移原理”建立模型(略写)。 模型求解的主要目标是“系统正常工作的可靠度 最直观的回答就是“画出可靠度曲线”。 函数文件cqsp57如下 function f-eqsp57(t, p) lp=0.0000lt=0.000l;gm=0.01 f=[(lp+lt)*p(1)+gm*p(2)t2*p(1)-(gm+lp+t)*p(2).(lp+ It2)*p(1)+(lp+t)*p(2) 再执行下面程序(文件名sp57)即可 clear t,p}=ode45(eqsp57,0:50:100001001; kkd=1-p(3) plot(t, kkd), grid
根据“马尔可夫状态转移原理”建立模型（略写）。 模型求解的主要目标是“系统正常工作的可靠度”， 最直观的回答就是“画出可靠度曲线”。 函数文件 eqsp57 如下： function f=eqsp57(t,p) lp=0.00001;lt=0.0001;gm=0.01; f=[-(lp+lt)*p(1)+gm*p(2);lt/2*p(1)-(gm+lp+lt)*p(2);(lp+ lt/2)*p(1)+(lp+lt)*p(2)]; 再执行下面程序（文件名 syp57）即可。 clear [t,p]=ode45('eqsp57',0:50:10000,[1;0;0]); kkd=1-p(:,3); plot(t,kkd),grid