数学建模讲义-微分方程模型

数学建模讲义微分方程模型动态模型

描述对象特征随时间(空间)的演变过程

分析对象特征的变化规律

预报对象特征的未来性态

研究控制对象特征的手段

根据函数及其变化率之间的关系确定函数微分方程建模

根据建模目的和问题分析作出简化假设

按照内在规律或用类比法建立微分方程主要内容生物单种群增长模型

3.1

人口增长模型

3.2

传染病模型生物多种群增长模型

3.3

正规战与游击战

3.4

捕食系统的Volterra方程

为了保持自然资料的合理开发与利用，人类必须保持并控制生态平衡，甚至必须控制人类自身的增长。本节将建立几个简单的单种群增长模型，以简略分析一下这方面的问题。一般生态系统的分析可以通过一些简单模型的复合来研究，大家若有兴趣可以根据生态系统的特征自行建立相应的模型。

美丽的大自然

种群的数量本应取离散值，但由于种群数量一般较大，为建立微分方程模型，可将种群数量看作连续变量，甚至允许它为可微变量，由此引起的误差将是十分微小的。离散化为连续，方便研究3.1

如何预报人口的增长

--Malthus模型与Logistic模型背景

年1625183019301960197419871999人口(亿)5102030405060世界人口增长概况中国人口增长概况

年19081933195319641982199019952000人口(亿)3.04.76.07.210.311.312.013.0研究人口变化规律控制人口过快增长指数增长模型——马尔萨斯提出(1798)常用的计算公式x(t):时刻t的人口基本假设

:人口(相对)增长率r

是常数,不考虑移民今年人口x0,年增长率rk年后人口随着时间增加，人口按指数规律无限增长模型检验

比较历年的人口统计资料，可发现人口增长的实际情况与马尔萨斯模型的预报结果基本相符，例如，1961年世界人口数为30.6（即3.06×109），人口增长率约为2%，人口数大约每35年增加一倍。检查1700年至1961的260年人口实际数量，发现两者几乎完全一致，且按马氏模型计算，人口数量每34.6年增加一倍，两者也几乎相同。模型预测

假如人口数真能保持每34.6年增加一倍，那么人口数将以几何级数的方式增长。例如，到2510年，人口达2×1014个，即使海洋全部变成陆地，每人也只有9.3平方英尺的活动范围，而到2670年，人口达36×1015个，只好一个人站在另一人的肩上排成二层了。故马尔萨斯模型是不完善的。几何级数的增长Malthus模型实际上只有在群体总数不太大时才合理，到总数增大时，生物群体的各成员之间由于有限的生存空间，有限的自然资源及食物等原因，就可能发生生存竞争等现象。所以Malthus模型假设的人口净增长率不可能始终保持常数，它应当与人口数量有关。指数增长模型的应用及局限性

与19世纪以前欧洲一些地区人口统计数据吻合

适用于19世纪后迁往加拿大的欧洲移民后代

可用于短期人口增长预测

不符合19世纪后多数地区人口增长规律

不能预测较长期的人口增长过程19世纪后人口数据人口增长率r不是常数(逐渐下降)阻滞增长模型(Logistic模型)人口增长到一定数量后，增长率下降的原因：资源、环境等因素对人口增长的阻滞作用且阻滞作用随人口数量增加而变大假设r~固有增长率(x很小时)xm~人口容量（资源、环境能容纳的最大数量）r是x的减函数dx/dtx0xmxm/2xmtx0x(t)~S形曲线,x增加先快后慢x0xm/2阻滞增长模型(Logistic模型)参数估计用指数增长模型或阻滞增长模型作人口预报，必须先估计模型参数r或r,xm

利用统计数据用最小二乘法作拟合例：美国人口数据（单位~百万）186018701880……196019701980199031.438.650.2……179.3204.0226.5251.4专家估计阻滞增长模型(Logistic模型)r=0.2557,xm=392.1继续最小二乘法设经实际测量已得到n组数据（xi,yi），i=1,…,n。将数据画在平面直角坐标系中，见图。如果建模者判断这n个点很象是分布在某条直线附近，令该直线方程为y=ax+b，进而利用数据来求参数a和b。由于该直线只是数据近似满足的关系式，故yi-(axi+b)=0一般不成立，但我们希望最小此式对a和b的偏导数均为0，解相应方程组，求得：y=ax+byO(xi,yi)x其中和分别为xi和yi的平均值

