第九讲常微分方程模型

第九讲常微分方程模型第1页，共75页，2023年，2月20日，星期一模型分析人口将按指数规律无限增长！人口将始终保持不变！人口将按指数规律减少直至绝灭！人口倍增时间模型求解第2页，共75页，2023年，2月20日，星期一Malthus模型预测美国人口第3页，共75页，2023年，2月20日，星期一Malthus模型预测美国人口第4页，共75页，2023年，2月20日，星期一Malthus模型预测的优缺点优点短期预报比较准确缺点不适合中长期预报原因预报时假设人口增长率r

为常数。没有考虑环境对人口增长的制约作用。第5页，共75页，2023年，2月20日，星期一2.阻滞增长模型假设人口增长率r(t)是t时刻人口x(t)的减函数：其中，xm

为考虑到受自然资源和环境条件限制所能容纳的最大人口数量（称最大人口容量）

模型假设模型建立第6页，共75页，2023年，2月20日，星期一模型分析（定性分析）人口将递减并趋向于xm！人口将始终保持xm不变！人口将递增并趋向于xm！

无论在哪种情况下，人口最终将趋向于最大人口容量！模型求解第7页，共75页，2023年，2月20日，星期一人口增长率达到最大值第8页，共75页，2023年，2月20日，星期一阻滞增长模型预测美国人口第9页，共75页，2023年，2月20日，星期一阻滞增长模型预测美国人口第10页，共75页，2023年，2月20日，星期一阻滞增长模型预测的优缺点优点中期预报比较准确缺点理论上很好，实用性不强原因预报时假设固有人口增长率r以及最大人口容量xm为定值。实际上这两个参数（特别是xm

）很难确定，而且会随着社会发展情况变化而变化。前面图中曲线末端分叉就是由于这个原因。第11页，共75页，2023年，2月20日，星期一6.2药物在体内的分布与排除药物进入机体形成血药浓度(单位体积血液的药物量)血药浓度需保持在一定范围内——给药方案设计

药物在体内吸收、分布和排除过程——药物动力学

建立房室模型——药物动力学的基本步骤房室——机体的一部分，药物在一个房室内均匀分布(血药浓度为常数)，在房室间按一定规律转移本节讨论二室模型——中心室(心、肺、肾等)和周边室(四肢、肌肉等)第12页，共75页，2023年，2月20日，星期一中心室周边室给药排出模型假设中心室(1)和周边室(2),容积不变药物在房室间转移速率及向体外排除

速率，与该室血药浓度成正比药物从体外进入中心室，在二室间

相互转移,从中心室排出体外模型建立第13页，共75页，2023年，2月20日，星期一线性常系数非齐次方程对应齐次方程通解模型建立第14页，共75页，2023年，2月20日，星期一几种常见的给药方式1.快速静脉注射t=0瞬时注射剂量d的药物进入中心室,血药浓度立即为d/V1给药速率f

(t)和初始条件第15页，共75页，2023年，2月20日，星期一2.恒速静脉滴注t>T时,c1(t)和c2(t)按指数规律衰减趋于零药物以恒定速率k进入中心室0Tt££第16页，共75页，2023年，2月20日，星期一吸收室中心室3.口服或肌肉注射相当于药物(剂量d)先进入吸收室，吸收后再进入中心室吸收室药量x0(t)第17页，共75页，2023年，2月20日，星期一参数估计各种给药方式下的c1(t),c2(t)取决于参数k12,k21,k13,V1,V2以快速静脉注射为例,在ti(i=1,2,n)测得c1(ti)由较大的用最小二乘法定A,由较小的用最小二乘法定B,第18页，共75页，2023年，2月20日，星期一参数估计法一进入中心室的药物全部排除第19页，共75页，2023年，2月20日，星期一参数估计法二%构造非线性拟合函数[TWOEXPS.M]functionE=twoexps(a,x,y)x=x(:);y=y(:);Y=a(1)*exp(-a(3)*x)+a(2)*exp(-a(4)*x);E=sum((y-Y).^2)第20页，共75页，2023年，2月20日，星期一a0=[10011];options=optimset('fminsearch');options.TolX=0.01;options.Display='off';a=fminsearch(@ps,a0,options,x,y)

