《第九讲 常微分方程模型2.ppt》由会员分享，可在线阅读，更多相关《第九讲 常微分方程模型2.ppt（75页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。
1、6.1 人口增长模型人口增长模型英国人口学家Malthus(1766-1834)模型假设模型假设人口自然增人口自然增长长率率 r 为常数为常数即即单单位位时间时间内人口的增内人口的增长长量与当量与当时时的人口呈正比的人口呈正比。模型建立模型建立1.指数增长模型指数增长模型人口以几何级数增加！人口以几何级数增加！模型分析模型分析人口将人口将按指数规律无限增长按指数规律无限增长！人口将人口将始终保持不变始终保持不变！人口将人口将按指数规律减少直至绝灭按指数规律减少直至绝灭！人口倍增时间人口倍增时间模型求解模型求解Malthus模型预测美国人口Malthus模型预测美国人口Malthus模型预测的优
2、缺点优点优点短期预报比较准确缺点缺点不适合中长期预报原因原因预报时假设人口增长率 r 为常数。没有考虑环境对人口增长的制约作用。2.阻滞增长模型阻滞增长模型假设人口增长率假设人口增长率 r(t)是是 t 时刻时刻人口人口 x(t)的减函数的减函数：其中，其中，xm 为为考考虑虑到受自然到受自然资资源和源和环环境条境条件限制所能容件限制所能容纳纳的最大人口数量的最大人口数量（称（称最大人口容量最大人口容量）模型假设模型假设模型建立模型建立模型分析（定性分析）模型分析（定性分析）人口将人口将递减并趋向于递减并趋向于xm！人口将人口将始终保持始终保持xm不变不变！人口将人口将递增并趋向于递增并趋向于
3、xm！无论在哪种情况下，人口最终将趋向于最大人口容量！无论在哪种情况下，人口最终将趋向于最大人口容量！模型求解模型求解人口增长率达到最大值人口增长率达到最大值阻滞增长模型预测美国人口阻滞增长模型预测美国人口阻滞增长模型预测的优缺点优点优点 中期预报比较准确缺点缺点 理论上很好，实用性不强原因原因预报时假设固有人口增长率 r 以及最大人口容量 xm 为定值。实际上这两个参数（特别是 xm）很难确定，而且会随着社会发展情况变化而变化。前面图中曲线末端分叉就是由于这个原因。6.2 药物在体内的分布与排除药物在体内的分布与排除 药物进入机体形成药物进入机体形成血药浓度血药浓度(单位体积血液的药物量单位
4、体积血液的药物量)血药浓度需保持在一定范围内血药浓度需保持在一定范围内给药方案设计给药方案设计 药物在体内吸收、分布和排除过程药物在体内吸收、分布和排除过程 药物动力学药物动力学 建立建立房室模型房室模型药物动力学的基本步骤药物动力学的基本步骤 房室房室机体的一部分，药物在一个房室内均匀机体的一部分，药物在一个房室内均匀分布分布(血药浓度为常数血药浓度为常数)，在房室间按一定规律转移，在房室间按一定规律转移 本节讨论本节讨论二室模型二室模型中心室中心室(心、肺、肾等心、肺、肾等)和和周边室周边室(四肢、肌肉等四肢、肌肉等)中心室中心室周边室周边室给药给药排出排出模型假设模型假设 中心室中心室(
5、1)和周边室和周边室(2),容积不变容积不变 药物在房室间转移速率及向体外排除药物在房室间转移速率及向体外排除 速率，与该室血药浓度成正比速率，与该室血药浓度成正比 药物从体外进入中心室，在二室间药物从体外进入中心室，在二室间 相互转移相互转移,从中心室排出体外从中心室排出体外模型建立模型建立线性常系数线性常系数非齐次方程非齐次方程对应齐次对应齐次方程通解方程通解模型建立模型建立几种常见的给药方式几种常见的给药方式1.快速静脉注射快速静脉注射t=0 瞬时瞬时注射剂量注射剂量d 的的药物进入中心室药物进入中心室,血药血药浓度立即为浓度立即为 d/V1给药速率给药速率 f(t)和初始条件和初始条件
6、2.恒速静脉滴注恒速静脉滴注t T时时,c1(t)和和 c2(t)按按指数规律衰减趋指数规律衰减趋于零于零药物以恒定速率药物以恒定速率k 进入中心室进入中心室0Tt 吸收室中心室3.口服或肌肉注射口服或肌肉注射相当于药物相当于药物(剂量剂量d)先进入吸收室，吸收后再进入中心先进入吸收室，吸收后再进入中心室室吸收室药量吸收室药量x0(t)参数估计参数估计各种给药方式下的各种给药方式下的 c1(t),c2(t)取决于参数取决于参数k12,k21,k13,V1,V2以以快速静脉注射快速静脉注射为例为例 ,在在ti(i=1,2,n)测得测得c1(ti)由由较大的较大的 用用最小二乘法最小二乘法定定A
7、A,由由较小的较小的 用用最小二乘法最小二乘法定定B,参数估计法一参数估计法一进入中心室的药物全部排除进入中心室的药物全部排除参数 估计法二n%构造非线性拟合函数构造非线性拟合函数nTWOEXPS.Mnfunction E=twoexps(a,x,y)nx=x(:);y=y(:);nY=a(1)*exp(-a(3)*x)+a(2)*exp(-a(4)*x);nE=sum(y-Y).2)na0=100 1 1;noptions=optimset(fminsearch);noptions.TolX=0.01;noptions.Display=off;na=fminsearch(ps,a0,opti
8、ons,x,y)a=112.2378 0.1823 2.1773 6.3 传染病模型传染病模型问题问题 描述传染病的传播过程描述传染病的传播过程 分析受感染人数的变化规律分析受感染人数的变化规律 预报传染病高潮到来的时刻预报传染病高潮到来的时刻 预防传染病蔓延的手段预防传染病蔓延的手段 按照传播过程的一般规律，按照传播过程的一般规律，用机理分析方法建立模型用机理分析方法建立模型 已感染人数已感染人数(病人病人)i(t)每个病人每天有效接触每个病人每天有效接触(足以使人致病足以使人致病)人数为人数为 模型模型1 1假设假设若有效接触的是病人，若有效接触的是病人，则不能使病人数增加则不能使病人数增
9、加必须区分已感染者必须区分已感染者(病人病人)和未感染者和未感染者(健康人健康人)建模建模?模型模型2 2区分已感染者区分已感染者(病人病人)和未感染者和未感染者(健康健康人人)假设假设1）总人数）总人数N不变，病人和健康不变，病人和健康 人的人的 比例分别为比例分别为 2）每个病人每天有效接触人数）每个病人每天有效接触人数为为,且且使接触的健康人致病使接触的健康人致病建模建模 日日接触率接触率SI 模型模型模型模型21/2tmii010ttm传染病高潮到来时刻传染病高潮到来时刻 (日接触率日接触率)tm Logistic 模型病人可以治愈！病人可以治愈！?t=tm,di/dt 最最大大模型模
10、型3传染病无免疫性传染病无免疫性病人治愈成病人治愈成为健康人，健康人可再次被感染为健康人，健康人可再次被感染增加假设增加假设SIS 模型模型3）病人每天治愈的比例为）病人每天治愈的比例为 日日治愈率治愈率建模建模 日接触率日接触率1/感染期感染期 一个感染期内一个感染期内每个病人的每个病人的有效接触人数，称为有效接触人数，称为接触数接触数。模型模型3i0i0接触数接触数 =1 阈值阈值感染期内感染期内有效接触感染的有效接触感染的健康者人数不超过病人数健康者人数不超过病人数1-1/i0模型模型2(SI模型模型)如何看作模型如何看作模型3(SIS模型模型)的特例的特例idi/dt01 10ti 1
11、1-1/i0t 1di/dt 1/i(t)先升后降至先升后降至0P2:s01/i(t)单调降至单调降至01/阈值阈值P3P4P2S0模型模型4SIR模型模型预防传染病蔓延的手段预防传染病蔓延的手段 (日接触率日接触率)卫生水平卫生水平 (日日治愈率治愈率)医疗水平医疗水平 传染病不蔓延的条件传染病不蔓延的条件s01/的估计的估计 降低降低 s0提高提高 r0 提高阈值提高阈值 1/降低降低 (=/),群体免疫群体免疫SIR模型模型被传染人数的估计法一被传染人数的估计法一 记被传染人数比例记被传染人数比例xs0i0P1i0 0,s0 1 小小,s0 1提高阈值提高阈值 1/降低降低被传染人数比例
12、被传染人数比例 xs0-1/=被传染人数的估计法二被传染人数的估计法二nX=fzero(x-1.2*log(x/0.96)-0.99,0.5)n X=0.86516.4多种群生态数学模型多种群生态数学模型 意大利生物学家DAncona曾致力于鱼类种群相互制约关系的研究，他从第一次世界大战期间,地中海各港口捕获的几种鱼类捕获量百分比的资料中，发现鲨鱼等的比例有明显增加（见下表），而供其捕食的食用鱼的百分比却明显下降.显然战争使捕鱼量下降，食用鱼增加，鲨鱼等也随之增加，但为何鲨鱼的比例大幅增加呢？他无法解释这个现象，于是求助于其岳父，著名的意大利数学家V.Volterra，希望建立一个食饵捕食系统
13、的数学模型，定性或定量地回答这个问题.该模型反映了在没有人工捕获的自然环境中食饵与捕食者之间的制约关系，没有考虑食饵和捕食者自身的阻滞作用，是Volterra提出的最简单的模型.模型（一）模型（一）不考虑捕获定理定理 Volterra微分方程组对应初值问题的解是周期函数，且解的周期平均值为 首先，建立m-文件shier.m如下：function dx=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,0 15,25 2);plot
14、(t,x(:,1),-,t,x(:,2),*)plot(x(:,1),x(:,2)To Matlab(shark)求解结果：左图反映了x1（t）与x2（t）的关系。可以猜测：x1（t）与x2（t）都是周期函数。模型（二）模型（二）考虑人工捕获 设表示捕获能力的系数为e，相当于食饵的自然增长率由a降为a-e，捕食者的死亡率由c增为 c+e设战前捕获能力系数e=0.3,战争中降为e=0.1,则战前与战争中的模型分别为:Volterra原理原理模型求解:1、分别用m-文件shier1.m和shier2.m定义上述两个方程2、建立主程序shark1.m,求解两个方程，并画出两种情况下鲨鱼数在鱼类总数中
15、所占比例 x2(t)/x1(t)+x2(t)To Matlab(shark1)实线为战前的鲨鱼比例，“*”线为战争中的鲨鱼比例结论：战争中鲨鱼的比例比战前高！结论：战争中鲨鱼的比例比战前高！function y=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)
16、,ylabel(x2)shiyan42注：ts 中终值(=15)和步长=(0.1)的确定计算结果（数值，图形）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.6.5其它生态数学模型其它生态数学模型存在一大类生态模型源于对Volterra模型的改造模型模型1考虑食饵种群与外界有迁入或迁出 G.R.
17、Walsh（1978）外界有食饵迁入外界有食饵迁出也可以表示人工干预，如投放或捕获模型讨论模型讨论模型模型2考虑食饵种群内部存在生存竞争 G.Bojadziev表示当没有捕食者存在时食饵种群的环境容纳量模型模型3考虑食饵和捕食者种群内部都存在生存竞争 张锦炎张锦炎（1979）表示当没有捕食者存在时食饵种群的环境容纳量模型模型4考虑双方内部都存在生存竞争，且捕食者另有食物来源 E.C.Pielou 平衡点平衡点：模型模型两种群竞争模型两种群竞争模型竞争排斥原理竞争排斥原理（Competition Exclution Law）多个种群依靠同一个生存资源而生活，如果生活在同一个地理空间，猎取相同食物
18、或营养物。在有限的相同生存资源条件下，如果存在竞争关系，它们必然相互排斥，展开激烈的生存竞争。结局是竞争力较弱的种群灭绝，竞争力最强的种群达到其环境容纳量。模型讨论模型讨论种群 y 最终将被灭绝,种群 x 最终趋于最大容量 种群 x 最终将被灭绝,种群 y 最终趋于最大容量 存在过正平衡点的一条分界线，将第一象限分成种群 x 和种群 y 的两个吸引域。种群 x 和种群 y 最终达到稳定的正平衡态竞争排斥原理竞争排斥原理 是针对前三种情形得出的结论，第四种情况极为罕见。n%M 函数nfunction dy=cwf1(t,y)n dy=zeros(2,1);n dy(1)=0.1*y(1)*(1-
19、0.001*y(1)-0.0008*y(2);n dy(2)=0.2*y(2)*(1-0.0012*y(1)-0.001*y(2);n 主程序T,X=ode45(cwf1,0 200,200 200);T,Y=ode45(cwf1,0 200,500 200);T,Z=ode45(cwf1,0 200,1200 500);Plot(X(:,1),X(:,2),(Y(:,1),Y(:,2),(Z(:,1),Z(:,2)n 模型模型两种群互惠模型两种群互惠模型研究多个种群之间相互依赖、共生现象。模型讨论模型讨论模型有三个平衡点，分模型有三个平衡点，分别为别为两种群最终达到稳定平衡态两种群共生P3
20、为正平衡点P3 稳定两种群最终达不到稳定平衡态P3 不稳定不稳定%m程序nfunction dy=cwf2(t,y)n dy=zeros(2,1);n dy(1)=0.1*y(1)*(1-0.001*y(1)+0.0005*y(2);n dy(2)=0.2*y(2)*(-1+0.0015*y(1)-0.001*y(2);nT,X=ode45(cwf2,0,100,1600,2800);nT,Y=ode45(cwf2,0,100,1200,2500);nT,Z=ode45(cwf2,0,100,2500,2200);nPlot(X(:,1),X(:,2),Y(:,1),Y(:,2),Z(:,1)
21、,Z(:,2)nText(2000,2600,2000,2600,sigma_11,sigma_1sigma_21)6.6常微分方程的数值解及实验常微分方程的数值解及实验（一）常微分方程数值解的定义一）常微分方程数值解的定义 在生产和科研中所处理的微分方程往往很复杂且大多得不出一般解。而在实际上对初值问题，一般是要求得到解在若干个点上满足规定精确度的近似值，或者得到一个满足精确度要求的便于计算的表达式。因此，研究常微分方程的数值解法是十分必要的因此，研究常微分方程的数值解法是十分必要的。欧拉公式1、用差商代替导数、用差商代替导数 若步长h较小，则有故有公式：此即欧拉法（欧拉法（向前欧拉法向前欧
22、拉法）。2、使用数值积分、使用数值积分对方程y=f(x,y),两边由xi到xi+1积分，并利用梯形公式，有：实际应用时，与欧拉公式结合使用：此即改进的欧拉法改进的欧拉法。故有公式：例1 求解初值问题龙格龙格库特方法库特方法 n考虑微分中值定理 2阶龙格库特公式阶龙格库特公式n n利用泰勒展式得到 4阶龙格库特公式阶龙格库特公式n n 一阶方程组与高阶方程一阶方程组与高阶方程n n Volterra方程解曲线方程解曲线 与与Volterra方程相轨线方程相轨线3、使用泰勒公式、使用泰勒公式 以此方法为基础，有龙格龙格-库塔法库塔法、线性多步法线性多步法等方法。4、数值公式的精度、数值公式的精度
23、当一个数值公式的截断误差可表示为O（hk+1）时（k为正整数，h为步长），称它是一个k阶公式阶公式。k越大，则数值公式的精度越高。欧拉法是一阶公式，改进的欧拉法是二阶公式。龙格-库塔法有二阶公式和四阶公式。线性多步法有四阶阿达姆斯外插公式和内插公式。（三）用三）用Matlab软件求常微分方程的数值解软件求常微分方程的数值解t，x=solver（f,ts,x0,options）ode45 ode23 ode113ode15sode23s由待解方程写成的m-文件名ts=t0，tf，t0、tf为自变量的初值和终值函数的初值ode23：组合的2/3阶龙格-库塔-芬尔格算法ode45：运用组合的4/5阶
24、龙格-库塔-芬尔格算法自变量值函数值用于设定误差限(缺省时设定相对误差10-3,绝对误差10-6),命令为：options=odeset（reltol,rt,abstol,at）,rt，at：分别为设定的相对误差和绝对误差.1、在解n个未知函数的方程组时，x0和x均为n维向量，m-文件中的待解方程组应以x的分量形式写成.2、使用Matlab软件求数值解时，高阶微分方程必须等价地变换成一阶微分方程组.注意注意:解解:令 y1=x，y2=y11、建立m-文件vdp1.m如下：function dy=vdp(t,y)dy=zeros(2,1);dy(1)=y(2);dy(2)=(1-y(1)2)*y(2)-y(1);2、取t0=0，tf=20，输入命令：T,Y=ode15s(vdp1,0 20,2 0);plot(T,Y(:,1),-)3、结果如图例例则微分方程变为一阶微分方程组：