如果建模者判断变量间的关系并非线性关系而是其他类型的函数，则可作变量替换使之转化为线性关系或用类似方法拟合。用MATLAB作线性最小二乘拟合1.作多项式f(x)=a1xm+…+amx+am+1拟合,可利用已有程序:a=polyfit(x,y,m)输出拟合多项式系数a=[a1,…am,

am+1](数组)）输入同长度的数组X，Y拟合多项式次数1.lsqcurvefit已知数据点：xdata=（xdata1，xdata2，…，xdatan），

ydata=（ydata1，ydata2，…，ydatan）

用MATLAB作非线性最小二乘拟合Matlab的提供了两个求非线性最小二乘拟合的函数：lsqcurvefit和lsqnonlin。两个命令都要先建立M-文件fun.m，在其中定义函数f(x)，但两者定义f(x)的方式是不同的,可参考例题.

lsqcurvefit用以求含参量x（向量）的向量值函数F(x,xdata)=（F（x，xdata1），…，F（x，xdatan））T中的参变量x(向量),使得

输入格式为:(1)x=lsqcurvefit(‘fun’,x0,xdata,ydata);(2)x=lsqcurvefit(‘fun’,x0,xdata,ydata,options);(3)x=lsqcurvefit(‘fun’,x0,xdata,ydata,options,’grad’);(4)[x,options]=lsqcurvefit(‘fun’,x0,xdata,ydata,…);(5)[x,options,funval]=lsqcurvefit(‘fun’,x0,xdata,ydata,…);(6)[x,options,funval,Jacob]=lsqcurvefit(‘fun’,x0,xdata,ydata,…);fun是一个事先建立的定义函数F(x,xdata)

的M-文件,自变量为x和xdata说明：x=lsqcurvefit(‘fun’,x0,xdata,ydata,options);迭代初值已知数据点选项见无约束优化

lsqnonlin用以求含参量x（向量）的向量值函数

f(x)=(f1(x),f2(x),…,fn(x))T

中的参量x，使得

最小。其中fi（x）=f（x，xdatai，ydatai）

=F(x,xdatai)-ydatai

2.lsqnonlin已知数据点：xdata=（xdata1，xdata2，…，xdatan）

ydata=（ydata1，ydata2，…，ydatan）输入格式为：

1）x=lsqnonlin（‘fun’，x0）；

2）x=lsqnonlin（‘fun’，x0，options）；

3）x=lsqnonlin（‘fun’，x0，options，‘grad’）；

4）[x，options]=lsqnonlin（‘fun’，x0，…）；

5）[x，options，funval]=lsqnonlin（‘fun’，x0，…）；说明：x=lsqnonlin

（‘fun’，x0，options）；fun是一个事先建立的定义函数f(x)的M-文件，自变量为x迭代初值选项有无约束优化

例用下面一组数据拟合中的参数a，b，k该问题即解最优化问题：

1）编写M-文件curvefun1.m

functionf=curvefun1(x,tdata)f=x(1)+x(2)*exp(-0.02*x(3)*tdata)%其中x(1)=a;x(2)=b；x(3)=k;2）输入命令tdata=100:100:1000cdata=1e-03*[4.54,4.99,5.35,5.65,5.90,6.10,6.26,6.39,6.50,6.59];x0=[0.2,0.05,0.05];x=lsqcurvefit('curvefun1',x0,tdata,cdata)f=curvefun1(x,tdata)

F(x，tdata)=

，x=(a，b，k)解法1.用命令lsqcurvefit3）运算结果为：f=0.00430.00510.00560.00590.00610.00620.00620.00630.00630.0063x=0.0063-0.00340.25424）结论：a=0.0063,b=-0.0034,k=0.2542返回模型检验用模型计算2000年美国人口，与实际数据比较实际为281.4(百万)模型应用——预报美国2010年的人口加入2000年人口数据后重新估计模型参数Logistic模型在经济领域中的应用(如耐用消费品的售量)阻滞增长模型(Logistic模型)r=0.2490,xm=434.0x(2010)=306.0模型检验