a=[112.23780.18232.1773]第21页，共75页，2023年，2月20日，星期一6.3传染病模型问题描述传染病的传播过程分析受感染人数的变化规律预报传染病高潮到来的时刻预防传染病蔓延的手段按照传播过程的一般规律，用机理分析方法建立模型第22页，共75页，2023年，2月20日，星期一已感染人数(病人)i(t)每个病人每天有效接触(足以使人致病)人数为模型1假设若有效接触的是病人，则不能使病人数增加必须区分已感染者(病人)和未感染者(健康人)建模?第23页，共75页，2023年，2月20日，星期一模型2区分已感染者(病人)和未感染者(健康人)假设1）总人数N不变，病人和健康人的比例分别为2）每个病人每天有效接触人数为,且使接触的健康人致病建模~日接触率SI模型第24页，共75页，2023年，2月20日，星期一模型21/2tmii010ttm~传染病高潮到来时刻(日接触率)tmLogistic模型病人可以治愈！?t=tm,di/dt最大第25页，共75页，2023年，2月20日，星期一模型3传染病无免疫性——病人治愈成为健康人，健康人可再次被感染增加假设SIS模型3）病人每天治愈的比例为~日治愈率建模~日接触率1/~感染期~一个感染期内每个病人的有效接触人数，称为接触数。第26页，共75页，2023年，2月20日，星期一模型3i0i0接触数=1~阈值感染期内有效接触感染的健康者人数不超过病人数1-1/i0模型2(SI模型)如何看作模型3(SIS模型)的特例idi/dt01>10ti>11-1/i0t1di/dt<0第27页，共75页，2023年，2月20日，星期一模型4传染病有免疫性——病人治愈后即移出感染系统，称移出者SIR模型假设1）总人数N不变，病人、健康人和移出者的比例分别为2）病人的日接触率,日治愈率,

接触数=/建模需建立的两个方程第28页，共75页，2023年，2月20日，星期一模型4SIR模型无法求出的解析解在相平面上研究解的性质第29页，共75页，2023年，2月20日，星期一模型4消去dtSIR模型相轨线的定义域相轨线11si0D在D内作相轨线的图形，进行分析第30页，共75页，2023年，2月20日，星期一si101D模型4SIR模型相轨线及其分析传染病蔓延传染病不蔓延s(t)单调减相轨线的方向P1s0imP1:s0>1/

i(t)先升后降至0P2:s0<1/

i(t)单调降至01/

~阈值P3P4P2S0第31页，共75页，2023年，2月20日，星期一模型4SIR模型预防传染病蔓延的手段(日接触率)卫生水平(日治愈率)医疗水平传染病不蔓延的条件——s0<1/的估计降低s0提高r0提高阈值1/降低(=/),群体免疫第32页，共75页，2023年，2月20日，星期一SIR模型被传染人数的估计法一记被传染人数比例x<<s0i0P1i00,s01小,s01提高阈值

1/

降低被传染人数比例xs0-1/=第33页，共75页，2023年，2月20日，星期一被传染人数的估计法二X=fzero(‘x-1.2*log(x/0.96)-0.99’,0.5)

X=0.8651第34页，共75页，2023年，2月20日，星期一6.4多种群生态数学模型意大利生物学家D’Ancona曾致力于鱼类种群相互制约关系的研究，他从第一次世界大战期间,地中海各港口捕获的几种鱼类捕获量百分比的资料中，发现鲨鱼等的比例有明显增加（见下表），而供其捕食的食用鱼的百分比却明显下降.显然战争使捕鱼量下降，食用鱼增加，鲨鱼等也随之增加，但为何鲨鱼的比例大幅增加呢？他无法解释这个现象，于是求助于其岳父，著名的意大利数学家V.Volterra，希望建立一个食饵—捕食系统的数学模型，定性或定量地回答这个问题.第35页，共75页，2023年，2月20日，星期一

