常微分方程模型

1、6.1 人口增长模型,英国人口学家Malthus (1766-1834,模型假设,人口自然增长率 r 为常数即单位时间内人口的增长量与当时的人口呈正比,模型建立,1.指数增长模型,人口以几何级数增加,模型分析,人口将按指数规律无限增长,人口将始终保持不变,人口将按指数规律减少直至绝灭,人口倍增时间,模型求解,Malthus模型预测美国人口,Malthus模型预测美国人口,Malthus模型预测的优缺点,优点,短期预报比较准确,缺点,不适合中长期预报,原因,预报时假设人口增长率 r 为常数。没有考虑环境对人口增长的制约作用,2.阻滞增长模型,假设人口增长率 r(t) 是 t 时刻人口 x(t)

2、的减函数,其中，xm 为考虑到受自然资源和环境条件限制所能容纳的最大人口数量,称最大人口容量,模型假设,模型建立,模型分析（定性分析,人口将递减并趋向于xm,人口将始终保持xm不变,人口将递增并趋向于xm,无论在哪种情况下，人口最终将趋向于最大人口容量,模型求解,人口增长率达到最大值,阻滞增长模型预测美国人口,阻滞增长模型预测美国人口,阻滞增长模型预测的优缺点,优点,中期预报比较准确,缺点,理论上很好，实用性不强,原因,预报时假设固有人口增长率 r 以及最大人口容量 xm 为定值,实际上这两个参数（特别是 xm ）很难确定，而且会随着社会发展情况变化而变化,前面图中曲线末端分叉就是由于这个原因

3、,6.2 药物在体内的分布与排除,药物进入机体形成血药浓度(单位体积血液的药物量,血药浓度需保持在一定范围内给药方案设计,药物在体内吸收、分布和排除过程 药物动力学,建立房室模型药物动力学的基本步骤,房室机体的一部分，药物在一个房室内均匀分布(血药浓度为常数)，在房室间按一定规律转移,本节讨论二室模型中心室(心、肺、肾等)和周边室(四肢、肌肉等,模型假设,中心室(1)和周边室(2),容积不变,药物在房室间转移速率及向体外排除 速率，与该室血药浓度成正比,药物从体外进入中心室，在二室间 相互转移,从中心室排出体外,模型建立,线性常系数非齐次方程,对应齐次方程通解,模型建立,几种常见的给药方式,1

4、.快速静脉注射,t=0 瞬时注射剂量d 的药物进入中心室,血药浓度立即为 d/V1,给药速率 f (t) 和初始条件,2.恒速静脉滴注,t T时, c1(t)和 c2(t)按指数规律衰减趋于零,3.口服或肌肉注射,相当于药物( 剂量d)先进入吸收室，吸收后再进入中心室,吸收室药量x0(t,参数估计,各种给药方式下的 c1(t), c2(t) 取决于参数k12, k21, k13, V1,V2,以快速静脉注射为例 ,在ti(i=1,2,n)测得c1(ti,由较大的 用最小二乘法定A,由较小的 用最小二乘法定B,参数估计法一,进入中心室的药物全部排除,参数 估计法二, 构造非线性拟合函数 TWOE

5、XPS.M function E=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,a0=100 1 1 ; options=optimset(fminsearch); options.TolX=0.01; options.Display=off; a=fminsearch(ps,a0,options,x,y,a= 112.2378 0.1823 2.1773,6.3 传染病模型,问题,描述传染病的传播过程,分析受感染人数的变化规律,预报传染病高潮到来的时刻,预防传染病蔓延的手段,

6、按照传播过程的一般规律，用机理分析方法建立模型,已感染人数 (病人) i(t,每个病人每天有效接触(足以使人致病)人数为,模型1,假设,若有效接触的是病人，则不能使病人数增加,建模,模型2,区分已感染者(病人)和未感染者(健康人,假设,1）总人数N不变，病人和健康 人的 比例分别为,2）每个病人每天有效接触人数为, 且使接触的健康人致病,建模,日 接触率,SI 模型,模型2,tm传染病高潮到来时刻,(日接触率) tm,病人可以治愈,t=tm, di/dt 最大,模型3,传染病无免疫性病人治愈成为健康人，健康人可再次被感染,增加假设,SIS 模型,3）病人每天治愈的比例为,日治愈率,建模,日接触

7、率,1/ 感染期,一个感染期内每个病人的有效接触人数，称为接触数,模型3,接触数 =1 阈值,感染期内有效接触感染的健康者人数不超过病人数,模型2(SI模型)如何看作模型3(SIS模型)的特例,模型4,传染病有免疫性病人治愈后即移出感染系统，称移出者,SIR模型,假设,1）总人数N不变，病人、健康人和移出者的比例分别为,2）病人的日接触率 , 日治愈率, 接触数 = /,建模,需建立 的两个方程,模型4,SIR模型,模型4,SIR模型,相轨线 的定义域,在D内作相轨线 的图形，进行分析,模型4,SIR模型,相轨线 及其分析,s(t)单调减相轨线的方向,P1: s01/ i(t)先升后降至0,P

8、2: s01/ i(t)单调降至0,1/ 阈值,模型4,SIR模型,预防传染病蔓延的手段,(日接触率) 卫生水平,日治愈率) 医疗水平,传染病不蔓延的条件s01,的估计,降低 s0,提高 r0,提高阈值 1,SIR模型,被传染人数的估计法一,记被传染人数比例,小, s0 1,提高阈值 1/ 降低被传染人数比例 x,s0 - 1/ =,被传染人数的估计法二,X=fzero(x-1.2*log(x/0.96)-0.99,0.5) X=0.8651,6.4多种群生态数学模型,意大利生物学家DAncona曾致力于鱼类种群相互制约关系的研究，他从第一次世界大战期间,地中海各港口捕获的几种鱼类捕获量百分比

9、的资料中，发现鲨鱼等的比例有明显增加（见下表），而供其捕食的食用鱼的百分比却明显下降.显然战争使捕鱼量下降，食用鱼增加，鲨鱼等也随之增加，但为何鲨鱼的比例大幅增加呢,他无法解释这个现象，于是求助于其岳父，著名的意大利数学家V.Volterra，希望建立一个食饵捕食系统的数学模型，定性或定量地回答这个问题,该模型反映了在没有人工捕获的自然环境中食饵与捕食者之间的制约关系，没有考虑食饵和捕食者自身的阻滞作用，是Volterra提出的最简单的模型,模型（一）不考虑捕获,定理 Volterra微分方程组对应初值问题,的解,是周期函数，且解的周期平均值为,首先，建立m-文件shier.m如下： func

10、tion 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(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增为

11、c+e,设战前捕获能力系数e=0.3, 战争中降为e=0.1, 则战前与战争中的模型分别为,Volterra原理,模型求解,1、分别用m-文件shier1.m和shier2.m定义上述两个方程,2、建立主程序shark1.m, 求解两个方程，并画出两种情况下鲨鱼数在鱼类总数中所占比例 x2(t)/x1(t)+x2(t,To Matlab(shark1,实线为战前的鲨鱼比例，“*”线为战争中的鲨鱼比例,结论：战争中鲨鱼的比例比战前高,注：ts 中终值(=15)和步长=(0.1)的确定,计算结果（数值，图形,x(t),y(t)是周期函数，相图(x,y)是封闭曲线,x(t),y(t)的周期约为10.

12、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.Walsh（1978,外界有食饵迁入,外界有食饵迁出,也可以表示人工干预，如投放或捕获,模型讨论,食饵捕食者模型,模型2,考虑食饵种群内部存在生存竞争 G.Bojadziev,表示当没有捕食者存在时食饵种群的环境容纳量,模型3,考虑食饵和捕食者种群内部都存在生存竞争 张锦炎（1979,表示当没有捕食者存在时食饵

13、种群的环境容纳量,模型4,考虑双方内部都存在生存竞争，且捕食者另有食物来源 E.C.Pielou,平衡点,模型,两种群竞争模型,竞争模型,竞争排斥原理（Competition Exclution Law,多个种群依靠同一个生存资源而生活，如果生活在同一个地理空间，猎取相同食物或营养物。在有限的相同生存资源条件下，如果存在竞争关系，它们必然相互排斥，展开激烈的生存竞争。结局是竞争力较弱的种群灭绝，竞争力最强的种群达到其环境容纳量,模型讨论,竞争排斥原理 是针对前三种情形得出的结论， 第四种情况极为罕见, M 函数 function dy=cwf1 (t,y) dy=zeros(2,1); dy(

14、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,主程序 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,模型,两种群互惠模型,互惠模型,研究多个种群之间相互依赖、共生现象,模型讨论,模型有三个平衡点，分别为,两种群最终达到稳定平衡态,两种群共生,P3

15、为正平衡点,P3 稳定,两种群最终达不到稳定平衡态,P3 不稳定,m程序,function dy=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,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(:

16、,2) Text(2000,2600, 2000,2600,sigma_11,sigma_1sigma_21,6.6常微分方程的数值解及实验,一）常微分方程数值解的定义,在生产和科研中所处理的微分方程往往很复杂且大多得不出一般解。而在实际上对初值问题，一般是要求得到解在若干个点上满足规定精确度的近似值，或者得到一个满足精确度要求的便于计算的表达式,因此，研究常微分方程的数值解法是十分必要的,欧拉公式,1、用差商代替导数,若步长h较小，则有,故有公式,此即欧拉法（向前欧拉法,2、使用数值积分,对方程y=f(x,y), 两边由xi到xi+1积分，并利用梯形公式，有,实际应用时，与欧拉公式结合使用,

17、此即改进的欧拉法,故有公式,例1 求解初值问题,龙格库特方法,考虑微分中值定理,2阶龙格库特公式,利用泰勒展式得到,4阶龙格库特公式,一阶方程组与高阶方程,Volterra方程解曲线 与Volterra方程相轨线,3、使用泰勒公式,以此方法为基础，有龙格-库塔法、线性多步法等方法,4、数值公式的精度,当一个数值公式的截断误差可表示为O（hk+1）时（k为正整数，h为步长），称它是一个k阶公式,k越大，则数值公式的精度越高,欧拉法是一阶公式，改进的欧拉法是二阶公式。 龙格-库塔法有二阶公式和四阶公式。 线性多步法有四阶阿达姆斯外插公式和内插公式,三）用Matlab软件求常微分方程的数值解,t，x=solver（f,ts,x0,options,1、在解n个未知函数的方程组