用Logistic模型来描述种群增长的规律效果如何呢？1945年克朗皮克（Crombic）做了一个人工饲养小谷虫的实验，数学生物学家高斯（E·F·Gauss）也做了一个原生物草履虫实验，实验结果都和Logistic曲线十分吻合。

大量实验资料表明用Logistic模型来描述种群的增长，效果还是相当不错的。例如，高斯把5只草履虫放进一个盛有0.5cm3营养液的小试管，他发现，开始时草履虫以每天230.9%的速率增长，此后增长速度不断减慢，到第五天达到最大量375个，实验数据与r=2.309，a=0.006157，N(0)=5的Logistic曲线：

几乎完全吻合，见图3.6。

图3-6Malthus模型和Logistic模型的总结

Malthus模型和Logistic模型均为对微分方程（3.7）所作的模拟近似方程。前一模型假设了种群增长率r为一常数，（r被称为该种群的内禀增长率）。后一模型则假设环境只能供养一定数量的种群，从而引入了一个竞争项。

用模拟近似法建立微分方程来研究实际问题时必须对求得的解进行检验，看其是否与实际情况相符或基本相符。相符性越好则模拟得越好，否则就得找出不相符的主要原因，对模型进行修改。Malthus模型与Logistic模型虽然都是为了研究种群数量的增长情况而建立的，但它们也可用来研究其他实际问题，只要这些实际问题的数学模型有相同的微分方程即可，下面我们来看两个较为有趣的实例。

年龄分布对于人口预测的重要性

只考虑自然出生与死亡，不计迁移人口发展方程人口发展方程一阶偏微分方程例2:传染病模型问题

描述传染病的传播过程

分析受感染人数的变化规律

预报传染病高潮到来的时刻

预防传染病蔓延的手段

按照传播过程的一般规律，用机理分析方法建立模型

已感染人数(病人)i(t)

每个病人每天有效接触(足以使人致病)人数为

模型1假设若有效接触的是病人，则不能使病人数增加必须区分已感染者(病人)和未感染者(健康人)建模?模型2区分已感染者(病人)和未感染者(健康人)假设1）总人数N不变，病人和健康人的比例分别为2）每个病人每天有效接触人数为

,且使接触的健康人致病建模

~日接触率SI模型模型21/2tmii010ttm~传染病高潮到来时刻

(日接触率)tm

Logistic模型病人可以治愈！?t=tm,di/dt最大模型3传染病无免疫性——病人治愈成为健康人，健康人可再次被感染增加假设SIS模型3）病人每天治愈的比例为

~日治愈率建模

~日接触率1/

~感染期

~一个感染期内每个病人的有效接触人数，称为接触数。1-1/

idi/dt01>1i0i00ti>11-1/

i0i0t

1di/dt<0接触数

=1~阈值如果感染期内有效接触感染的健康者人数不超过病人数,患者就会全部治愈。模型4传染病有免疫性——病人治愈后即移出感染系统，称移出者SIR模型假设1）总人数N不变，病人、健康人和移出者的比例分别为2）病人的日接触率

,日治愈率

,

接触数

=/建模需建立的两个方程模型4SIR模型无法求出的解析解在相平面上研究解的性质模型4消去dtSIR模型相轨线的定义域相轨线11si0D在D内作相轨线的图形，进行分析si101DSIR模型相轨线及其分析传染病蔓延传染病不蔓延s(t)单调减

相轨线的方向P1s0imP1:s0>1/

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

i(t)单调降至01/~阈值P3P4P2S0模型4SIR模型预防传染病蔓延的手段

(日接触率)卫生水平

(日治愈率)

医疗水平传染病不蔓延的条件——s0<1/

降低s0提高r0

提高阈值1/

降低

(=

/

)

,

群体免疫

的估计模型4SIR模型被传染人数的估计记被传染人数比例x<<s0i0P1

i0

0,s0

1

小,s0

1提高阈值1/

