第四章 种群数量的状态转移----微分方程人口数量，物种数量，网络系统元件的运行状态。
状态转移：变化率（即：某个物理量导数）。
微分方程模型：导数与各个变量之间的某种平衡关系。
本章内容：微分方程（组）的模型建立、数值解、图形解；用Matlab几何直观地展示各种求解方法的求解结果。
§4—1 人口问题
马尔萨斯（Malthus），英 1766--1834 关于人口数量的指数增长模型：
基本假设：人口增长率（即：单位时间内增长的人口数量与总人口数量之比）为常数。
模型建立：时刻的人口数量为已知，记为，人口增长率常数记为，人口数量是时间的函数，，在短时间内，人口增长数量为
，
得 ，
解得 （这就是马尔萨斯人口模型）
模型求解：为了确定模型中的参数，不妨借用美国从1790年至1840年的人口数据：
年份,1790 1800 1810 1820 1830 1840
人口/：3929 5308 7240 9638 12866 17069
代入得
模型检验：用此结果来预测美国从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
预测值：582000 780900 1047700 1405800 1886300
实际值：179323 203212 226505 248710 281416
相对误差：225% 284% 363% 465% 570%
---------------------------------------------------------------------
该模型作长期预测不合理。
§4—2 微分方程的数值解法
略。
§4—3 微分方程图解法
1．斜率场表示微分方程的解给定，则过平面直角坐标系上任意一点，都可以按照斜率做一条短直线。这样的短直线布满整个坐标平面，形成的图形称为的斜率场。
斜率场的构造方法：
（1）明确矩形区域；
（2）沿着横向（x方向）、纵向（y方向）分别以小步长分割该区域；
（3）过每个结点分别做短直线，如过结点做直线：若则代入解得，连接两点即可；否则代入解得，连接两点即可。
例1：画出的斜率场。
解：考虑在区域上，以步长来分割，,
程序清单：
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,y1])
else
x1=x+h/f;
plot([x,x1],[y,y+h])
end
end
end
hold off
读书P45 中下段。
上机练习：的斜率场。
2．相平面轨迹表示微分方程的解
关于两个函数的微分方程组
如何直观地表示与的关系呢？
法一：相平面轨迹图
求出数值解，得到一系列的离散解
再在X0Y面上连接这些点。如书P48 图4.3
法二：利用斜率场作图
，做斜率场即可。如书P48 图4.4
§4—4 Matlab软件求解
1．解析解
（1）一阶微分方程求的通解：dsolve('Dy=1+y^2','x')
求的通解：dsolve('Dy=1+x^2-y','x')
求的特解：dsolve('Dy=1+y^2',’y(0)=1’,'x')
（2）高阶微分方程求解其中，，命令为：
dsolve('x^2*D2y+x*Dy+(x^2-0.5^2)*y=0','y(pi/2)=2,Dy(pi/2)=-2/pi','x')
求的通解，命令为：
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）一阶微分方程组求的通解：
[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,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）一阶微分方程求的数值解，现以步长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）一阶微分方程组
只须向量化，即可用前面方法。
function f=eqs2(x,y)
f=[cos(x)+2*y(1)-y(2);sin(x)-y(1)+2*y(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）高阶微分方程先化成一阶微分方程组，再用前面方法。
上机练习：
准备：令，化成
用机器,函数文件eqs3内容？
命令？
画图？
§4—4 微分方程的应用
一．逻辑斯谛人口模型
基本假设：地球资源限制，人口数有最大值M，人口增长率k不再是常数，而是随人口数量p的增大而线性递减，且当p达到最大值M时k=0，记 k=r(M-p),
，
得 ，其中为参数。
二．弱肉强食的双种群模型
一个孤立的海岛上：青草；兔子；狐狸。
青草足够多，兔子吃草大量繁殖，兔子多则狐狸容易捕食，狐狸数量增加，导致兔子减少，狐狸也减少。这样，兔子、狐狸数量交替增减，无休止的循环，形成生态平衡。
数学家Volterra（意）建立了这样的“双种群生态模型”：
基本假设：（1）兔子的繁殖率为常数；（2）任何一只兔子与任何一只狐狸相遇的概率是常数，相遇时兔子被杀的概率是常数；（3）狐狸的死亡率为常数；（4）狐狸的繁殖率与“兔子数乘狐狸数”成正比。
建立模型：时刻兔子数，狐狸数，据基本假设得
用Volterra模型做数值实验：
1．利用数值解作相平面轨迹图
参数赋值：
给初值：兔子；
狐狸数分别取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．利用斜率场作图模型中，两个方程相除得
,(*)式代入数据在区域内作斜率场图，可得“图4.7”（程序 略）
3．利用等值线作图
用分离变量法，可得(*)式的通解为
，其中C为任意常数。
初始条件决定C的值。每给定一个初值就可在XOY面上画出一条等值线。现用前面的参数值作等值线图：
clear
[x,y]=meshgrid(100:10:12000,10:1:200);
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无关，分别为：
（从（1）到（2））； （从（2）到（1））；
（从（1）到（3））； （从（2）到（3））。
在任意t时刻，系统分别以概率处于这三种状态。显然有
初始条件：
根据“马尔可夫状态转移原理”建立模型（略写）。
模型求解的主要目标是“系统正常工作的可靠度”，最直观的回答就是“画出可靠度曲线”。
函数文件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
状态转移：变化率（即：某个物理量导数）。
微分方程模型：导数与各个变量之间的某种平衡关系。
本章内容：微分方程（组）的模型建立、数值解、图形解；用Matlab几何直观地展示各种求解方法的求解结果。
§4—1 人口问题
马尔萨斯（Malthus），英 1766--1834 关于人口数量的指数增长模型：
基本假设：人口增长率（即：单位时间内增长的人口数量与总人口数量之比）为常数。
模型建立：时刻的人口数量为已知，记为，人口增长率常数记为，人口数量是时间的函数，，在短时间内，人口增长数量为
，
得 ，
解得 （这就是马尔萨斯人口模型）
模型求解：为了确定模型中的参数，不妨借用美国从1790年至1840年的人口数据：
年份,1790 1800 1810 1820 1830 1840
人口/：3929 5308 7240 9638 12866 17069
代入得
模型检验：用此结果来预测美国从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
预测值：582000 780900 1047700 1405800 1886300
实际值：179323 203212 226505 248710 281416
相对误差：225% 284% 363% 465% 570%
---------------------------------------------------------------------
该模型作长期预测不合理。
§4—2 微分方程的数值解法
略。
§4—3 微分方程图解法
1．斜率场表示微分方程的解给定，则过平面直角坐标系上任意一点，都可以按照斜率做一条短直线。这样的短直线布满整个坐标平面，形成的图形称为的斜率场。
斜率场的构造方法：
（1）明确矩形区域；
（2）沿着横向（x方向）、纵向（y方向）分别以小步长分割该区域；
（3）过每个结点分别做短直线，如过结点做直线：若则代入解得，连接两点即可；否则代入解得，连接两点即可。
例1：画出的斜率场。
解：考虑在区域上，以步长来分割，,
程序清单：
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,y1])
else
x1=x+h/f;
plot([x,x1],[y,y+h])
end
end
end
hold off
读书P45 中下段。
上机练习：的斜率场。
2．相平面轨迹表示微分方程的解
关于两个函数的微分方程组
如何直观地表示与的关系呢？
法一：相平面轨迹图
求出数值解，得到一系列的离散解
再在X0Y面上连接这些点。如书P48 图4.3
法二：利用斜率场作图
，做斜率场即可。如书P48 图4.4
§4—4 Matlab软件求解
1．解析解
（1）一阶微分方程求的通解：dsolve('Dy=1+y^2','x')
求的通解：dsolve('Dy=1+x^2-y','x')
求的特解：dsolve('Dy=1+y^2',’y(0)=1’,'x')
（2）高阶微分方程求解其中，，命令为：
dsolve('x^2*D2y+x*Dy+(x^2-0.5^2)*y=0','y(pi/2)=2,Dy(pi/2)=-2/pi','x')
求的通解，命令为：
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）一阶微分方程组求的通解：
[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,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）一阶微分方程求的数值解，现以步长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）一阶微分方程组
只须向量化，即可用前面方法。
function f=eqs2(x,y)
f=[cos(x)+2*y(1)-y(2);sin(x)-y(1)+2*y(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）高阶微分方程先化成一阶微分方程组，再用前面方法。
上机练习：
准备：令，化成
用机器,函数文件eqs3内容？
命令？
画图？
§4—4 微分方程的应用
一．逻辑斯谛人口模型
基本假设：地球资源限制，人口数有最大值M，人口增长率k不再是常数，而是随人口数量p的增大而线性递减，且当p达到最大值M时k=0，记 k=r(M-p),
，
得 ，其中为参数。
二．弱肉强食的双种群模型
一个孤立的海岛上：青草；兔子；狐狸。
青草足够多，兔子吃草大量繁殖，兔子多则狐狸容易捕食，狐狸数量增加，导致兔子减少，狐狸也减少。这样，兔子、狐狸数量交替增减，无休止的循环，形成生态平衡。
数学家Volterra（意）建立了这样的“双种群生态模型”：
基本假设：（1）兔子的繁殖率为常数；（2）任何一只兔子与任何一只狐狸相遇的概率是常数，相遇时兔子被杀的概率是常数；（3）狐狸的死亡率为常数；（4）狐狸的繁殖率与“兔子数乘狐狸数”成正比。
建立模型：时刻兔子数，狐狸数，据基本假设得
用Volterra模型做数值实验：
1．利用数值解作相平面轨迹图
参数赋值：
给初值：兔子；
狐狸数分别取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．利用斜率场作图模型中，两个方程相除得
,(*)式代入数据在区域内作斜率场图，可得“图4.7”（程序 略）
3．利用等值线作图
用分离变量法，可得(*)式的通解为
，其中C为任意常数。
初始条件决定C的值。每给定一个初值就可在XOY面上画出一条等值线。现用前面的参数值作等值线图：
clear
[x,y]=meshgrid(100:10:12000,10:1:200);
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无关，分别为：
（从（1）到（2））； （从（2）到（1））；
（从（1）到（3））； （从（2）到（3））。
在任意t时刻，系统分别以概率处于这三种状态。显然有
初始条件：
根据“马尔可夫状态转移原理”建立模型（略写）。
模型求解的主要目标是“系统正常工作的可靠度”，最直观的回答就是“画出可靠度曲线”。
函数文件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