该模型反映了在没有人工捕获的自然环境中食饵与捕食者之间的制约关系，没有考虑食饵和捕食者自身的阻滞作用，是Volterra提出的最简单的模型.模型（一）不考虑捕获第36页，共75页，2023年，2月20日，星期一定理

Volterra微分方程组对应初值问题的解是周期函数，且解的周期平均值为

第37页，共75页，2023年，2月20日，星期一首先，建立m-文件shier.m如下：functiondx=shier(t,x)dx=zeros(2,1);dx(1)=x(1)*(1-0.1*x(2));dx(2)=x(2)*(-0.5+0.02*x(1));其次，建立主程序shark.m如下：[t,x]=ode45('shier',[015],[252]);plot(t,x(:,1),'-',t,x(:,2),'*')plot(x(:,1),x(:,2))ToMatlab(shark)第38页，共75页，2023年，2月20日，星期一求解结果：左图反映了x1（t）与x2（t）的关系。可以猜测：x1（t）与x2（t）都是周期函数。第39页，共75页，2023年，2月20日，星期一模型（二）考虑人工捕获设表示捕获能力的系数为e，相当于食饵的自然增长率由a降为a-e，捕食者的死亡率由c增为c+e设战前捕获能力系数e=0.3,战争中降为e=0.1,则战前与战争中的模型分别为:Volterra原理第40页，共75页，2023年，2月20日，星期一模型求解:1、分别用m-文件shier1.m和shier2.m定义上述两个方程2、建立主程序shark1.m,求解两个方程，并画出两种情况下鲨鱼数在鱼类总数中所占比例x2(t)/[x1(t)+x2(t)]ToMatlab(shark1)实线为战前的鲨鱼比例，“*”线为战争中的鲨鱼比例结论：战争中鲨鱼的比例比战前高！第41页，共75页，2023年，2月20日，星期一functiony=shier(t,x)r=1;d=0.5;a=0.1;b=0.02;y=diag([r-a*x(2),-d+b*x(1)])*x;shier.mts=0::0.1:35;x0=[25,2];[t,x]=ode45('shier',ts,x0);[t,x],plot(t,x),grid,gtext('x1(t)'),gtext('x2(t)'),pause,plot(x(:,1),x(:,2)),grid,xlabel('x1'),ylabel('x2')shiyan42注：ts中终值(=15)和步长=(0.1)的确定第42页，共75页，2023年，2月20日，星期一第43页，共75页，2023年，2月20日，星期一计算结果（数值，图形）x(t),y(t)是周期函数，相图(x,y)是封闭曲线；观察，猜测x(t),y(t)的周期约为10.7;xmax=99.3,xmin=2.0,ymax=28.4,ymin=2.0.用数值积分可算出x(t)一周期的平均值为25,y(t)一周期的平均值为10.第44页，共75页，2023年，2月20日，星期一6.5其它生态数学模型存在一大类生态模型源于对Volterra模型的改造模型1考虑食饵种群与外界有迁入或迁出

G.R.Walsh（1978）外界有食饵迁入外界有食饵迁出也可以表示人工干预，如投放或捕获模型讨论食饵－捕食者模型第45页，共75页，2023年，2月20日，星期一模型2考虑食饵种群内部存在生存竞争G.Bojadziev表示当没有捕食者存在时食饵种群的环境容纳量模型3考虑食饵和捕食者种群内部都存在生存竞争张锦炎（1979）表示当没有捕食者存在时食饵种群的环境容纳量第46页，共75页，2023年，2月20日，星期一模型4考虑双方内部都存在生存竞争，且捕食者另有食物来源E.C.Pielou平衡点：第47页，共75页，2023年，2月20日，星期一模型两种群竞争模型竞争模型竞争排斥原理（CompetitionExclutionLaw）

多个种群依靠同一个生存资源而生活，如果生活在同一个地理空间，

猎取相同食物或营养物。在有限的相同生存资源条件下，如果存在竞争关系，

它们必然相互排斥，展开激烈的生存竞争。结局是竞争力较弱的种群灭绝，

竞争力最强的种群达到其环境容纳量。