降低被传染人数比例xs0-1/

=

例3：正规战与游击战战争分类：正规战争，游击战争，混合战争只考虑双方兵力多少和战斗力强弱兵力因战斗及非战斗减员而减少，因增援而增加战斗力与射击次数及命中率有关建模思路和方法为用数学模型讨论社会领域的实际问题提供了可借鉴的示例第一次世界大战Lanchester提出预测战役结局的模型一般模型

每方战斗减员率取决于双方的兵力和战斗力

每方非战斗减员率与本方兵力成正比

甲乙双方的增援率为u(t),v(t)f,g

取决于战争类型x(t)~甲方兵力，y(t)~乙方兵力模型假设模型正规战争模型

甲方战斗减员率只取决于乙方的兵力和战斗力双方均以正规部队作战

忽略非战斗减员

假设没有增援f(x,y)=

ay,a~乙方每个士兵的杀伤率a=rypy,ry~射击率，

py~命中率0正规战争模型为判断战争的结局，不求x(t),y(t)而在相平面上讨论x与y的关系平方律模型乙方胜游击战争模型双方都用游击部队作战

甲方战斗减员率还随着甲方兵力的增加而增加

忽略非战斗减员

假设没有增援f(x,y)=

cxy,c~乙方每个士兵的杀伤率c=rypyry~射击率py~命中率py=sry/sxsx~甲方活动面积sry~乙方射击有效面积0游击战争模型线性律模型0混合战争模型甲方为游击部队，乙方为正规部队乙方必须10倍于甲方的兵力设x0=100,rx/ry=1/2,px=0.1,sx=1(km2),sry=1(m2)微分方程的解析解

求微分方程（组）的解析解命令:dsolve(‘方程1’,‘方程2’,…‘方程n’,‘初始条件’,‘自变量’)

结果：u=tg(t-c)

解

输入命令:y=dsolve('D2y+4*Dy+29*y=0','y(0)=0,Dy(0)=15','x')结果为:y=3e-2xsin（5x）解输入命令：

[x,y,z]=dsolve('Dx=2*x-3*y+3*z','Dy=4*x-5*y+3*z','Dz=4*x-4*y+2*z','t')；

x=simple(x)%将x化简

y=simple(y)z=simple(z)结果为：x=(c1-c2+c3+c2e-3t-c3e-3t)e2t

y=-c1e-4t+c2e-4t+c2e-3t-c3e-3t+c1-c2+c3)e2tz=(-c1e-4t+c2e-4t+c1-c2+c3)e2t

微分方程的数值解（一）常微分方程数值解的定义

在生产和科研中所处理的微分方程往往很复杂且大多得不出一般解。而在实际上对初值问题，一般是要求得到解在若干个点上满足规定精确度的近似值，或者得到一个满足精确度要求的便于计算的表达式。因此，研究常微分方程的数值解法是十分必要的。（二）建立数值解法的一些途径1、用差商代替导数

若步长h较小，则有故有公式：此即欧拉法。2、使用数值积分对方程y’=f(x,y),两边由xi到xi+1积分，并利用梯形公式，有：实际应用时，与欧拉公式结合使用：此即改进的欧拉法。故有公式：3、使用泰勒公式

以此方法为基础，有龙格-库塔法、线性多步法等方法。4、数值公式的精度

当一个数值公式的截断误差可表示为O（hk+1）时（k为正整数，h为步长），称它是一个k阶公式。k越大，则数值公式的精度越高。欧拉法是一阶公式，改进的欧拉法是二阶公式。龙格-库塔法有二阶公式和四阶公式。线性多步法有四阶阿达姆斯外插公式和内插公式。（三）用Matlab软件求常微分方程的数值解[t，x]=solver（’f’,ts,x0,options）ode45ode23ode113ode15sode23s由待解方程写成的m-文件名ts=[t0，tf]，t0、tf为自变量的初值和终值函数的初值ode23：组合的2/3阶龙格-库塔-芬尔格算法ode45：运用组合的4/5阶龙格-库塔-芬尔格算法自变量值函数值用于设定误差限(缺省时设定相对误差10-3,绝对误差10-6),命令为：options=odeset（’reltol’,rt,’abstol’,at）,rt，at：分别为设定的相对误差和绝对误差.

