常微分方程和人口模型实验.docx
- 文档编号:12664841
- 上传时间:2023-06-07
- 格式:DOCX
- 页数:20
- 大小:59.40KB
常微分方程和人口模型实验.docx
《常微分方程和人口模型实验.docx》由会员分享,可在线阅读,更多相关《常微分方程和人口模型实验.docx(20页珍藏版)》请在冰点文库上搜索。
常微分方程和人口模型实验
常微分方程和人口模型实验
实验8常微分方程和人口模型实验目的:
(1)了解常微分方程的基本概念。
(2)了解常微分方程的解析解。
(3)了解常微分方程的数值解。
⑷学习、掌握MATLAB^件有关的命令。
实验内容:
常微分方程模型的建立及求解。
一、常微分方程的解析解
调用格式:
dsolve('S','s1','s2',...,'x')
其中,S为方程或方程组,方程S中用D表示求导数,D2,D3,…表示二阶、三阶等高阶导数,方程间用逗号分隔;s1,s2,…为初始条件,x为自变量,初始条件和自变量都可以省略;初始条件缺省时给出带任意常数C1,C2,…的通解;自变量缺省时设定为t(表8.1给出了这个命令的基本用法,请读者运行表中的命令,并观察所得结果。
表8.1dsolve命令用法
常微分方程问题求解命令
2yy'1,,y=dsolve(,Dy=1+yA2?
)
2,yy'1,,y=dsolve(,Dy=1+yA2?
?
y(0)=1?
?
x?
),y(0)1,,
txxxe''2'2,,,x=dsolve(,D2x+2*D1x+2*x=exp(t)?
?
x(0)=1?
?
Dx(0)=0?
)
2,xyyx''3',,y=dsolve(,x*D2y-
3*Dy=xA2?
?
y
(1)=0,y(0)=0?
?
x?
),yy
(1)(0)0,,,
xxy'34,,,[x,y]=dsolve(,Dx=3*x+4*y?
?
Dy=-4*x+3*y?
),yxy'43,,
xxy',,,
yxy',,[x,y]=dsolve(,Dx=x+y,Dy=y-x?
?
Dx(0)=1,Dy(0)=1?
),
xy'(0)'(0)1,,,
二、常微分方程的数值解
在自然科学的众多领域都会遇到常微分方程的求解问题,然而,我们知道,只是有少数十分简单的常微分方程可用简单的方法给出解析解,大部分微分方程只能用数值解法求近似解。
下面用一个例子简要介绍常微分方程数值解法的基本思想。
考虑一阶微分方程初值问题
dy,(,()),,,,ftytttt,0f(8.1)dt,
()yty,00,
t,t,?
t,t所谓数值解法,就是寻求上述初值问题的解在一系列离散节点处的
01nf
t««
yt()yyy,,,h,t,thttn,,()值的的近似值。
称为步长,通常取为常量。
kk,
1kk01nf0最简单的数值解法是欧拉法。
欧拉法的思路很简单:
在节点处用差商近似
代替导数
1ythytyt'()[()()],,kkk,1
将上式代入(8.1)就得到近似计算公式(称为欧拉公式)
t««
yyhftyk,,,(,),0,1,2,kkkk,1
根据初值问题对步长的要求,把常微分方程进一步划分为刚性方程和非刚性方程两类。
所谓刚性方程,就是只有步长充分小时才能得到满意近似解的方程。
MATLA提供了若干直接求解一阶常微分方程(组)初值问题的命令,调用格式如
[x,y]=odeN(,funname?
[t0,tf],y0,tol,trace)
其中,odeN可以是ode23,ode45,ode113,ode15s,ode23s中的任意一个命令;,funname?
是定义
设置解的精函数的文件名;[t0,tf]是常微分方程的求解区间;y0是初始状态列向量;tol用于度,tol可以缺省,缺省时,在ode23中tol=1E-3,在ode45中为tol=1E-6;trace设定是否显示求解的中间结果,缺省值为trace=0,表示不显示中间结果。
注意:
命令ode23,ode45,ode113用于求解非刚性方程,命令ode15s,ode23s用于求解刚性方程;这五个命令采用不同的算法,是用于求解不同的方程,表8.2简要给出了这些命令介绍。
表8.2odeN命令的简要介绍
ode23基于Runge-Kutta(2,3)法,单步法,低阶,精度较低,速度较快。
Ode45基于Runge-Dutta(4,5),单步法,中阶,精度较高,速度较慢。
Ode113采用Adams-Bashforth-moullton法,各(可变)阶方法,多步法。
Ode15s采用多步法,精度较低。
Ode23s采用二阶改进的Rosenbrock法,单步法,速度较快。
对于高阶微分方程,必须先用换元法将其化为化为一阶微分方程组,即“状态方程”,再加以调用Matlab命令求解。
练习1求下列微分方程的数值解,并作解函数的曲线图。
2x,y',y,,0,x,20,y,
y(0),1,
1)建立M函数文件dyex801.m
functiondy=dyex801(x,y)dy=[y-2*x/y];
2)调用程序文件ex801.m
[x,y]=ode45('dyex801',[0,20],1);
Plot(x,y)
xlabel('x')
ylabel('y')
执行结果(图8-1)
>>ex801
为清楚地观察解的变化趋势,图形范围可限制在14?
x?
20,0?
y?
4*10A6,于是
可调整如下
>>plot(x,y,'r-'),axis([142004*10八6])
执行结果如图8-2
图8-1图8-2
练习2求下列微分方程的解,并作解函数的曲线图。
2y?
?
-2(1-y)y?
+y=0,0?
x?
30,y(0)=1,y?
(0)=0
第一步:
先把问题转化为一阶问题。
利用代换法,令y仁y,y2=y?
2所以y?
?
-2(1-y)y?
+y=0等同于下面一阶微分方程组
y1?
=y2
2y2?
=2(1-y1)y2-y1
第二步:
编写M文件。
程序dyex802.m文件
functionz=dyex802(x,y);z=[y
(2);2*(1-y⑴T)*y
(2)-y
(1)];第三步:
编
写执行程序文件ex802.m
clear
y0=[1;0];%设定y的初始值[x,y]=ode23s('dyex902',[0,30],y0);%用
ode23调用ODE文件a6.mplot(x,y(:
1),'r:
',x,y(:
2),'k-');%作解函数y1和
导函数y2的图形xlabel('x')
ylabel('thesolutionofequation')legend('y1','y2')
执行结果如图8-3所示:
图8-3图8-4
练习3求下列微分方程的解,先求解析解,再求数值解,并作解函数的曲线图
进行比y',,y,x,1,y(0),1较。
。
1)求解析解,
>>y=dsolve('Dy=-y+x+1','y(0)=1','x')
y=
x+exp(-x)
tyte,,可得解析解为
2)求数值解,先编写M函数文件dyex803.m
functiondy=dyex803(x,y)
dy=-y+x+1;
3)将解析解和数值解画在同一坐标面上,以比较数值解的误差,M文件
ex803.m[x1,y1]=ode45('dyex803',[0,1],1);%求数值解
x=0:
0.1:
1;
y=x+exp(-x);%解析解
plot(x,y,x1,y1,'ro');%将解析解和数值解画在同一坐标面上legend('解
析解','数值解');
title('微分方程解析解与数值解对比图');
xlabel('x');ylabel('y');
执行结果(图8-4)
>>ex803
从图形上可以看出,此微分方程解析解与数值解是基本吻合的。
三、人口问题
1.问题的提出
长期以来,人类的繁殖一直在自发地进行着(只是由于人口数量的迅速膨胀和环境质量
的急剧恶化,人们才猛然醒悟,开始研究人类和自然的关系、人口数量的变化规律,以及如
何进行人口控制等问题(
表8.3给出1790,1990年两个多世纪间美国的人口统计数据。
请利用这组数据,建立
一个人口增长的数学模型,并用它来它预报2000年、2010年美国的人口(
表8.31790,1990年美国人口统计数据
年(公元)1790180018101820183018401850
人口(百万)3.95.37.29.612.917.123.2
年(公元)1860187018801890190019101920
人口(百万)31.438.650.262.976.092.0106.5
年(公元)1930194019501960197019801990
人口(百万)123.2131.7150.7179.3204.0226.5251.4
2(指数增长模型分析数据中隐藏的规律往往可以帮助我们发现建立模型的线索。
因此,我们利用表8.3
中的数据计算1790,1870年80年间美国的人口十年增长率。
>>x=[3.95.37.29.612.917.123.231.4];%1790-1860年的人口统计数
>>fori=1:
7
r(i)=(x(i+1)-x(i))/x(i);%从1790,1860年的十年增长率
end
>>r
r=0.35900.35850.33330.34380.32560.35670.3534从计算结果中我们发现,在这长达80年的时期内,美国人口的十年增长率的平均值为0.3472,标准差为0.0050。
换言之,在1790,1870年80年间,美国的人口十年增长率基本为一个常数。
因此,我们提出如下假设:
假设1、美国人口的十年增长率为一个常数。
txt()x以1790为0时刻,记时刻0美国的人口数为,时刻美国人口数为。
0
xt()假设2视为连续可微函数。
xt()假设2为利用微积分这一数学工具而做出的。
由于美国人口数的规模很大,的取值都是很大的整数(假设2所带来的误差远小于人口统计工作中不可避免的误差。
因此,这条假设是适当的。
[,]ttt,,考虑时间区间上人口的增量,我们有
xttxtrxttot()()()(),,,,,,,
xt(),t,,t0注意到的可微性,上式两边同除以,并令得
dx,rxdt
加上0时刻的人口数条件,我们得如下一阶常微分方程初值问题
dx(8.2),,rxxx,(0)0dt
在Matlab中,方程
(2)可以通过以下命令求解:
>>x=dsolve('Dx=r*x','x(0)=x0','t')
x=
exp(r*t)*x0
所以,方程
(2)的解为
rt,xtxe(),(8.3)0
(8.3)告诉我们,在人口增长率为常数的假设下,人口将按指数函数的形式增长。
它最早由英国经济学家、统计学家和人口学家马尔萨斯提出,因此也被称之为马尔萨斯模型。
根据1790,1990年美国的人口统计数据,我们可以统计方法估计出用(8.3)中的参数xr,。
将这两个参数值代入(8.3)就可以预测未来某一时刻的人口数。
0
历史上,这个模型与19世纪以前欧洲一些地区人口统计数据可以很好地吻合,迁往加拿大的欧洲移民后代人口也大致符合这个模型(另外,用它作短期人口预测可以得到较好的结果。
但是长期来看,这个模型的预测效果很差,预测结果根本无法采用。
根据表8.3中的数据,我们在MATLAB^计算出美国人口从1790,1990年的每十年增长率:
>>x=[3.95.37.29.612.917.123.231.438.650.262.976.092.0...
106.5123.2131.7150.7179.3204.0226.5251.4];%1790-1990年的人口
统计数
>>fori=1:
20
r(i)=(x(i+1)-x(i))/x(i);%从1790,1990年的每十年增长率
end
>>r
r=0.35900.35850.33330.34380.32560.35670.35340.2293
0.30050.25300.20830.21050.15760.15680.06900.1443
0.18980.13780.11030.1099
从以上结果可以看到,尽管在1790,1870年的70年间,美国人口的每十年增长率的波动不大,基本保持为一个常数(大约0.35)增长率,但是,1870开始,十年增长率就呈现明显的下降趋势,随时间的推移,十年增长率偏离0.35的程度越来越大.这说明,用马尔萨斯模型来预测当代美国人口的变化规律是不恰当的.
3.阻滞增长模型(Logistic模型)
通过前面的分析,我们看到:
人口增长率是准确刻画人口增长规律一个关键的因素,且“人口增长率为一个常数”的假设不合理,那么,我们应当怎样刻画人口增长率的变化规律呢,为此,我们作出1790,1990年美国十年人口增长率随人口数变化而变化的图像。
>>x=[3.95.37.29.612.917.123.231.438.650.262.976.092.0...
106.5123.2131.7150.7179.3204.0226.5251.4];%1790-1990年的人口
统计数
>>fori=1:
20
r(i)=(x(i+1)-x(i))/x(i);%从1790,1990年的每十年增长率
end
>>x1=x(2:
21);
>>plot(x1,r)
图8.51790,1990年美国十年人口增长率关于人口的变化趋势
从图上我们可以清楚地看到,随人口规模的增加,十年人口增长率总体上呈现出线性下降的趋势。
那么,为什么会出现这种趋势呢,著名的草履虫试验告诉我们:
在环境条变化不大的情
x况下,既定环境对生活在其中的生物种群都有一个最大供养量,当种群数量接近这个最m
大供养量时,就会出现因生存资源短缺、生活空间紧张而降低增长率的现象,当种群数量达到最大供养量时,种群数量停止增长。
人类也是一种生物种群,因此,在一个国家或地区内
rx()rx()0,的人口增长率实际上是该地区人口规模的减函数,并且
综合以上分析,我们提出如下新的假设:
rx()假设1?
、美国人口的十年增长率是人口数的线性函数
x(8.4)()
(1)rxr,,,xm
xx,0这里,称为固有增长率,表示人口很少时(理论上是)的增长率(是自然资源和rm
xx,环境条件所能容纳的最大人口数量,称为人口容量,当时人口不再增长,即增长率m
rx()0,。
m
假设1?
考虑到环境对人口增长的阻滞作用,相应地把模型称为阻滞增长模型。
rx()用(8.4)中的代替(8.2)中的就得到如下人口的阻滞增长模型r
dxx(8.5)
(1),(0),,,rxxx0dtxm
在Matlab中,方程(6)可以通过以下命令求解:
>>x=dsolve('Dx=r*x*(1-x/xm)','x(0)=x0','t')
x=
xm/(1+exp(-r*t)*(xm-x0)/x0)
所以,方程(6)的解为
xm(),xt(8.6),,x,rtm,,1,,1e,,x0,,
由方程(8.6)表示的阻滞增长模型,是荷兰生物数学家Verhulst19世纪中叶提出的(它不仅能够大体上符合人口及许多生物数量的变化规律,而且在社会经济领域也有广泛的应用,例如耐用消费品的售量就可以用它来描述(基于这个模型能够描述一些事物符合逻辑的客观规律,人们常称它为逻辑斯蒂(Logistic)模型(
4.模型的参数估计
用指数增长模型或阻滞增长模型进行人口预报,先要作参数估计(除了初始人
口x外,0指数增长模型要估计r,阻滞增长模型则要估计r和x(它们可以用人口
统计数据拟合得到(m
通常使用的数据拟合方法是最小二乘法,这里只给出利用表9-1中1790—1980年的数据对r和x的拟合结果.m
在Matlab中,非线性最小二乘的数据拟合用以下命令:
g=lsqnonlin(,f?
g0)f为拟合参数的函数文件,g0为拟合参数的初值
下面我们通过⑺式,利用表13-1中1790—1990年的数据对r和x进行拟合,首先作变换:
m
k=1790+10t,k=1790,1800,…,1990;t=0,1,…,20
xmx(t),于是x=x(0)=3.9,(7)式变为0x,,,rtme1,,1,,3.9,,
xmx(t),令f=,有两个参数x和r.mx,,,rtme1,,1,,3.9,,
现在,分两步对r和x进行拟合.m
第一步:
定义参数r和x的函数文件(.M文件).文件名为xt.mm
functionf=xt(g)%函数定义,g为参数名,函数名为xt,必须与函数文件名一
致
t=0:
20;%对应于1790-1990年的t值
x=[3.95.37.29.612.917.123.231.438.650.262.976.092.0...
106.5123.2131.7150.7179.3204.0226.5251.4];%1790-1990年的人口
统计数
f=x-g
(1)./(1+(g
(1)./3.9-1).*exp(-g
(2).*t));%g有两个参数,一个是g
(1)
对应于xm,另一个是
g
(2)对应于r
第二步:
编写执行文件ex904.m
g0=[300,0.2];%参数xm,r的初值,也可以取另外的初值,有可能产生不同的拟
合结果g=leastsq('xt',g0)%调用参数函数xt按初值进行拟合g
(1),g
(2)
t=0:
20;%相应于1790-1990年
y=1790+10.*t%转化为相应的年份
f=g
(1)./(1+(g
(1)./3.9-1).*exp(-g
(2).*t))%对拟合的参数与人口统计数据
相比是否吻合
于是在Matlab空间,有
>>ex904
g=
311.95410.2798(即xm=311.9541,r=0.2798)y=
Columns1through10
179018001810182018301840
1850186018701880
Columns11through20
189019001910192019301940
1950196019701980
Column21
1990
f=
Columns1through12
3.90005.13856.76178.882811.644015.220219.822425.6956
33.111742.351753.673767.2672
Columns13through21
83.1946101.3315121.3253142.5932164.3747185.8325206.1781
224.7812
241.2348
5.模型检验
我们可以用模型预报1990年美国的人口,与已知的实际数据比较,来检验模型是否合适.
注意,这样做的前提是,在估计模型参数时,不能再用1990年的数据.
首先,我们在文件xt.m中修改数据如下:
t=0:
19;%对应于1790-1980年的t值
x=[3.95.37.29.612.917.123.231.438.650.262.976.092.0...
106.5123.2131.7150.7179.3204.0226.5];%1790-1980年的人口统计数f=x-g
(1)./(1+(g
(1)./3.9-1).*exp(-g
(2).*t));%g有两个参数,一个是g
(1)对应于xm,另一个是
g
(2)对应于r
在文件ex904.m中修改数据如下:
g0=[300,0.2];%参数xm,r的初值,也可以取另外的初值,有可能产生不同的拟合结果g=lsqnonlin('xt',g0);%调用参数函数xt按初值进行拟合g
(1),g
(2)xm=g
(1)
r=g
(2)
t=0:
20;%相应于1790-1990年
y=1790+10.*t%转化为相应的年份
f=g
(1)./(1+(g
(1)./3.9-1).*exp(-g
(2).*t))%对拟合的参数与人口统计数据
相比是否吻合于是,
xm=
285.8958
r=
0.2858
y=
Columns1through6
179018001810182018301840
Columns7through12
185018601870188018901900
Columns13through18
191019201930194019501960
Columns19through21
197019801990f=
Columns1through7
3.90005.16696.83549.025311.887015.604920.3990
Columns8through14
26.521534.244643.836555.522469.429985.5277103.5721
Columns15through21
123.0849143.3828163.6632183.1261201.0957217.1036230.9158可以看出,以上结果与实际数据相差较大,为探究其原因,我们重新审视图8.5。
从图上可以看出,一方面,从1890年以后,每十年增长率基本上是随人口的增长而减少,
另一方面,每十年增长率基本是人口数的折线函数,拐点大约发生在1亿人口附近,因此,
这一时期美国人口的变化基本符合逻辑斯蒂(Logistic)模型,为保证预测的准
确定,把参数
值的样本区间限定为1900,1980年间的数据,再次估计,所以我们以1900年的人口为初值进
行检验,此时x=x(0)=76.0.0
我们在文件xt.m中修改数据如下:
functionf=xt(g)%函数定义,g为参数名,函数名为xt,必须与函数文件名一
致t=0:
8;%对应于1900-1980年的t值
x=[76.092.0106.5123.2131.7150.7179.3204.0226.5];%1900-1980年的人口统计数f=x-g
(1).
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 微分方程 人口 模型 实验