第48页，共75页，2023年，2月20日，星期一模型讨论种群y最终将被灭绝,种群x最终趋于最大容量种群x最终将被灭绝,种群y最终趋于最大容量存在过正平衡点的一条分界线，将第一象限分成

种群x和种群y的两个吸引域。种群x和种群y最终达到稳定的正平衡态竞争排斥原理

是针对前三种情形得出的结论，

第四种情况极为罕见。第49页，共75页，2023年，2月20日，星期一%M函数functiondy=cwf1(t,y)dy=zeros(2,1);dy(1)=0.1*y(1)*(1-0.001*y(1)-0.0008*y(2));dy(2)=0.2*y(2)*(1-0.0012*y(1)-0.001*y(2));第50页，共75页，2023年，2月20日，星期一

主程序[T,X]=ode45(‘cwf1’,[0200],[200200]);[T,Y]=ode45(‘cwf1’,[0200],[500200]);[T,Z]=ode45(‘cwf1’,[0200],[1200500]);Plot(X(:,1),X(:,2),(Y(:,1),Y(:,2),(Z(:,1),Z(:,2))第51页，共75页，2023年，2月20日，星期一

第52页，共75页，2023年，2月20日，星期一模型两种群互惠模型互惠模型研究多个种群之间相互依赖、共生现象。模型讨论模型有三个平衡点，分别为第53页，共75页，2023年，2月20日，星期一两种群最终达到稳定平衡态两种群共生P3为正平衡点P3稳定两种群最终达不到稳定平衡态P3不稳定第54页，共75页，2023年，2月20日，星期一%m程序functiondy=cwf2(t,y)

dy=zeros(2,1);dy(1)=0.1*y(1)*(1-0.001*y(1)+0.0005*y(2));dy(2)=0.2*y(2)*(-1+0.0015*y(1)-0.001*y(2));第55页，共75页，2023年，2月20日，星期一[T,X]=ode45(‘cwf2’,[0,100],[1600,2800]);[T,Y]=ode45(‘cwf2’,[0,100],[1200,2500]);[T,Z]=ode45(‘cwf2’,[0,100],[2500,2200]);Plot(X(:,1),X(:,2),Y(:,1),Y(:,2),Z(:,1),Z(:,2))Text(2000,2600,2000,2600,'{\sigma_1<1,\sigma_2>1,\sigma_1\sigma_2<1}')第56页，共75页，2023年，2月20日，星期一第57页，共75页，2023年，2月20日，星期一6.6常微分方程的数值解及实验（一）常微分方程数值解的定义在生产和科研中所处理的微分方程往往很复杂且大多得不出一般解。而在实际上对初值问题，一般是要求得到解在若干个点上满足规定精确度的近似值，或者得到一个满足精确度要求的便于计算的表达式。因此，研究常微分方程的数值解法是十分必要的。第58页，共75页，2023年，2月20日，星期一欧拉公式1、用差商代替导数若步长h较小，则有故有公式：此即欧拉法（向前欧拉法）。第59页，共75页，2023年，2月20日，星期一2、使用数值积分对方程y’=f(x,y),两边由xi到xi+1积分，并利用梯形公式，有：实际应用时，与欧拉公式结合使用：此即改进的欧拉法。故有公式：第60页，共75页，2023年，2月20日，星期一例1求解初值问题第61页，共75页，2023年，2月20日，星期一第62页，共75页，2023年，2月20日，星期一龙格—库特方法

考虑微分中值定理

第63页，共75页，2023年，2月20日，星期一2阶龙格－库特公式

第64页，共75页，2023年，2月20日，星期一利用泰勒展式得到第65页，共75页，2023年，2月20日，星期一4阶龙格－库特公式

第66页，共75页，2023年，2月20日，星期一

第67页，共75页，2023年，2月20日，星期一一阶方程组与高阶方程第68页，共75页，2023年，2月20日，星期一

第69页，共75页，2023年，2月20日，星期一第70页，共75页，2023年，2月20日，星期一

Volterra方程解曲线与Volterra方程相轨线第71页，共75页，2023年，2月20日，星期一3、使用泰勒公式以此方法为基础，有