1、在解n个未知函数的方程组时，x0和x均为n维向量，m-文件中的待解方程组应以x的分量形式写成.

2、使用Matlab软件求数值解时，高阶微分方程必须等价地变换成一阶微分方程组.注意:解:令y1=x，y2=y1’1、建立m-文件vdp1000.m如下：

functiondy=vdp1000(t,y)dy=zeros(2,1);dy(1)=y(2);dy(2)=1000*(1-y(1)^2)*y(2)-y(1);

2、取t0=0，tf=3000，输入命令：

[T,Y]=ode15s('vdp1000',[03000],[20]);plot(T,Y(:,1),'-')3、结果如图解

1、建立m-文件rigid.m如下：

functiondy=rigid(t,y)dy=zeros(3,1);dy(1)=y(2)*y(3);dy(2)=-y(1)*y(3);dy(3)=-0.51*y(1)*y(2);2、取t0=0，tf=12，输入命令：

[T,Y]=ode45('rigid',[012],[011]);plot(T,Y(:,1),'-',T,Y(:,2),'*',T,Y(:,3),'+')3、结果如图图中，y1的图形为实线，y2的图形为“*”线，y3的图形为“+”线.稳定性问题

在研究许多实际问题时，人们最为关心的也许并非系统与时间有关的变化状态，而是系统最终的发展趋势。例如，在研究某频危种群时，虽然我们也想了解它当前或今后的数量，但我们更为关心的却是它最终是否会绝灭，用什么办法可以拯救这一种群，使之免于绝种等等问题。要解决这类问题，需要用到微分方程或微分方程组的稳定性理论。在下两节，我们将研究几个与稳定性有关的问题。一般的微分方程或微分方程组可以写成：定义称微分方程或微分方程组为自治系统或动力系统。（3.28）

若方程或方程组f(x)=0有解Xo，X=Xo显然满足（3.28）。称点Xo为微分方程或微分方程组（3.28)的平衡点或奇点。例7Logistic模型

共有两个平衡点：N=0和N=K，分别对应微分方程的两两个特殊解。前者为No=0时的解而后者为No=K时的解。

当No<K时，积分曲线N=N(t)位于N=K的下方；当No>K时，则位于N=K的上方。从图3-17中不难看出，若No>0，积分曲线在N轴上的投影曲线（称为轨线）将趋于K。这说明，平衡点N=0和N=K有着极大的区别。图3-17

定义1自治系统的相空间是指以（x1,…,xn）为坐标的空间Rn。

特别，当n=2时，称相空间为相平面。空间Rn的点集{(x1,…,xn)}|xi=xi(t)满足(3.28)，i=1,…,n}称为系统的轨线，所有轨线在相空间的分布图称为相图。定义2设x0是（3.28）的平衡点，称：

（1）x0是稳定的，如果对于任意的ε>0，存在一个δ>0，只要|x(0)-x0|<δ，就有|x(t)-x0|<ε对所有的t都成立。（2）x0是渐近稳定的，如果它是稳定的且。

微分方程平衡点的稳定性除了几何方法，还可以通过解析方法来讨论，所用工具为以下一些定理。（3）x0是不稳定的，如果（1）不成立。根据这一定义，Logistic方程的平衡点N=K是稳定的且为渐近稳定的，而平衡点N=0则是不稳定的。解析方法定理1设xo是微分方程的平衡点：若,则xo是渐近稳定的若,则xo是渐近不稳定的证由泰勒公式，当x与xo充分接近时，有：

由于xo是平衡点，故f(xo)=0。若，则当x<xo时必有f(x)>0，从而x单增；当x>xo时，又有f(x)<0，从而x单减。无论在哪种情况下都有x→xo，故xo是渐进稳定的。的情况可类似加以讨论。高阶微分方程与高阶微分方程组平衡点的稳定性讨论较为复杂，大家有兴趣可参阅微分方程定性理论。为了下两节的需要，我们简单介绍一下两阶微分方程组平衡点的稳定性判别方法。考察两阶微分方程组：（3.29）

令，作一坐标平移，不妨仍用x记x’，则平衡点xo的稳定性讨论转化为原点的稳定性讨论了。将f(x1,x2)、g(x1,x2)在原点展开，（3.29）又可写成:考察（3.29）的线性近似方程组:（3.30）其中：记λ1、λ2为A的特征值则λ1、λ2是方程：det（A-λI）=λ2-(a+b)λ+(ad–bc)=0的根令p=a+d,q=ad-bc=|A|，则，记。讨论特征值与零点稳定的关系（1）若△>0，可能出现以下情形：

①若q>0，λ1λ2>0。当p>0时，零点不稳定；当p<0时，零点稳定

若q<0，λ1λ2<0

零点为不稳定的鞍点③q=0，此时λ1=p，λ2=0，零点不稳定。（2）△=0，则λ1=λ2:

当p>0时，零点不稳定当p<0时，零点稳定（2）△<0，此时若a<0，零点稳定,a>0,零点不稳定若a=0，有零点为中心的周期解

综上所述：仅当p<0且q>0时，（3.30）零点才是渐近稳定的；当p=0且q>0时（3.30）有周期解，零点是稳定的中心（非渐近稳定）；在其他情况下，零点均为不稳定的。

非线性方程组（3.29）平衡点稳定性讨论可以证明有下面定理成立:定理2若（3.30）的零点是渐近稳定的，则（3.29）的平衡点也是渐近稳定的；若（3.30）的零点是不稳定的，则（3.29）的平衡点也是不稳定的。捕食系统的Volterra方程

问题背景：

意大利生物学家D’Ancona曾致力于鱼类种群相互制约关系的研究，在研究过程中他无意中发现了一些第一次世界大战期间地中海沿岸港口捕获的几种鱼类占捕获总量百分比的资料，从这些资料中他发现各种软骨掠肉鱼，如鲨鱼、鳐鱼等我们称之为捕食者（或食肉鱼）的一些不是很理想的鱼类占总渔获量的百分比。在1914~1923年期间，意大利阜姆港收购的鱼中食肉鱼所占的比例有明显的增加：年代19141915191619171918百分比11.921.422.121.236.4年代19191920192119221923百分比27.316.015.914.810.7

他知道，捕获的各种鱼的比例近似地反映了地中海里各种鱼类的比例。战争期间捕鱼量大幅下降，但捕获量的下降为什么会导致鲨鱼、鳐鱼等食肉鱼比例的上升，即对捕食者有利而不是对食饵有利呢？他百思不得其解，无法解释这一现象，就去求教当时著名的意大利数学家V.Volterra，希望他能建立一个数学模型研究这一问题。Volterra将鱼划分为两类。一类为食用鱼（食饵），数量记为x1(t)，另一类为食肉鱼（捕食者），数量记为x2(t)，并建立双房室系统模型。1、模型建立

大海中有食用鱼生存的足够资源，可假设食用鱼独立生存将按增长率为r1的指数律增长（Malthus模型），既设：

由于捕食者的存在，食用鱼数量因而减少，设减少的速率与两者数量的乘积成正比（竞争项的统计筹算律），即：对于食饵（Prey）系统：λ1反映了捕食者掠取食饵的能力对于捕食者（Predator）系统：捕食者设其离开食饵独立存在时的死亡率为r2，即：但食饵提供了食物，使生命得以延续。这一结果也要通过竞争来实现，再次利用统计筹算律，得到：综合以上分析，建立P-P模型（Volterra方程）的方程组：（3.31）方程组（3.31）反映了在没有人工捕获的自然环境中食饵与捕食者之间的相互制约关系。下面我们来分析该方程组。2、模型分析

方程组（3.31）是非线性的，不易直接求解。容易看出，该方程组共有两个平衡点，即：Po(0,0)是平凡平衡点且明显是不稳定，没必要研究方程组还有两组平凡解：

和和所以x1、x2轴是方程组的两条相轨线。

当x1(0)、x2(0)均不为零时，，应有x1(t)>0且x2(t)>0，相应的相轨线应保持在第一象限中。求（