python中绘制&方向场和积分曲线族 大多使用matlab或Maple绘制的。 A Course in Ordinary Differential Equations &作者：Randall J. Swift，Stephen A. &&& from scipy import * &&& from scipy import integrate &&& from scipy.integrate import ode &&& import numpy as np &&& import matplotlib.pyplot as plt &&& fig=plt.figure(num=1) &&& ax = fig.add_subplot(111) &&& def vf(t,x): dx = np.zeros(2) dx[1]=x[0]**2-x[0]-2 &&& t0=0; tEnd=10; dt=0.01; r=ode(vf).set_integrator("vode",method='bdf',max_step=dt) &&& ic=[[-3.5,-10],[-3,-10],[-2.5,-10]] &&& color=['r','b','g'] &&& for k in range(len(ic)): Y=[];T=[];S=[]; r.set_initial_value(ic[k],t0).set_f_params() while r.successful() and r.t +dt & tEnd: ...&&&&&&&& r.integrate(r.t+dt) ...&&&&&&&& Y.append(r.y) S=np.array(np.real(Y)) ax.plot(S[:,0],S[:,1],color=color[k],lw=1.25) 另一个微分方程y'=y+x 的积分曲线族。y=C*e^x-x-1，但是其积分曲线与这个原函数无关。 #生态学的Lokta-Volterra模型 以上网友发言只代表其个人观点，不代表新浪网的观点或立场。当前位置： > 微分方程方向场MATLAB仿真工具箱设计 【摘要】 围绕微分方程方向场这个核心概念，利用MATLAB求解微分方程的优势，得到了求解微分方程初值解、微分方程方向场等各种算法的程序实现，实现了MATLAB仿真工具箱设计，利用图形用户界面（GUI）方法，设计了良好的人-机交互系统的主界面，最后给出了实际例子的程序运行结果，对推动微分方程在具体实践中的应用和普及，具有实际意义。　　【关键词】　　　本文出自：中国原创论文网 微分方程；方向场；MATLAB ；仿真；工具箱　　【　　【　　0 引言　　微分方程是近代数学的一个十分重要的分支，无论在工程技术上、自动控制理论、物理等自然科学领域，还是在经济、金融、保险等社会科学领域中，都有着广泛的应用。随着互联网和计算机技术的迅速发展，数学也需要与时俱进，特别是把计算机技术应用到微分方程体系中，是微分方程发展的一个新方向。　　微分方程方向场是微分方程的一个重要概念，它给出了微分方程解的一种几何解释。因为绝大多数的微分方程是不能用初等积分法求解的，所以由这些几何解释，可以从微分方程本身的特性了解到它的任一解所应具有的某些几何特征，无疑在理论研究和实际应用中，都有着重要的意义。　　MATLAB是当前最流行的、功能强大的高效率的数值计算的可视化科技应用软件和编程语言之一。尽管MATLAB具有很强大的计算功能，但是，MATLAB里仍然没有微分方程方向场方面的函数和工具箱，考虑到MATLAB具有许多现成的求解微分方程及其他数学运算的函数及其解释性语言的特点，结合考虑微分方程方向场的自身特点，利用MATLAB开发微分方程方向场仿真工具箱能够发挥二者各自的优势，它比用VC++等其他计算机语言具有事半功倍的效果。微分方程方向场的工具软件国内尚不多见，国外有DFIELD等，但在国内并不普及。为了促进微分方程理论更广泛地应用到具体实践中，利用现有的算法，开发一个成功的微分方程方向场系统具有重大实际意义。　　1 微分方程方向场简介　　2 MDFIELD系统实现　　微分方程方向场主要是借助一种简单的几何方法来研究已知的微分方程y'=f(x,y)的解，直观上给出了微分方程解曲线的一般形状，通过每一点解曲线的切线方向都与方向场中的线段接近平行，从初始点(x,y)开始，可以通过方向场画出近似的解曲线与可见的斜率线段尽可能平行，主要函数如下。　　2.1 MATLAB微分方程求解函数介绍　　MATLAB的70余个工具箱覆盖面极广，包括丰富的数值分析、矩阵运算、图形绘制、数据处理、信号处理、图像处理、小波分析、鲁棒控制、系统辨别、非线性控制、模糊逻辑、神经网络、优化理论、统计分析等。由于它的功能强大且应用越来越广泛，加之不断开发出来的大量的应用于不同学科的工具箱，使MATLAB越来越受到人们的重视，成为国际控制界广泛使用的语言之一。下面对本开发平台中所使用和涉及的主要的函数或指令的含义进行描述，如表1所示。　　2.2 主要MATLAB程序　　实际应用问题通过数学建模所归纳而得到的方程，绝大多数都是微分方程，真正能得到代数方程的机会很少。另一方面，能够求解的微分方程也是十分有限的，特别是高阶方程和偏微分方程（组）。这就要求我们必须研究微分方程（组）的解法，既要研究微分方程（组）的解析解法（精确解），更要研究微分方程（组）的数值解法（近似解）。MATLAB强大的数据分析和高级可视化软件使其成为许多科研与应用开发的首选平台，本文对常微分方程向量场理论利用MATLAB的函数实现其功能仿真，下面介绍主要MATLAB程序。　　3 图形用户界面　　利用MATLAB程序特点进行微分方程方向场仿真大大减少了编程工作量，对于一个成功的软件来说，其内容和基本功能固然应该是第一位的，但除此之外，图形界面的优劣往往也决定着软件的档次，MATLAB图形用户界面的优势使得系统工具箱形象直观，容易对输入输出、函数方程、取值范围进行修改，也可以通过存储在workspace或文本文件上的数据进行分析和处理。界面结构的设计利用MATLAB的GUIDE功能。　　MDFIELD系统应用平台的界面程序主要内容就是每个控件的回调函数，利用回调函数，把编制的功能分散的微分方程向量场函数整合在如图1的主界面中，下面以控件Proceed为例，给出回调函数Callback的语句。　　4 结论　　本文介绍了如何利用MATLAB的微分方程函数及图形用户界面开发和设计微分方程方向场仿真工具箱的方法和步骤。将MATLAB和微分方程理论有机结合起来，实现了复杂微分方程应用问题系统的设计和高效仿真及图形输出，对普及微分方程向量场理论，推广微分方程在实际中应用和进一步提高教学效果都具有明显的促进作用。　　参考文献　　[1] 东北师范大学微分方程教研室.常微分方程[M].北京:高等教育出版社,.　　[2] 王高雄.常微分方程[M].北京:高等教育出版社,.　　[3] C.Henry Edwards,David E.Penney. Differential Equations and Boundary Value Problems: Computing and Modeling [M]. New York:Pearson College Division, .　　[4] 张雪峰,张庆灵.粗糙集数据分析系统MATLAB仿真工具箱设计[J].东北大学学报,):40-43.　　[5] 罗华飞.MATLAB GUI 设计学习笔记[M].北京:北京航空航天大学出版社,.　　[6] 孙燮华.关于MATLAB符号工具箱的若干问题[J].计算机应用,):14-16.　　作者简介：　　席伟（1978-），男，东北大学硕士毕业，沈阳化工大学计算机学院讲师；主要研究方向是模式识别、微分方程应用等。 [浏览次数： ] 评论内容只代表网友观点，与本站立场无关！微分方程方向场MATLAB仿真工具箱设计.doc_百度文库 两大类热门资源免费畅读 续费一年阅读会员，立省24元！ 您可以上传图片描述问题 联系电话： 请填写真实有效的信息，以便工作人员联系您，我们为您严格保密。 微分方程方向场MATLAB仿真工具箱设计.doc ||暂无简介 北京龙源网通电子商务有限公司| 总评分0.0| 试读已结束，如果需要继续阅读，敬请购买 注：购买后，该文档仅支持在线阅读 定制HR最喜欢的简历 你可能喜欢Mathematica 与常微分方程---- 方向场和积分曲线_百度文库 两大类热门资源免费畅读 续费一年阅读会员，立省24元！ Mathematica 与常微分方程---- 方向场和积分曲线 上传于||文档简介 &&数​学​给​我​们​的​印​象​是​,​沿​定​义​→​公​理​→​定​理​→​推​论​→​证​明​这​么​一​条​演​绎​道​路​进​行​的​、​一​个​十​分​严​格​的​数​学​推​理​王​国​和​一​个​充​满​美​感​的​抽​象​世​界​。​然​而​,​我​们​却​不​知​道​,​也​许​也​没​有​想​过​,​这​些​如​此​严​密​、​完​整​、​美​妙​的​结​论​是​怎​么​来​的​？​数​学​家​是​通​过​什​么​样​的​方​式​发​现​它​们​的​？​我​们​从​这​些​可​爱​结​论​本​身​看​不​到​数​学​家​发​现​它​们​的​艰​辛​,​也​体​会​不​到​数​学​家​在​发​现​它​们​之​后​的​一​种​喜​悦​。 阅读已结束，如果下载本文需要使用5下载券 想免费下载本文？ 定制HR最喜欢的简历 下载文档到电脑，查找使用更方便 还剩3页未读，继续阅读 定制HR最喜欢的简历 你可能喜欢[微分方程模型]微分方程模型 微分方程模型 第三章 微分方程模型当我们描述实际对象的某些特性随时间（或空间）而演变的过程、分析它的变化规律、预 测它的未来性态，研究它的控制手段时，通常要建立对象的动态模型.建模时首先要根据建模目 的和对问题的具体分析作出简化假设，然后按照对象内在的或可以类比的其他对象的规律列出 微分方程，求出方程的解并将结果翻译回实际对象，就可以进行描述、分析、预测或控制了.
3.1 经济增长模型本节的模型将首先建立产值与资金、劳动力之间的关系，然后研究资金与劳动力的最佳分 配，设投资效益最大，最后讨论如何调节资金与劳动力的增长率，使劳动生产率得到有效增长。3.1.1.道格拉斯（Douglas）生产函数 道格拉斯（ 道格拉斯 ） 用 Q(t ) , K (t ), L(t ) 分别表示某一地区或部门在时刻 t 的产值、资金和劳动力，它们的关系可以一般地记作Q(t ) = F ( K (t ), L(t ))其中 F 为待定函数。对于固定的时刻 t，上述关系可写作（1） （2） （3）Q = F ( K , L)为寻求 F 的函数形式，引入记号Z = Q / L, y = K / LZ 是每个劳动力的产量， y 是每个劳动力的投资，如下的假设是合理的： Z 随着 y 的增加而增长，但增长速度递减。进而简化地把这个假设表示为Z = cg ( y ), g ( y ) = y a0 & a &1（4）显然函数 g（y）满足上面的假设，常数 c&0 可看成技术的作用。由（3） （4）即可得到(2） 式中 F 的具体形式为 由（5）式容易知道 Q 有如下性质Q =cKαL1?α,0 & α & 12（5）记 Qk =?Q ?Q , Q 表示单位资金创造的产值； Q , 表示单位劳动力 K L ?K ?LKQ QKQ Q ?Q ?Q , & 0, ? 2 , ? 2 & 0 ?K ?L ?K ?L2（6）创造的产值，则从（5）式可得= α,LQ QL= 1 ? α , K Q + L Q = Q （7）K L（7）式可解释为： a 是资金在产值中占有的份额，1- a 是劳动力在产值中占有的 份额。于是 a 的大小直接反应了资金、劳动力二者对于创造产值的轻重关系。 (5）式是经济学中著名的 Cobb-Douglas 生产函数，它经受了资本主义社会 一些实际数据的检验，更一般形式的生产函数表为Q=c K 3.1.2 资金与劳动力的最佳分配αLβ,.0 & α , β & 1（8）这里将根据生产函数（5）式讨论，怎样分配资金和劳动力，使生产创造的效益最大。假定 资金来自贷款，利率为 r ，每个劳动力需付工资 w，于是当资金 K 、劳动力 L 产生产值 Q 时， 得到的效益为 (9) S = Q ? rK ? ω L, 0 & α , β & 1问题化为求资金与劳动力的分配比例 K / L （即每个劳动力占有的资金） ，使效益 S 最大。 这个模型用微分法即可解得Q QK L=rω公式（10）再利用（7）式，有K α ω = L 1?α r（11）这就是资金与劳动力的最佳分配。从（11）式可以看出，当 a, w 变大、 r 变小时， 分配比例 K / L 变大，这是符合常识的。 3.劳动生产率增长的条件 常用的衡量经济增长的指标，一是总产值 Q(t ) ，二是每个劳动力的产值Z (t ) = Q(t ) / L(t ) ，这个模型讨论 K (t ) ， L(t ) 满足什么条件才能使 Q(t ) ， Z (t ) 保持增长。首先需要对资金和劳动力的增加作出合理的简化假设： 1） 投资增长率与产值成正比，比例系数 λ&0，即用一定比例扩大再生产； 2） 劳动力的相对增长率为常数 ?，? 可以是负数，表示劳动力减少。 这两个条件的数学表达式分别为dK = λQ, λ & 0 公式（12） dt dL = ?L 公式（13） dt方程（13）的解是*L(t ) =Le0?t公式（14）将（4）（.5）带入（12）式得 ，α dK 公式（15） = cLλ y dt dK dy = L + ?Ly (16) dt dt比较（15）（16）得到关于 y（t）的方程 ，dy + ?y = cλ dtyα公式（17）这是著名的 Bernoulli 方程，它的解是? ? cλ ? Y(t)= ? ?1 ? 1 ? ? ?? ? ?1） Q（t）增长，即K e K0 . 0? (1?α ) ?t? ? ? ? ?}1 1?α(18)一下根据（3.2.18）式研究 Q（t） ，z（t）保持增长的条件。α dQ ， & 0,由Q = cL y ，由 Q=cLy 及（13）（17）式可算得 dt α ?1 dy α 2α ?1 1?α dQ = cLα y + c?L y = cL y ?cαλ + ? (1 ? α )c?L y ? 公式（19） ? ? ? ? dt dt将其中的 y 以（18）式代入，可知条dQ & 0 等价于 dt1? ?K e K0 . 0? (1?α ) ?t&1 公式（20） 1?α因为上式右端大于 1，所以当 ?≥0（即劳动力不减少）时（20）式恒成立；而当 ? &0 时， （20）式成立的条件是? 1 t& ln ?(1 ? α ) 1 ? ? 1?α ? ? ?K K.00? ? 公式（21） ? ? ?说明如果劳动力减少，Q（t）只能在有限时间内保持增长。但应注意，若（21）式中 的 (1 ? α ) 1 ? ?K K.0≥1，则不存在这样的增长时段。0? dz dy & 0, 由z = c y 知相当于 & 0 由方程（17）知，当 ?≤0 时该条件 dt dt dy 恒成立；而当 ?&0 时由（18）式可得 & 0 等价于 dt2）z（t）怎长，即1? ?K e K0 . 0? (1?α ) ?t显然，此式成立的条件为 ?K K.0、&1，即0?&K K.0公式（23）0这个条件的含义是，劳动力增长率小于初始投资增长率。 评注 Douglas 生产函数是计量经济学中重要的数学模型，本节给出它的一 种简洁的建模过程，在此基础上讨论的资金与劳动力的最佳分配，是一个静态模 型。而利用微分方程研究的劳动生产率增长的条件，是一个动态模型，虽然它的 导过程稍繁，但其结果却相当简明，并且可以给出合理的解释。3.2 捕鱼业的持续收获可持续发展是一项基本国策，对于像渔业，林业这样的再生资源，一定要注意适度开发， 不能为一时的高产去“涸泽而渔” ，应该在持续稳产的前提下追求产量或效益的最优化。 考察一个鱼肠，其中的鱼量在天然环境下按一定规律增长，如果捕捞量恰好等于增长量， 那么渔场鱼量将保持不变，这个捕捞量就可以持续下去。本节要建立在捕捞情况下渔场鱼量遵 从的方程，分析鱼量稳定的条件，并且在稳定的前提下讨论如何控制捕捞使持续产量或经济效 益达到最大，最后研究所谓捕捞过度的问题。 产量模型 记时刻 t 渔场中鱼量为 x(t)，关于 x(t)自然增长和人工捕捞作如下假设：1.在无捕捞条件下，x(t)的自然增长服从 Logistic 规律，即& x(t ) = f ( x) = rx(1 ?x ) N(1)R 是故有增长率，N 是环境容许的最大鱼量，用 f(x)表示单位时间的增长率。 2.单位时间的捕捞量（即产量）与渔场鱼量 x(t)成正比，比例常数 E 表示单位时间捕捞率， 又称捕捞强度 捕捞强度，可以用比如捕鱼网眼的大小或出海渔船数量来控制其大小。于是单位时间的捕 捕捞强度 捞量为 (2) h( x ) = Ex 根据以上假设并记F ( x ) = f ( x ) ? h( x )得到捕捞情况下渔场鱼量满足的方程& x(t ) = F ( x) = rx(1 ?x ) ? Ex N(3)我们并不需要解方程(3)以得到 x(t)的动态变化过程，只希望知道渔场的稳定鱼量和保持稳 定的条件，即时间 t 足够长以后渔场鱼量 x(t)的趋向，并由此确定最大持续产量。为此可以直接 求方程(3)的平衡点并分析其稳定性。 令F ( x) = rx(1 ? x / N ) ? Ex = 0得到两个平衡点x0 = N (1 ?不难算出E )，x1 = 0 r(4)F ′( x0 ) = E ? r , F ′( x1 ) = r ? E所以若 (5) 有 F ( x0 ) & 0, F ( x1 ) & 0 ,故 x0 点稳定， x1 点不稳定(判断平衡点稳定性的准则见 6.6 节);若E&rE & r ,则结果相反。 E 是捕获率， r 是最大增长率，上述分析表明只要捕捞适度( E & r )，就可以使渔场鱼量稳 定在 x0 ,从而获得持续产量 h( x0 ) = Ex0 ;而当捕捞过度时( E & r ),渔场鱼量将趋向 x1 =0,当然谈不上获得持续产量了。 进一步讨论渔场鱼量稳定在 x0 的前提下，如何控制捕捞强度 E 使持续产量最大的问题。用 图解法可以非常简单地得到结果。 根据(1),(2)式作抛物线 y = f ( x ) 和直线 y = h( x ) = E ( x) ,如图 1.注意到 y = f ( x ) 在原点 的切线为 y = rx ,所以在条件(5)下 y = Ex 必与 y = f ( x ) 有交点 p , p 的横坐标就是稳定平衡点x0 .yy = rxhm hP*y = h ( x ) = ExP y = f (x )O* x 0 = N / 2 x0Nx图 3-1 最大持续产量的图解法 根据假设 2，p 点的纵坐标 h 为稳定条件下单位时间的持续产量，有图 3-1 立刻知道，当 y = Ex 与 y = f ( x ) 在抛物线顶点 p*相交时可获得最大持续产量，此时的稳定平衡点为x0 = N / 2*(6) (7)而单位时间的最大持续产量为hm = rN / 4而由(4)式不难算出保持渔场鱼量稳定在 x0 的捕捞率为(8) 综上所述，产量模型的结论是将捕捞率控制在故有增长率 r 的一半，更容易些，可以说使渔 场鱼量保持在最大量 N 的一半时，能够获得最大的持续产量。 从经济角度看不应追求产量最大， 而应考虑效益用从捕捞所得的收入中扣除开支 效益模型 后的利润来衡量，并且简单地假设：鱼的销售单价为常数 p ，单位捕捞率(如每条出海渔船)的 费用为常数 c,那么单位时间的收入是 T 和支出 S 分别为 单位时间的利润为E? = r / 2T = ph( x) = pEx, S = cER = T ? S = pEx ? cE(9) (10) (11)在稳定条件 x = x0 下，以(4)代入(10)式得R( E ) = T ( E ) ? S ( E ) = pNE (1 ? E / r ) ? cE用微分法容易求出使利润 R(E)达到最大的捕捞强度为 (12) ER = (r / 2)(1 ? c / p ( N )) 将 ER 代入(4 )式可得最大利润下的鱼场稳定鱼量 xR 以及单位时间的持续产量 hR 为 (13) XR = N / 2+ c / 2pxR rN c2 hR = rx R (1 ? ) = (1 ? 2 2 ) N 4 p N(14)将(12)~(14)式与产量模型中的(6)~(8)式比较可以看出， 在最大效益原则下捕捞率和持续产 量均有所减少，而渔场应保持稳定鱼量有所增加，并且减少或增加的部分随着捕捞成本 c 的增 长而变大，随着销售价格 p 的增长而变小。请读者分析这些结果是符合实际情况的. 捕捞过度 上面效益模型是以计划捕捞(或称为封闭式捕捞)为基础的，即渔场有单独的经营者 有计划地捕捞，可以追求最大的利润。如果渔场向众多的盲目的经营者开放，比如在公海上无规则的捕捞，那么即使只有微薄的利润，经营者也会蜂拥而去，这种情况称为盲目捕捞(或开放 式捕捞)。这种捕捞方式将导致捕捞过度，下面讨论这个模型。 (11)式给出了利润与捕捞强度的关系 R( E ) ,令 R( E ) = 0 的解为 ES ,可得ES = r (1 ? c / ( pN ))当然要减少强度。所以 ES 是盲目捕捞下的临界强度。（15）当 E & ES 时利润 R( E ) &0,盲目的经营者们会加大捕捞强度；若 E & ES 利润 R( E ) &0,他们T ( E ), S ( E )pNES ( E ) = cES ( E ) = cE T (E )OE s1E* = r /2E s2rE图 3-2 盲目捕捞强度的图解法 ES 也有图解法确定。在图 3-2 中以 E 横坐标按(11)式画出 T ( E ) 和 S ( E ) ,他们焦点横坐标即为ES (图 3-2 中的 ES1 或 ES 2 ).由(15)式或图 3-2 容易知道 ES 存在的必要条件(即 ES &0)是 (16) p&c/ N 及销售大于(相当于总量而言)成本。并且由(15)式可知，成本越低，售价越高，则 ES 越大。将(15)代入(4)式得到盲目捕捞下的渔场稳定捕鱼量为XS = c / p（17）随着价格的上升和成本的下降，X S 将迅速减少， 出现捕捞过度。 X S 完全由成本－价格比决定， 比较(12)和(15)式可知 ES = 2 ER ,即盲目捕捞强度比最大效益下捕捞强度大一倍。 从(15)式和图 3-2 还可以得到，当 c / N & p & 2c / N 时 ( ER &) ES & E??,如图 3-2 中 ES 1 ,称经济学捕捞过度：当 p & 2c / N 时 ES & E ,如图 3-2 中 ES 2 ，称生态学捕捞过度。 为了研究渔业的产量、效益即捕捞过度问题，首先在对鱼的自然增长和捕捞过度情况 评注 的合理假设下， 建立渔场鱼量的基本方程(3),并利用平衡点稳定性分析确定了保持渔场鱼量稳定 的条件。产量、效益和捕捞过度 3 个模型在稳定的前提下步步深入，数学推导过程十分简单， 却得到了在定性关系上与实际情况完全符合的结果。 如果改变对鱼的自然增长和人工捕捞的假设，模型及结果将随之变化。建模案例： 建模案例：最优捕鱼策略问题简介 生态学表明，对可再生资源的开发策略应事先可持续收获的前提下追求最大经济效益。 考虑具有 4 个年龄组：1 龄鱼，...4 龄鱼的某种鱼。该鱼类在每年后 4 个月季节性集中产卵繁殖。而 按规定，捕捞作业只允许在前 8 个月进行，每年投入的捕捞能力不变，单位时间捕捞量与各年龄组鱼群条 数的比例成为捕捞强度系数。 使用只能捕捞 3`4 龄鱼的 13mm 网眼的拉网， 其两个捕捞强度系数比为 0.42： 1.渔业上称这种方式为固定努力量捕捞。 该鱼群本身有如下数据： 该鱼群本身有如下数据：1 各年龄组鱼的自然死亡率为 0.8（1/年） ，其平均重量分别为 50.7，11.55，17.86，22.99(单位：g) （n 为产卵总量） 3 卵孵化的成活率为 1.22 × 10 /（1.22 × 10 +n） 有如下问题需要解决： ，并在次前提下得到最优捕 1） 分析如何实现可持续捕获（即每年开始捕捞时渔场中各年龄组鱼群不变） 捞量 承包的各年龄组鱼群 2） 合同鱼要求某渔业公司在 5 年合同期满后鱼群的生产能力不能受到太大的破坏，11 112 1 龄鱼和 2 龄鱼不产卵，产卵期间，平均每条 4 龄鱼产卵量为 1.109 × 105(个)，3 龄鱼为其一半；数量为 122，29.7，10.1，3.29（ × 10 条 ） ，在固定努力量的捕捞方式下，问该公司应采取作怎 样的捕捞策略，才能使总量获得最高。9基本假设 （略）α ――单位时间的自然死亡率 c ――年存活率， c = 1 ? 0.8 = 0.2 .X = ( X 1 X 2 X 3 X 4 ) ――鱼群数向量Tk ――单位时间 4 龄鱼捕捞强度系数 β ――孵化卵成活率， β = 1.22 × 1011 / (1.22 ×1011 + n )m ――4 龄鱼的平均产卵量， m 为 1.109 × 105 (个)，3 龄鱼为其一半。模型建立 这里讨论问题 2） ，即可持续捕获策略模型。以一年为一个离散化的单位时间。 记年初鱼群为 X ( t ) = X 1 ( t )X ( t + 1) = ( X 1 ( t + 1) X 2 ( t + 1) X 3 ( t + 1) X 4 ( t + 1) )显 然 ，4(X 2 ( t ) X 3 ( t ) X 4 ( t ) ) ，下半年的鱼群数为TT( i = 1, 2,3, i = 4时X ( t + 1)中还包括X ( t )中存活数。X ( t ) 指上一年由卵孵化而得到1龄鱼) ，据此可建立4 0X i ( t + 1)是X i ?1 ( t )到 年 底 存 活 下 来 的 鱼 群 数如下差分方程：?m ? X 1 ( t + 1) = β ? ( c ? k3 ) X 3 ( t ) + m ( c ? k4 ) X 4 ( t ) ? 2 ? ? X 2 ( t + 1) = cX 1 ( t ) X 3 ( t + 1) = cX 2 ( t ) X 4 ( t + 1) = ( c ? k3 ) X 3 ( t )( c ? k4 ) X 4 ( t )X ( t + 1) = PX ( t )其中因为 3 龄鱼与 4 龄鱼捕捞强度系数比为 0.42:1,故有 k3 = 0.42k 4 = 0.42k ， 写成矩阵的形式：? ?0 ? P = ?c ?0 ? ?0 ?m ? β ( c ? 0.42k ) mβ ( c ? k ) ? 2 ? 0 0 0 ? ? c c 0 ? c?k ? 0 c ? 0.42k ? 0c 0.2 = = 0.76 0.42 0.42时，不论上一年鱼群数仔细考虑矩阵 P ，当 4 龄鱼捕捞强度系数 k &目如何，下一年的鱼群将出现负数。这个结论显然是荒谬的。事实上，只要 3 龄鱼和 4 龄鱼不被 同时捕光，下一年 4 龄鱼存在存活，即鱼群数不会出现负数。 造成这种现象的原因的单位时间离散化程度不够精细。假设单位时间为一个月，定义月死亡 率为 α ，月存活率为 (1 ? α ) ，月捕捞系数为 k ，则年存活率 (1 ? α )12应为 c = 0.2 ，从而得α = 0.1255一个月实际存活率： (1 ? α ? k ) 二个月实际存活率: (1 ? α ? k ) 三个月实际存活率: (1 ? α ? k ) 八个月实际存活率: (1 ? α ? k )，考虑一年中各月鱼群数目的分布，不难得到如下分析：2 38九个月实际存活率: (1 ? α ? k ) (1 ? α )8一年后实际存活率: (1 ? α ? k ) (1 ? α )84同理可得第 i 月捕捞率： (1 ? α ? k )i ?1k , i = 1, 2,L ,8.8 4 8 4因此可得一年后 3 龄鱼实际存活数： (1 ? α ? k3 ) (1 ? α ) X 3 一年后 4 龄鱼实际存活数： (1 ? α ? k4 ) (1 ? α ) X 4 该 3 龄鱼总捕捞数：3 k3 ?1 ? (1 ? α ? k3 ) ? ?X k3 X 3 = ? 3 α + k33 k4 ?1 ? (1 ? α ? k4 ) ? ? ?X k4 X 4 = 4 α + k4∑ (1 ? α ? k )i =1 33 i =1 43i ?1该 4 龄鱼总捕捞数：∑ (1 ? α ? k )i ?1m 8 (1 ? α ? k3 ) X 3 2 m 8 该 4 龄鱼产卵总量 n4 = (1 ? α ? k4 ) X 4 2该 3 龄鱼产卵总量 n3 = 因此矩阵应当为：? 0 ?0 ? 12 ? 1? α ) 0 P = ?( 12 ?0 (1 ? α ) ? ?0 0 ?m 3 β (1 ? α ? k3 ) 2 0 0mβ (1 ? α ? k4 )(1 ? α ) (1 ? α ? k3 )43? ? ? 0 ? ? ? 0 ? 4 3 (1 ? α ) (1 ? α ? k4 ) ? ?3关于鱼群的差分方程为： X ( t + 1) = PX ( t )为实现持续捕获（1）式，必须存在稳定解：X ( t ) = PX ( t )，有由差分方程稳定性理论知其充分性为：对 P 的所有特征根 λiλi ≤ 1用 Mathematica 软件包按上述步骤得最优解如下： MAX : f = 5.94389 × 1010 g = 59438.9t 最 佳 月 捕 捞 强 度 系 数 ： 4 龄 鱼 k4 = k = 0.777697 ； 3 龄 鱼k3 = 0.42k = 0.32195；可持续最佳捕获下，渔场各年龄组鱼群数：X * = (115.222 ×109 23.0444 × 109 4.60888 ×109 2.21149 ×107 )这说明该类鱼群不论开始鱼群数目如何，经过一定时间持续捕获，总能使鱼群数目稳定下来， 且在这种稳定生长的情形下，我们可用 Mathematica 软件给出捕获量 f 与月捕获强度系数（4 龄 鱼）k 的图形（图 3-3）图 3-33.3 导弹系统的改进海军方面要求改进现有的舰对舰导弹系统.目前的电子系统能迅速测出敌舰的种类、位置以 及敌舰行驶速度和方向,且导弹自动制导系统能保证在发射后任一时刻都能对准目标.根据情报, 这种敌舰能在我军舰发射导弹后 T 分钟作出反应并摧毁导弹.现在要求改进电子导弹系统使能自 动计算出敌舰是否在有效打击范围之内. 设我舰发射导弹时位置在坐标原点,敌舰在 x 轴正向 d km 处,其行驶速度为 a km/h,方向与 x 轴夹角为 θ,导弹水平飞行线速度为 b km/h.问题的关键是求出导弹击中敌舰的时间. 设 t 时刻时导弹位置为(x(t),y(t)).那么? dx ? + ? dy ? ? ? ? ? ? dt ? ? dt ?22=b(1)易知 t 时刻敌舰位置为(d+atcosθ,atsinθ),为了保持对准目标,导弹轨迹切线方向应为dy at sin θ ? y (t ) = , dx d + at cos θ ? x(t )=(2)由(1)式和(2 式得下列微分方程 dx b = dtb 1+ ?1+?? dy ? ? ? dx ?b2? at sin θ ? y(t ) ? ? d + at cos θ ? x(t ) ? ? ? ?b2,(3)dy = dt1+ ?? dx ? ? dy ? ? ? ?2=1+ ?? d + at cosθ ? x(t ) ? ? at sin θ ? y(t ) ? ? ? ?2,(4)初始条件 x(0)=0,y(0)=0.对于给定的 a,b,d,θ 进行计算.当 x(t)满足 X(t)≥d+atcosθ, (5) 则认为已击中目标.如果 t&T,则敌舰在打击范围内,可以发射. 例 在导弹系统中设 a=90 km/h,b=450km/h,T=0.1h.现要求 d, θ 的有效范围. 解 有两个极端情形容易算出.若 θ=0,即敌舰正好背向行驶,即 x 轴正向.那么导弹直线飞行, 击中时间 t=d/(b-a)&T, 得 d=T(b-a)=36 km。若 θ=π，即迎面驶来，类似有 d=T(a+b)=54 km.一般地，应有 36&d&54. 下面我们考虑三种算法解上例。 (1)在线算法 对于测定的 d 和 θ,可用(3)式和(4)式计算出 t.比如 d=50, θ=π2,写 M 函数missile.m 为了防止分母为 0,加了一个小正数 1e-8.并且使用附加参数 a,b,d,theta 传递. function dy=missilefun(t,y,a,b,d,theta) dydx=(a*t*sin(theta)-y(2)+1e-8)/(abs(d+a*t*cos(theta)-y(1))+1e-8); dy(1)=b/(1+dydx^2)^0.5; dy(2)=b/(1+dydx^(-2))^0.5; dy=dy(:); launch1. a=90;b=450;d=50;theta=pi/2;[t,y]=ode45(@missile,[0 0.1],[0 0],[],a,b,d,theta); plot(y(:,1),y(:,2)); if max(y(:,1)-d-a*t*cos(theta))&0 answer='Please Wait' end 由于在 T=0.1h 内,(5)式不成立,所以敌舰不在有效打击范围,应等近一些再发射. answer = Please Wait 图形见图 3-4图 3-4 在线算法 (2)离线算法 首先对于所有可能的 d 和 θ,计算击中所需时间,从而对不同 θ,得 d 的临界值.具体 应用时直接查表判断.编写 m 脚本文件 launch2. %如果*在临界曲线下等待 %否则发射 answer='如果*在临界曲线下等待否则发射' a=90;b=450; i=1; for d=54:-1:36 for theta=0:0.1:pi [t,y]=ode45(@missile,[0 0.1],[0 0],[],a,b,d,theta); if max(y(:,1)-d-a*t*cos(theta))&0 range(i,:)=[d ,theta]; i=i+1; plot(range(:,1),range(:,2)); xlabel('d');ylabel('theta'); hold on d=50;theta=pi/2; plot(d,theta,'*') 运行得临界曲线.见图 3-5answer = 如果 * 在临界曲线下，等待！否则，发射！图 3-5 离线算法 图 3-5 中,曲线上方为打击范围.由于 θ=1.57,d=50 在曲线下方,这样即可知不在打击范围内. (3)计算机模拟 一个较基本但形象的方法.对于任意选定的参数 a,b,d, θ,T,下面的 M 函数 提供一个导弹追击敌舰演示工具.其中使用了 MATLAB 动画制作指令 getframe 和动画播放指令 movie. function m=launch3(a,b,d,theta,T) [t,y]=ode45(@missile,[0,T],[0 0],[],a,b,d,theta); x=[d+a*t*cos(theta),a*t*sin(theta)]; n=length(t);j=n; for i=1:n; plot(x(i,1),x(i,2),'o',y(i,1),y(i,2),'r.'); axis([0 max(x(:,1)) 0 max(x(:,2))]); m(i)= if(y(i,1)&x(i,1)) j=i; brmovie(m); legend('敌舰'，'导弹',2); if j& plot(y(j,1),y(j,2),'rh','markersize',18); hold off title(['导弹将在第',num2str(t(j)),'小时击中敌舰'); else title(['导弹在',num2str(T),'小时内不能击中敌舰']); end 对于敌舰速度 a=90 km/h,导弹速度 b=450 km/h,距离 d=30 Km,敌舰行驶角度 θ=0.3π,反应 时间 T=0.1h,在指令窗口执行 && launch3(90,450,30,0.3*pi,0.1)得图 3-6 可见导弹约在 t=0.08 时击中敌舰,位置约在(34,5.5).图 3-6 模拟算法 应该说,三种算法各有千秋.在线算法灵活,容易调整参数和模型,但速度慢.离线算法 事先计算好,实时使用查询方式,不需计算,速度极快.模拟算法比较直观、生动3.4 传染病模型随着卫生设施的改善、医疗水平的提高以及人类文明的不断发展，诸如霍乱、天花等曾经 肆虐全球的传染性疾病已经得到有效的控制.但是一些新的、不断变异着的传染病毒却悄悄向人 类袭来.20 世纪 80 年代十分险恶的艾滋病毒开始肆虐全球，至今仍在蔓延；2003 年春来历不 明的 SARS 病毒突袭人间，给人们的生命财产带来极大的危害.长期以来，建立传染病的数学模 型来描述传染病的传播过程，分析受感染人数的变化规律，探索制止传染病蔓延的手段等，一 直是各国有关专家和官员关注的课题. 不同类型传染病的传播过程有其各自不同的特点，弄清这些特点需要相当多的病理知识， 这里不可能从医学的角度一一分析各种传染病的传播，而只是按照一般的传播机理建立几种模型[11, 28, 39 ].模型 1 在这个最简单的模型中，设时刻 t 的病人人数 x(t)是连续、可微函数，并且每天每 个病人有效接触（足以使人致病的接触）的人数为常数 λ，考察 t 到 t+?t 病人人数的增加，就 有 x(t+?t)-x(t)=λx(t)?t 再设 t=0 时有 0 个病人，即得微分方程xdx = λ x, x(0) = x0 dt方程(1)的解为（1）x (t ) = x0 e λt（2）结果表明，随着 t 的增加，病人人数 x (t ) 无限增长，这显然是不符合实际的. 建模失败的原因在于：在病人有效接触的人群中，有健康人也有病人，而其中只有健康人 才可以被传染为病人，所以在改进的模型中必须区别这两种人. 模型 2（SI 模型） 假设条件为 （ 模型） 1.在疾病传播期内所考察地区的总人数 N 不变，既不考虑生死，也不考虑迁移.人群分为易 和已感染者 （Infective） 两类 （取两个词的第 1 个字母， 称之为 SI 模型） ， 感染者 （Susceptible） 以下简称健康者和病人.时刻 t 这两类人在总人数中所占的比例分别记为 s(t)和 i(t). 2.每个病人每天有效接触的平均人数是常数 λ，λ 称为日接触率.当病人与健康者有效接触 时，使健康者受感染变为病人. 根据假设,每个病人每天可使 λs(t)个健康者变为病人， 因为病人数为 Ni(t)， 所以每天共有 λNs(t)i(t) 个健康者被感染，于是 λNsi 就是病人数 Ni 的增加率，即有 N 又因为 S(t)+i(t)=1 再记初始时刻（t=0）病人的比例为 （4）di =λNsi dt（3）io，则 （5）di = λi（1-i） ，i（0）= io dt方程（5）是 1.5 节中出现过的 Logistic 模型.它的解为 i(t)=1+ (i(t)和1 1i? 1)e? λt(6)0di 的图形如图 3-7 和图 3-8 所示. dt图 3-7 SI 模型的 i~t 曲线图 3-8SI 模型的di dti 曲线由（5）（6）式及图 3-7 可知，第一，当 i=1/2 时 ，di ，这个时刻为 di 达到最大值 ( ) dt dtmt =λm?1ln（1i? 1）（7）0这时病人增加的最快，可以认为是医院的门诊量最大的一天，预示着传染病高潮的到来，是医 疗卫生部门关注的时刻. 水平越高.所以改善保健设施、 提高卫生水平可以推迟传染病高潮的到来.第二， t → ∞时 i → 1， 当 即所有人终将被传染，全变为病人，这显然不符合实际情况.其原因是模型中没有考虑到病人可 以治愈，人群中的健康者只能变为病人，病人不会再变成健康者. 为了修正上述结果必须重新考虑模型的假设，下面两个模型中我们讨论病人可以治愈的情 况. 模型） 模型 3（SIS 模型） 有些传染病如伤风、痢疾等愈后免疫力很低，可以假定无免疫性， （ 于是病人被治愈后变成健康者，健康者还可以被感染再变成病人，所以这个模型称 SIS 模型. SIS 模型的假设条件 1，2 与 SI 模型相同，增加的条件为 3.每天被治愈的病人数占病人总数的比例为常数 ? ，称为日治愈率.病人治愈后成为仍可被 感染的健康者.显然 1/ ? 是这种传染病的平均传染期. 不难看出，考虑到假设 3， SI 模型的（3）式应修正为tm与 λ 成反比，因为日接触率 λ 表示该地区的卫生水平，λ 越小卫生N（4）式不变，于是（5）式应改为di = λ Nsi ? ? Ni dt（9）（8）di = λi (1 ? i ) ? ? i, i (0) = i0 dt我们不去求解方程（9） （虽然它的解可以解析地表出） ，而是通过图形分析 i（t）的变化规 律.定义 (10) σ =λ/? 注意到 λ 和 1/? 的含义，可知 σ 是整个传染期内每个病人有效接触的平均人数，称为接触数. 利用 σ，方程（9）可以改写作di 1 = ? λi[i ? (1 ? )] （11） dt σ di 由方程 （11） 容易先画出 ~i 的图形 （图 3-9， 3-11） 再画出 i~t 的图形 图 ， （图 3-10， 3-12） 图 . dt图 3-9 SIS 模型的di ~i 曲线（σ&1） dt图 3-10SIS 模型的 i~t 曲线（σ&1）图 3-11 SIS 模型的di ~i 曲线（σ≤1） 图 3-12 SIS 模型的 i~t 曲线（σ≤1） dt 不难看出，接触数 σ = 1 是一个阀值.当 σ & 1 时 i (t ) 的增减性取决于 io 的大小（见图 3-7） ， 1但其极限值 i (∞ ) = 1 ?σ随 σ 的增加而增加（试从 σ 的含义给以解释） ；当 σ ≤1 时病人比例i（t）越来越小，最终趋于零，这是由于传染期内经有效接触从而使健康者变成的病人数不超过 原来病人数的缘故. SI 模型可视为本模型的特例，请读者考虑它相当于本模型中 ? 或 σ 取何值的情况. 模型） 模型 4（SIR 模型） 大多数传染病如天花、流感、肝炎、麻疹等治愈后均有很强的免疫力， （ 所以病愈的人既非健康者（易感染者） ，也非病人（已感染者）,他们已经退出传染系统.这种情 况比较复杂，下面将详细分析建模过程. 模型假设 1.总人数 N 不变.人群分为健康者、病人和病愈免疫的移出者（Removed）三类，称 SIR 模 型.三类人在总人数 N 中占的比例分别记作 s (t ) ， i (t ) 和 r (t ) . ，传染期接触数为 σ = λ / ? . 2.病人的日接触率为 λ，日治愈率为 ? （与 SI 模型相同） 模型构成 由假设 1 显然有 (12) s (t ) + i (t ) + r (t ) = 1 根据条件 2 方程（8）仍成立.对于病愈免疫的移出者而言应有Ndr = ? Ni dt（13）再记初始时刻的健康者和病人的比例分别是s (s00&0)和i （i00&0） （不妨设移出者的初始值r0=0） ，则由（8）（12）（13）式，SIR 模型的方程可以写作 ， ，? di ? dt = λsi ? ?i, i (0) = i0 ? ds ? = ?λsi, s (0) = s0 ? dt（14）方程（14）无法求出 s（t）和 i（t）的解析解，我们先做数值计算. 数值计算 在方程（14）中设 λ =1， ? =0.3， i (0) = 0.02, s (0) = 0.98 ，用 MATLAB 软 件编程。 输出的简明计算结果列入表 1， i (t ), s (t ) 的图形见图 3-13，图 3-14 是 i ~ s 的图形,称为相 轨线,初值 i (0) = 0.02, s (0) = 0.98 相当于图 3-14 中的p0点,随着 t 的增加, ( s, i ) 沿轨线自右向左运动.由表 1、图 3-13、图 3-14 可以看出， i (t ) 由初值增长至约 t =7 时达到最大值，然后 减少， t → ∞ ， i → 0 ； s (t ) 则单调减少， t → ∞ ， s → 0.0398 .t i(t) s(t) t i(t) s(t)0 0.7 0.7 9 31 0.5 10 0.5表 3-1 i（t） ，s（t）的数值计算结果 2 3 4 5 6 0.5 0.5 0.9 0.87 0.58 0.915 20 25 30 0.3 0.7 0.3 0.835 40 45 0.1 0 0.9 0.0398图 3-13 i（t） ，s（t）图形 图 3-14 i~s 图形（相轨线） 为了分析 i（t） ，s（t）的一般规律，需要进行相轨线分析. 相轨线分析 我们在数值计算和图形观察的基础上， 利用相轨线讨论解 i t） （ ， s （t） 的性质. s~i 平面称为相平面，相轨线在相平面上的定义域（s，i）为 D={（s，i）|s≥0，i≥0，s+i≤1} （15） 在方程（14）中消去 dt 并注意到 σ 的定义（10） ，可得di 1 = ? 1, i | = i0 s=s0 ds σs容易求出方程（16）的解为（16）i = ( s 0 + i 0) ? s +1σlnss(17)0在定义域 D 内， （17）式表示的曲线即为相轨线，如图 3-15 所示.其中表示了随着时间 t 的增加 s（t）和 i（t）的变化趋势. 下面根据（14）（17）式和图 3-15 分析 s（t）和 r（t）的变化情况（t→∞时它们的极限值 ， 分别记作s ， 和r i∞ ∞∞）.001.不论初始条件 其证明如下.s ，i如何，病人终将消失，即i首先，由（14） ，∞=0（18）ds dr ≤0，而 s（t）≥0，故 s∞ 存在；由（13） ， ≥0，而 r（t）≤1，故 r ∞ dt dt存在；再由（12）知 其次，若 盾.i∞存在.r∞=ε&0，则由（12） ，对于充分大的 t 有dr ε &? ,这将导致 r ∞ =∞,与 r ∞ 存在相矛 dt 2从图形上看,不论相轨线从p 或从 p 点出发,它终将与 s 轴相交(t 充分大).1 2图 3-15 SIR 模型的相轨线 2.最终未被感的健康者的比例是s∞,在(17)式中令 i=0 得到,s∞是方程 (19)s +i ? s0 0∞+1σln s∞ = 0s0在(0,1/σ)内的根.在图形上 3.若s∞是相轨线与 s 轴在(0,1/σ)内交点的横坐标.s0&1/σ,则 i(t)先增加,当 s=1/σ 时,i(t)达到最大值im= s0 + i0 ?1σ∞(1 + ln σs)0(20)100然后 i(t)减小且趋于零,s(t)则单调减小至s4.若 S0 ≤ 1/ σ ,则 i (t ) 单调减小至零,p ( s , i )出发的轨线. s (t ) 单调减小至 s ,如图 3-15 中 p ( s , i,如图 3-15 中由∞200)出发的轨线. 可以看出,如果仅当病人比例 i(t)有一段增长的时期才认为传染病在蔓延,那么 1/σ 是一个阀 值,当 S0 & 1/ σ (即 σ & 1/ S0 )时传染病就会蔓延.而减小传染期接触数 σ ,即提高阈值 1/ σ ,使得S0 ≤ 1/ σ (即 σ ≤ 1/ S0 ),传染病就不会蔓延(健康者比例的初始值 s0 是一定的,通常可认为 s0 接近 1). 并且,即使 S0 & 1/ σ ,从(19),(20)式可以看出,σ 减小时,也控制了蔓延的程度.我们注意到在 σ = λ / ? 中，人们的卫生水平越高，日接触率 λ 越小；医 疗水平越高，日治愈率 ? 越大，于是 σ 越小，所以提高卫生水平和医疗水平有助于控制传染病 的蔓延. 从另一方面看, σ s = λ s ? 1/ ? 是传染期内一个病人传染的健康者的平均数,称为交换数,其 含义是一个病人被 σ s 个健康者交换.所以当 S0 ≤ 1/ σ ,即 σ S0 ≤1 时,必有 σ s ≤1.既然交换数不 超过 1,病人比例 i (t ) 绝不会增加,传染病不会蔓延. 群体免疫和预防 根据对 SIR 模型的分析,当 S0 ≤ 1/ σ 时传染病不会蔓延.所以为制止蔓 延,除了提高卫生和医疗水平,使阈值 1/ σ 变大以外,另一个途径是降低 接种使群体免疫的办法做到. 忽略病人比例的初始值s∞增加(通过作图分析),im降低,s0,这可以通过比如预防i0,有s0=1-r0.于是传染病不会蔓延的条件 S0 ≤ 1/ σ 可以表为r0 ≥ 1 ?1σ(21)这就是说,只要通过群体免疫使初始时刻的移出者比例(即免疫者比例)r0满足(21)式,就可以制止出传染病的蔓延. 这种办法生效的前提条件是免疫者要均匀分布在全体人口中,实际上这是很难做到的.据估 计当时印度等国天花传染病的接触数 σ=5,由(21)式至少要有 80G的人接受免疫才行.据世界卫 生组织报告,即使花费大量资金提高 0 ,也因很难做到免疫者的均匀分布,使得天花直到 1977 年r才在全世界根除.而有些传染病的 σ 更高,根除就更加困难. 数值验证与估计 根据上面的分析,制止传染病蔓延有两种手段,一是提高卫生水平和医 疗水平,即降低日接触率 λ,提高日治愈率 ?,二是群体免疫,即提高移出者比例的初值 0 (相当于降r低健康者比例的初值 染的健康者的比例s0).下面做一点数值计算,验证并估计这两种办法的效果.不妨用最终未感s∞和病人比例的最大值im,作为传染蔓延程度的度量指标. ,用(20)式计算给定不同的 λ,?,s ,i00,用(19)式计算sm∞im(当s0&1/σ),结果列入表 2.表2 Λ 1.0 0.6 0.5 0.4 1.0 0.6 0.5 0.4 ? 0.3 0.3 0.5 0.5 0.3 0.3 0.5 0.50s1/σ∞和i的计算结果s0i0s∞im0.3 0.5 1.0 1.25 0.3 0.5 1.0 1.250.98 0.98 0.98 0.98 0.70 0.70 0.70 0.70∞0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02m0.5 0.2 0.6 0.50.5 0.0 0.8 0.00s ,降低 λ,提高 ?,会使 s 变大, i 变小;对于一定的 λ,?,降低 s (即提 高 r ),也会使 s 变大(但是 s ≤1/σ 时 s 反而小了,你能解释吗?), i 变小.当然 s ≤1/σ 时 i 始 终等于 i ,即传染病不会蔓延.0可以看出,对于一定的∞0∞m0m0我们看到在 SIR 模型中 σ=λ/? 是一个重要参数,实际上 λ,? 很难估计,而当一次传染病结束以 后,可以获得s 和s0∞,在(19)式中略去很小的σ=ln s0 ? ln s∞i0,即有 (22)s ?s0∞当同样的传染病到来时,如果估计 λ,? 没有多大变化,那么就可以用上面得到的 σ 分析这次传染病 的蔓延过程. 模型验证 上世纪初在印度孟买发生的一次瘟疫中几乎所有病人都死亡了.死亡相当于 移出传染系统,有关部门记录了每天移出者的人数,即有了 数据对 SIR 模型做了验证.dr 的实际数据,Kermack 等人用这组 dt首先,由方程(12),(14)可以得到s (t ) = s 0 e?σr ( t )(23) (24)dr ?σr = ? (1 ? r ? s0 e ) dt当 r≤1/σ 时,取(24)式右端e2?σrTaylor 展开的前 3 项,在初始值r (t ) =其中1α2=(s 0σ ?1) + 2 s i σ2 0 0sσ0? ? α?t ?? ?( s0 σ ? 1) + αth? 2 ? φ ?? ? ?? ?2r0=0 下的解为 (25)thφ =s σ ? 1 .从(25)式容易算出0dr = dt被传染比例的估计 的初始值α2 s0σ2α2ch? 2 ? α?t? ?φ ? ? 2 ? ?(26)在一次传染病的传播过程中，被传染人数的比例是健康者人数比例s 与s0∞之差，记作 x，即x = S0 ? S∞? x? ln?1 ? ? ≈ 0 σ ? s0 ? ? ? 1（27）当i0很小，s0接近于 1 时，由（19）式可得x+（28）取对数函数 Taylor 展开的前两项有? 1 x x ?1 ? ? ? s σ 2 2σ 0 s0 ?记? ?≈0 ? ?（29）s0=1σ+δ，δ 可视为该地区人口比例超过阈值 1/σ 的部分.当 δ&&1/σ 时（29）式给出1? ? x ≈ 2 s0 σ ? s0 ? ? ≈ 2δ σ? ?（30）这个结果表明，被传染人数比例约为 δ 的 2 倍.对一种传染病，当该地区的卫生和医疗水平 不变，即 σ 不变时，这个比例就不会改变.而当阈值 1/σ 提高时，δ 减小，于是这个比例就会降 低. 评注 本节介绍的传染病模型从几个方面很好的体现了模型的改进、建模的目的性，以 及方法的配合. 第一，最初建立的模型 1 基本上不能用，修改假设后得到的模型 2 虽有所改进，但仍不符 合实际.进一步修改假设，并针对不同情况建立的模型 3，4 才是比较成功的. 第二，模型 3，4 的可取之处在于它们比较全面地达到了建模的目的，即描述传播过程、分 析感染人数的变化规律，预测传染病高潮的到来时刻，度量传染病蔓延的程度并探索制止蔓延 的手段. 第三，对于比较复杂的模型 4，采用了数值计算，图形观察与理论分析相结合的方法，先 有感性认识（表 3-1，图 3-13，图 3-14） ，再用相轨线作理论分析，最后进行数值验证和估算. 可以看作计算机技术与建模方法的巧妙配合.3.5 微分方程稳定性理论微小的干扰因素对于物质系统运动的影响，对不同的运动是不一样的。对于一些运动，这 种影响并不显著，因而受干扰的运动与不受干扰的运动相差很小；反之，对于另外一些运动， 干扰的影响就可能很显著，以至于干扰力无论多么小，受干扰的运动与不受干扰的运动随着时 间的推移而可能相差很大。这就提出了干扰对运动的影响程度的研究，从而建立了运动的稳定 性理论。李雅普诺夫给运动的稳定性以精确的数学定义并系统地解决了运动稳定性问题，创造 了解决稳定性问题的两种方法。第一种方法是寻求扰动微分方程组的通解或特解，以级数形式 将他们表达出来，最简单的级数形式是任意常数的正整数幂级数，而在该基础上研究稳定性问 题。第二种方法则不需考虑扰动微分方程组的级数形式，而仅借助于一个称为李雅普诺夫函数 的辅助函数与扰动微分方程组所计算出来的全导数的符号性质来直接推断稳定性问题，亦称李 雅普诺夫直接法。以下简单介绍利用第二种方法来判断方程组的稳定性问题。3.5.1 平衡点假设我们考虑动力系统有下列方程组描述， 其中 y = ( y1 , y2 , L , yn )Tdy （1） = g (t, y ) dt ∈ R n , g (t , y ) 在区域 G* = ( t , y ) y ∈ D* ? R n , t ≥ t0 内连续且{}满足解的唯一性条件。 如果对于某个常数点 y0 ∈ D* ,有 g (t , y0 ) = 0 ,则称点 y0 是系统（1）的平衡点 平衡位置， 平衡点或平衡位置 平衡点 平衡位置， 且称 y = y0 是系统（1）的平凡解 平凡解。 平凡解 设在 t = t0 时刻过初始点 y0 的解为y = ? (t , t0 , y0 ) = ? (t )（2）则我们关心的是：在时刻 t = t0 过 y0 邻域内任一点 y1 的解 y (t , t0 , y1 ) 与解（2）当 t 增加时是否 还 在 解 (2) 的 邻 域 内 ； 换 句 话 说 ， 即 是 否 有 ； 任 给 ε & 0, ?δ & 0 , 当 y1 ? y0 & δ 时 ， 有 | y (t , t0 , y1 ) ? ? (t ) & ε (t ≥ t0 ) 。如果对一切 t ≥ t0 ,上面的等式均成立，则称解 y = ? (t ) 是稳定 稳定 不稳定的. 的；否则称其是不稳定的 不稳定的 为了讨论方便，先作变换，令 则有x = y ? ? (t )dx dy d? (t ) = ? = g (t , y ) ? g (t , ? (t )) = g[t , x + ? (t )] ? g[t , ? (t )] = f (t , x ) dt dt dt 显然有 f ( t ,0 ) ≡ 0 ，故系统（1）变为 dx = f (t , x) dt（3）于是可知 y = ? (t ) 是系统（1）的解对应于（3）的解。这样，研究方程组（1）的任一特 y= ? (t ) 的稳定性问题可转化为研究系统（3）的平凡解 x = 0 的稳定性问题。3.5.2 稳定性的定义对于系统（3）的零解的稳定性定义，通常是指李雅普诺夫下的稳定性意义下的稳定性定义。 假设系统（3）右端函数满足条件： （1）f ( t, 0) ≡ 0 ,（2） 存在某个 H ，使的在 GH =在 f (t , x ) 的假设下， x = 0 是系统（3）的平凡解且系统（3）满足解的存在唯一性、延拓、 解对参数的连续性、可微性等条件。下面给出稳定性的定义。 定义 1{(t, x )x ≤ H .t ≥ t0 , x ∈ R n 上有连续的偏导数。}?ε & 0, ?δ = δ (t 0 , ε ) & 0, 若当 x0 & δ 时，有x(t , t 0 , x0 ) & ε (?t ≥ t 0 )则称方程组（3）的零解是稳定 稳定的。 稳定 定义 2 如果方程组（3）的零解是稳定的，且有这样的 δ 0 & 0 ,使当 x0 & δ 0 时，满足初t →+∞始条件 x (t0 ) = x0 的解 x(t , t0 , x0 ) ,均有 lim x (t , t0 , x0 ) = 0 ，则称方程组（3）的零解是渐近进稳 渐近进稳 定的。 定义 3?ε 0 & 0, ?δ & 0, 若 ?xo & 0, 满足 x0 & δ ,且 ?t1 &t 0 ，有 x(t1 , t 0 , x0 ) ≥ ε 0则称方程组（3）的零解是不稳定 不稳定的。 不稳定 定义 4 若 x = 0 渐近稳定，且存在域 D0 当且仅当 x0 ∈ D0 时满足初始条件 x (t0 ) = x0 的解 则称零解 x = 0 为全局渐近稳定 全局稳定。 全局渐近稳定或简称全局稳定 全局稳定 一致稳定的。 定义 5 在定义 1 中，若 δ 与 ε 的选取无关，则称零解是一致稳定的 一致稳定的t →+∞x ( t ) 均有 lim x (t , t0 , x0 ) = 0 ，则域 D0 称为（渐近）稳定域。若稳定域为全空间，即 δ 0 = +∞ ，3.5.3 关于李雅普诺夫函数（V 函数） 关于李雅普诺夫函数（ 函数）李雅普诺夫第二方法的特点是不必求出微分方程的解，而借助于构造一个具有特是性质的 函数，结合方程组本身来讨论稳定性。此函数即李雅普诺夫函数。 考虑纯量函数 V ( x ) V ( t , x ) 满足如下条件： （1） V ( 0 ) = 0 V ( t , x ) = 0 ， （2） 在定义域(())DH = {x | x & H }(GH = {(t , x) | x ∈ DH , t ≥ t0 })上有定义、连续和有连续的一阶偏导数；如果要求全局稳定，只需要取 H=+∞即可。 下满给出定号函数的定义。 是定负函数 定负函数，则称 V ( x ) 是定负函数；定正和定负函数统称为定号函数 定号函数。 定负函数 定号函数 负函数，则称 V ( x ) 是常负函数；常正和常负函数统称为常号函数 常号函数。 负函数 常号函数 变号函数。 定义 8 称不是常号的函数为变号函数 变号函数 定义 9 定正函数； 定义 6 如果对 ? x ∈ DH \ {0} ,均有 V ( x ) & 0 ,则称 V ( x ) 是定正函数 如果 ?V ( x ) 定正函数 常正函数； 定义 7 如果对 ?x ∈ DH ,均有 V ( X ) ≥ 0 ,则称 V ( x ) 是常正函数 如果 ?V ( x ) 是常 常正函数 常?V ( t , x ) 是定正函数，则 V ( t , x ) 称定负函数。 定负函数。 定负函数定义 10设 ω (x ) 是定正函数，若 V ( t , x ) ≥ ω ( x) ，则称 V ( t , x ) 是定正函数；如果 定正函数； 定正函数 若对一切 ( t , x ) ∈ GH ,均有 V ( t , x ) ≥ 0 ，则称 V ( t , x ) 是常正函数；如果 常正函数； 常正函数?V ( t , x ) 是常正函数，则 V ( t , x ) 称常负函数。 常负函数。 常负函数{( x x V )} 中，V = V ( x x ) 是一个位于坐标面 Ox x 上方的曲面：在 x1, 2,1, 2定号函数的几何意义。仅考虑二维空间中的定正函数 V = V x1, x2 .在三维空间1 2()2 12 + x2 ≤ H 2 内它，用水平面 V = c 切割曲面 V = V x1, x2 ，当 c 是适当 与坐标面只有在原点接触（图 3-16） 小 的 正 数 时 ， 截 线 是 一 条 封 闭 曲 线 V x1, x2 = c ， 其 投 影 Ox1 x2 平 面 上 就 是 曲 线V ( x1, x2 ) = c 。由于 ( x1, x2 ) 连续可微，且 V ( 0, 0 ) = 0 ，所以当 c 取一系列足够小的正数时， V x1, x2 = c 是包围原点的封闭曲线族。且当 c1 & c2 & 0 时，闭曲线 V x1, x2 = c2 整 个包含在闭曲线 V x1, x2 = c1 之内，因为 V x1, x2 是单值的，所以这族曲线不相交。()()()()()()图 3-163.5.4 稳定性的基本定理1. 李雅普诺夫稳定性定理 考虑驻定系统性条件，并且有 f ( 0 ) = ( 0 ) ,dx （4） = f (x) dt 其中 f ( x ) 在 DH = {x | x & H } （H&0 常数）上有定义，且系统（4）的解满足存在唯一f ( x ) = ( f1 ( x1 , x2 , L, xn ) , f 2 ( x1 , x2 ,L, xn ) ,L, f n ( x1 , x2 , L, xn ) )T设纯量函数 V ( X ) 在 DH 上有定义、 连续且有一阶偏导数和 V ( 0 ) = 0 .取函数 V ( x ) 沿着系 统（4）解的全导数为n n dV ?V dx i ?V =∑ =∑ f i ( x1 , x 2 , L, x n ) = GrandV ( x) ? f ( x) = ω ( x) dt i =1 ?x i dt i i =1 ?x i（5） 显然 ω (x ) 在 DH 上有定义、连续且满足 ω ( x ) = 0 . （1） V ( x ) 定正（定负） ，定理 6 对于方程组（4）的解 x = 0 ，如果存在 V ( x ) ，满足 （2） 由（5）式定义的 ω ( x) ≤ 0 ( ≥ 0 ) ， 定理 7 对于方程组（4）的解 x = 0 ,如果存在 V ( x ) ，满足 （1） V ( x ) 定正（定负） ， （2）由（5）式定义的 ω (x ) 定负（定正） ， 则系统（4）的解 x = 0 是渐进稳定的。 2 李雅普诺夫不稳定性定理 定理 8 对于系统（4）的解 x = 0 ,如果存在 V 函数，满足 （1） V 是定号函数， 则方程组（4）的零解是稳定的。（2）由（5）式定义的 ω ( x) ≤ 0 ( ≥ 0 ) ，但在原点的任一邻域内至少有一点 x0 ，使V ( x0 ) & 0 ( & 0 )则称系统（4）的解 x = 0 是不稳定的。 定理 9 对于系统（4）的解 x = 0 ,如果存在 V 函数，满足 （1） V = λV + ω ( x ) ，其中 λ & 0 是常数， （2）ω ( x ) ≡ 0 或 ω (x ) 常号，且存在任意小的 x ，使得 Vω & 0 ，则称系统（4）的解是不稳定的。 3 非驻定系统的稳定性定理 考虑非驻定系统dx （6） = f (t , x) dt 其中 f ( t , x ) 在 I 0 × DH ( I 0 = (t0 , +∞)) 上有定义、 连续且满足解的唯一存在性条件， 对一切 t ≥ t0 ,有 f ( t, 0) ≡ 0 . 考虑纯量函数 V ( t , x ) ,使得 V ( t , x ) 在 I 0 × DH 上有定义，连续和有连续的一阶偏导数，且对一切 t ≥ t0 均有 V ( t , x ) ≡ 0 .（1） V ( t , x ) 定正（定负） ， （2） V ( t , x ) ≤ 0 ( ≥ 0 ) ， 则系统（6）的解是稳定的。定理 10 对于系统（6）的解 x = 0 ,如果存在函数 V ( t , x ) ，满足定理 11 对于系统（6）的解 x = 0 ,如果存在函数 V ( t , x ) ，满足 （1） V ( t , x ) 定正（定负） ， （3） V ( t , x ) ≤ 0 ( ≥ 0 ) ， （2） 当 x → 0 时 V ( t , x ) 一致趋于零， 则系统（6）的解是一致稳定的。3.6 一阶微分方程定性解的图示 一阶微分方程定性解的对一阶微分方程 y`=f(x,y)利用积分求得解析解的种类少之又少，而对解 y=y(x)的定性（而非 定量）描述，可以得到非常有用的信息。如他能使你考察解得单调，极值等状态。在许多情形 下，我们并不需要明确解，尤其在求解存在技术困难或在初等函数范围内，解根本不存在时。 定性的近似解就成为解决实际问题的一种有效方法。当然，定性解应为一族解。 如果一个一阶方程可写为 y`=f(x,y)，那么我们就能够确定其解 y=y(x)在任意点（x，y）的斜 率，从图像上看，给出平面上一些离散点，通过每一点，可以画出一条“短”线，斜率为 f(x,y)， 这个图称为斜率为 f(x,y)的方向场或矢量场。 例 1 绘制dy = xy 的矢量场。 dx先说明手工绘制的方法，获得一定的感性认识，之后给出一个 MATLAB 程序，名为 vector.m， 其程序编写按照手工绘制原理。 起初，要做一些简单的运算，先确定 x,y 的变化范围，如 x ∈ [?2,2], y ∈ [?1,1] ，并给出 x，y 的基本步长，对给定的节点 ( xi , y i ) ∈ [ ?2,2] * [ ?1,1] ，计算该点的切线斜率值 f ( xi , y i ) 之后， 以该点为起点，以 f ( xi , y i ) 为斜率画一条带箭头的短线段 lij 。为保证所画线段斜率确为f ( xi , yi ) ，可以利用 dyi = dxi ? f ( xi , y i ) ，其中 dxi 为短线段 lij 在 x 轴的投影长度，dyi 为 lij 在y 轴 的 投 影 长 度 ， 你 不 妨 将 所 有 dxi 取 为 同 一 长 度 ， 如 dxi =0.1 ， i = 1, 2, L , m ， 其 中x1 = ?2, x m = 2. dxi 的大小应参考步长的大小，手工计算还应该注意一下事实。（1） 是否存在关于原点及两坐标轴的对称性，若存在，则只需考虑第一象限。 （2） 对固定的 x，斜率随 y 变化而变化的趋势， （3） 对 y 考虑 的问题 向量场的绘制可由计算机来完成。程序见附录 3.1。 运行程序时， 首先要求输入 f ( x, y ) = ? ,如输入 x + y, ? x.* y, exp( x. ^ 2 + y. ^ 2) 等等。 关于 x 和y 的变换范围及取值步长，由语句 [ x, y ] = meshgrid () 给出，对不同的函数 f ，相应绘制的矢量场范围亦应该作调整。 例 2 利用已有方向场，试试定性的给出几条解曲线，应注意的是过（ xi , y i ）点的那条解曲线 在（ xi , y i ）处的切线应为短线 lij ，即解曲线必须与方向场相符合。 在画解曲线的过程中，你会遇到一些问题。其一，是否会有两条解曲线过同一个点 ( x, y ) 。我 们以下的讨论基于如下的事实，对大多数的解而言（当 f 连续因而 y 连续时） ，唯一性告诉我 们，任意两个解不会相交。其二，即使如上所论，你仍会感到困难，你会发现，从一点开始画 解曲线，到下一点，曲线是上升还是下降，而曲线的凹凸方向亦是产生绘图困难的主要原因之 一。下面就解曲线的绘制作进一步的讨论。 考虑最简单的一类人口增长模型。 若记时刻 t 时， 人口数量为 N (t ) ， 又假设人口增长速率与 t 时 刻人口数量 N (t ) 成正比，比例系数为 r0 ，则有dN = r0 N (t ) dtN(t)（r0 &0）.o 图 3-17t显然，对这一具体问题，只需讨论解在第一象限的分布，即 N & 0, t & 0 . 可以看出，当 N 增加时，斜率亦增加，故解曲线不是直线。你应猜出，解曲线类似于图 3-17 它显示了对所有解， N 随时间增加而不断加快上升而无止境的趋势。 进一步考虑由于人口增加而过分拥挤时人口增长率递减时产生的情况，给出一个改进的模型 为dN N = r0 N (1 ? ) = 斜率(k & 0) dt k 可以看出，当 N ＝０或 N ＝ k 时，斜率＝０；当 N ＜ k 时，斜率＞０；当 N ＞ k 时，斜率＜０。据此，还不能做出不相交的图解，参看图 3-18 NKt 图 3-18再利用二阶导数确定解曲线的凹凸。2N d 2N N 2 )。 = r0 N (1 ? )(1 ? 2 k k dt d 2N 因为在 中，不显含 t ，当固定 r0 与 k 时，对给定的 N 值，对一切 t ，在（ t ， N ）处有 dt 2 相同的斜率值，因此全部解曲线可由一条解曲线沿 t 轴平移得到。上述例子引导我们给出一般图示近似解的方法。 （1） 将一阶导数与二阶导数的有关信息结合起来，绘制解图像草图。 （2） 斜率中不显含自变量，解图像可由一条解曲线沿自变量轴左右平移得到，若不含 因变量，解图像可由一条解曲线沿变量轴上下平移得到。 （3） 结合草图与计算机给出向量场绘制较精确的图示解（此一步常可省略） ，有时亦可 给出一些特殊点与线上的矢量图，以帮助确定解的形式。3.7 地中海鲨鱼问题意大利生物学家 D ' Ancona Яτ谟憷嘧苋合嗷ブ圃脊叵档难芯俊４拥谝淮问澜绱笳 期间地中海各港口捕获的几种鱼类捕获量百分比的资料中，发现鲨鱼等的比例有明显增加（见 下表） ，而提供其捕食的食用鱼的百分比却明显下降。显然战争使捕鱼量下降，食用鱼-增加， 鲨鱼等也随之增加，但为何鲨鱼的比例大幅增加呢？ 16 年代 11.9 21.4 22.1 21.2 36.4 百分比 21 年代 27.3 16.0 15.9 百分比 他无法解释这个现象，于是求助于著名的意大利数学家 V .Volterra ，希望能建立一个食饵 ―捕食系统的数学模型，定量的回答这个问题。Volterra 模型食饵和捕食者在时刻 t 的数量分别记作 x1 ( t ) 和 x2 ( t ) 。因为大海中资dx1 = r1 x1 。捕食 dt 者的存在使鱼饵的增长率降低，设降低的程度与捕食者数量成正比，于是 x1 ( t ) 满足方程源丰富，可以假设如果食饵独立生存则将以增长率 r1 按指数规律增长，即有其中，比例系数 λ1 反应捕食者掠取食饵的能力。dx1 = x1 ( r1 ? λ1 x2 ) dt（1）dx2 = ?r2 x2 ，而食饵为它提供食 dt 物的作用相当于使死亡率降低，或使之增长。设这个作用与食饵数量成正比于是 x2 ( t ) 满足捕食者离开食饵无法生存，若设它独自存在时死亡率为 r2 ，即dx2 = x2 ( ?r2 + λ2 x1 ) dt（2）其中，比例系数 反映食饵对捕食者的供养能力。 方程（1） （2）是人工捕获情况下自然环境中食饵与捕食者之间的制约关系，是 Volterra 提 出的最简单的模型。可以看出这个模型没有考虑自身的阻滞作用，即未引入了 Logistic 项。 模型分析 仍然通过平衡点的稳定分析，研究 x1 ( t ) 和 x2 ( t ) 的变化规律。容易得到方程 （1）(2)的平衡点为?r r ? P0 ? 2 , 1 ? , P0 ' ( 0, 0 ) ? λ2 λ1 ?下面用分析相轨线的方法解决这个问题 为了确定相轨线，从（1）(2)中消去 dt，得x (r ? λ x ) dx1 = 1 1 1 2 dx2 x2 ( ?r2 + λ2 x1 )r2 1（3）将变量分离，再积分得到 上式可以改写成?r2 ln x1 + λ2 x1 = r1 ln x2 ? λ1 x2 + c1(xe ? λ2 x1 )( x2 r1 e ? λ1x2 ) = c（4）其中，c 是任意常数。为了研究（4）式确定的相轨线的图形，记ψ ( x2 ) = ( x2 r e? λ x1φ ( x1 ) = ( x1r e ? λ x22 1)( 5)1 2)( 6)图 3-19图 3-20利用数学分析的方法可以作出 ? 和 ψ 的图形(图 3-21)，若他们的极大值记作 ? m 和ψ m ， 则不难确定x x100 2满足? ( x10 ) = ?m , x10 =λ2r1r2(7) (8)λ1 显然仅当(4)式右端常数 c ≤ ? mψ m 时相轨线才有意义。0 0ψ ( x20 ) = ψ m , x2 0 =当 c = ? mψ m 时， x1 = x1 , x2 = x2 , 将（7） （8）与（3）式比较可知 x1 , x20(0) 正式平衡点P0 ，所以 P0 是相轨线的退化点。为 了 考 察 c ≤ ? mψ mx2 = x20，则由（4）―（8）式得 ? ( x1 ) = α 。而从图 3-21 知道，必存在 x1 ' 和 x1 '' 使时 ( c & 0 ) 轨 线 的 形 状 ， 首 先 设 c = αψ m ( 0 & α & ψ m ) 。 若 令? ( x1 ') = ? ( x1 '') = α 且 x1 ' & x10 & x1 '' 于是这条轨线应通过点。 接着， 分析区间 ( x1 ', x1 '' ) 内的任意一点0Q1 ( x1 '', x2 0 )和 Q2 x1 ', x2 0()知 ψ ( x2 ) & ψ m ， 记 ψ ( x2 ) = β ( & ψ m ) ， 从 图 3-21 知 道 ， 存 在 x1 ' （图 3-22） 。注意到 x1 是 点 ( x1 ', x1 '' ) 示的封闭曲线，同时它绝不会越出区间 平衡点 P0。 因为 ? ( x1 ) & α ， 代人 ? ( x1 )ψ ( x2 ) = αψ m 可 和 x2 '' ， 使ψ ( x2 ' ) = ψ ( x2 '') = β ， x2 ' & x2 & x2 '' 。 且 于是这条轨线又通过 Q3 ( x1 , x2 ') 和 Q4 ( x1 , x2 '' ) 点内的任一点，立即可知这条轨线是如图图 8―13 所 。这样，对于不同的 c 值 ( 0 & c & ? mψ m ) ，方程（1） （2）的解（4）式确定的轨线是一族以 为中心的封闭曲线，称为闭轨线族。当 c 由 ? mψ m 变小时闭轨线向外扩展。图 3-21P0 点和一条相轨线0图 3-22 闭轨线族及其方向0闭轨线的方向很容易确定。考察相平面上被 x1 = x1 和 x2 = x2 域两条直线分成的 4 个区dx1 dx2 、 的正负。 dt dt程 的 周 期 解 周期为 T ， 解的示意图， 其增减由方程（1） （2）可得如图 8―14 所示的结果，因而决定了 闭轨线如箭头方向。 闭轨线对应着方x1 ( t ) 、x2 ( t ) ，记图 3-23 画出了周期 性是由图 3-22 闭 轨线的方向决定的。图 3-23 可以看出，食饵 的变化比捕食者 提前了 T 。闭轨线 （周期解） 的存在说明 P0 x10 , x2 0 别在 x0 1和 x20上下振动。 我们只能用 x1 ( t )()1 4点不是 （渐进） 稳定的。x1 ( t ) 和 x2 ( t )和 x2 ( t )分在一周期 T 内的平均值作为食饵和捕食者的数量的近似度量。记这两个平均值分别为 x1.和 x2 。因为方程（2）可写作1 x x1 (t ) = ( 2 + r2 ) λ2 x2所 以 容 易 由 此 算 出 x1 ( t ) 在 T 内 的 平 均 值 （ 利 用 x2 (T ) = x2 ( 0 ) ）r 1 T ∫0 x1 ( t )dt = λ22 T r r r 0 0 类似地可得 x2 = 1 ，于是 x1 = x1 = 2 , x2 = x2 = 1 x1 =λ1λ2λ1（9）这表明食饵和捕食者在平衡点 P0的值正好代表了它们的（平均）数量。 和 λ2 ， 和 λ1 。当食饵的自然增长率 下降模型解释 （9）式表明，食饵的数量取决于捕食者方程（2）中的两个参数 r2 而捕食者的数量取决于食饵方程（1）中的两个参数 r1时捕食者的数量将减少，这就是说，在弱肉强食情况下降低弱者的繁殖率可以使强者减少。而 当捕食者掠取食饵的能力 λ1 提高时也会使捕食者减少。一方面，捕食者死亡率 r2 的下降， 或者食饵对捕食者供养能力 λ2 为 了 用 型解释本节 期间捕获量 者）比对食用 的问题，需要 下所得的结 人工捕获的 的提高，都会导致食饵的减少。 模 开头提出的战争 下降对鲨鱼 （捕食 鱼（鱼饵）更有利 在上述自然环境 果的基础上考虑 影响（图 3-24）Volterra图 3-24 捕获系数的改变对食饵、捕食者数量的影响 设表示捕获能力的系数为 e，相当于食饵的自然增长率由 r1 降为 r1 ? e 亡率由 r2 增为 r2 + e，捕食者的死P0变为 P0 '，可利用（3）式的结果得到 y1 ( t )。用 y1 ( t ) 和 y2 ( t )表示这种情况下食饵和捕食者的数量，平衡点由 和 y2 ( t ) 的平均值为y1 =r2 + eλ2, y2 =r1 ? eλ1战争期间捕获系数由 e 下降为 e1，食饵 z1 ( t )和捕食者 z2 ( t )的平均值为z1 =r2 + e1λ2, z2 =r1 ? e1λ1因为 e1 & e，所以显然有 z1 & y1 , z2 & y2平衡点又变为 P0 ''，这就是说，战争时期捕获能力的下降使食用鱼（鱼饵）数量减少，而鲨鱼（捕食者）的数量增加， Volterra 用他的模型解释了 D ' Ancona 提出的问题。 下面针对一组具体的数据用 Matlab 软件进行计算。设 鱼 饵 和 捕 食 者 的 初 始 数 量 分 别 为 x1 ( 0 ) = x10 , x2 ( 0 ) = x20 。 对 于 数 据 方程（1） （2）成为以下形式：r1 = 1, λ1 = 0.1, r2 = 0.5, λ2 = 0.02, x10 = 25, x20 = 2, t 的终值试验后确定为 15，以便于观察。? x1 ' = x1 (1 ? 0.1x2 ) ? ? ? ? x2 ' = x2 ( ?0.5 + 0.02 x1 ) ? ? ? ? x1 ( 0 ) = 25, x2 ( 0 ) = 2 ? 首先，建立 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) ) ;[t , x] = ode45 ( ' shier ', [0 15] , [ 25 2]) ; plot ( t , x (:,1) , '? ', t , x (:, 2 ) , '? ') plot ( x (:,1) , x (:, 2 ) )得到数值解 x1 ( t ) 、x2 ( t ) 的图形（图 3-25） ，其中 x1 ( t ) 为实线； x2 ( t ) 图如图 3-26 所示。 为“*”线。相图 图 3-263-25最小值为 2.0； x1 ( t ) 的最大值为 2804，最小值为 2.0.求 x1 ( t ) 、x1 ( t ) 在一个周期内的平均值得解近似定出周期为 10.7. x1 ( t ) 的最大值为 99.3，相图 x1 ( t ) 、x2 ( t ) 是封闭曲线， 且从数值x1 = 25, x2 = 10 与理论值相符。考虑人工捕获时，设表示捕获能力的系数为 e ，相当于食饵自然生长率由 r1 降为 r1 ? e ，捕食 者的死亡率由 r2 增为 r2 + e 方程变为：? dx1 ? ? ? ? ? dt = x1 ?( r1 ? e ) ? λ1 x2 ? ? ? ? ? ? dx2 = x ? ? ( r + e ) + λ x ? ? 2 ? 2 2 1? ? dt ? ? ?仍 取r1 = 1, λ1 = 0.1, r2 = 0.5, λ2 = 0.02, x1 ( 0 ) = 25, x2 ( 0 ) = 2.设战前捕获能力系数 e = 0.3 ， 战争中降为 e = 0.1 ，则战前与战争中的模型分别为：? dx1 ? dt = x1 ( 0.7 ? 0.1x2 ) ? ? dx2 = x2 ( ?0.8 + 0.02 x1 ) ? dt ? ? x1 ( 0 ) = 25, x2 ( 0 ) = 2 ? ? ? dx1 ? dt = x1 ( 0.9 ? 0.1x2 ) ? ? dx2 = x2 ( ?0.6 + 0.02 x1 ) ? ? dt ? x1 ( 0 ) = 25, x2 ( 0 ) = 2 ? ? 分别用 m ? 文件 shier1.m 和 shier 2.m 定义上述两个方程，与 shier.m 类似。建立主程序 shark1.m , 求 解 两 个 方 程 ， 并 画 出 两 种 情 况 下 鲨 鱼 在 鱼 类 总 数 中 所 占 比 例 x2 ( t ) / ? x1 ( t ) + x2 ( t ) ? 。 shark1.m 如下： ? ?[t1 , x ] = ode45 ( ' shier1', [ 0 15] , ?[ 25 2]? ) ; ? ? [t2 , y ] = ode45 ( ' shier 2 ', [0 15] , ?[ 25 2]? ) ; ? ? x1 = x (:,1) ; x2 = x (:, 2 ) ; x3 = x2. / ( x1 + x2 ) ; y1 = y (:,1) ; y2 = y (:, 2 ) ; y3 = y2. / ( y1 + y2 ) ; plot ( t1 , x3 , '? ', t2 , y3 , '? ' )运行 shark1.m ，得图 3-27，其中，实现为战前的鲨鱼比例： “*”线为战争中的鲨鱼比例。 可以看出：战争中鲨鱼的比例比战前高，图 3-273.8 状态空间分析与泛函分析简介3.8.1 状态空间分析1、线性系统的状态空间描述 （1）状态空间概念 状态 反映系统运动状况，并可用以确定系统未来行为的信息集合。 确定系统状态的一组独立（数目最少）变量，它对于确定系统的运动状态是 以状态变量为元素构成的向量。 以状态变量为坐标所张成的空间。系统某时刻的状态可用状态空间上的点来 状态变量的一阶导数与状态变量、输入变量之间的数学关系，一般是关于系 状态变量 状态向量 状态空间 表示。 状态方程 统的一阶微分（或差分）方程组。 输出方程 输出变量与状态变量、输入变量之间的数学关系。 状态方程与输出方程合称为状态空间描述或状态空间表达式。线性定常系统状态空间表达 式一般用矩阵形式表示：必需的，也是充分的。& ? x = Ax + Bu ? ? y = Cx + Du（2）状态空间表达式的建立。系统状态空间表达式可以由系统微分方程、结构图、传递 函数等其他形式的数学模型导出。 （3）状态空间表达式的线性变换及规范化。描述某一系统的状态变量个数（维数）是确 定的，但状态变量的选择并不唯一。某一状态向量经任意满秩线性变换后，仍可作为状态向量 来描述系统。状态变量选择不同，状态空间表达式形式也不一样。利用线性变换的目的在于使 系统矩阵 A 规范化，以便于揭示系统特性，利于分析计算。满秩线性变换不改变系统的固有特 性。 根据矩阵 A 的特征根及相应的独立特征向量情况， 可将矩阵 A 化为三种规范形式： 对角形、 约当形和模式矩阵。 2. 线性系统的可控性与可观测性& （1） 系统的 （状态） 可控性。 设系统状态方程为 x = Ax + Bu ， 若在有限时间间隔 t ∈ [t 0 , t f ]内存在无约束的分段连续控制函数 u (t ) ，能使系统从任意初始状态 x (t 0 ) 转移到任意的终止状 态 x (t f ) ，则称系统是状态完全可控的，简称可控。 （2）系统输出可控性。设系统状态空间表达式为式（9.1）,若在有限时间间隔 t ∈ [t 0 , t f ] 内，存在无约束的分段连续控制函数 u (t ) ，能使系统从任意初始输出 y (t 0 ) 转移到最终内测量 到的输出 y (t f ) ，则称系统是输出完全可控的，简称输出可控。 状态可控性与输出可控性是两个不同的概念，其间没有必然联系。单输入单输出系统，若输出不可控，则系统或不可控或不可观测。 （3）系统状态可观测性。已知输出 u (t ) 及有限时间间隔 t ∈ [t 0 , t f ] 内测量到的输出 y (t ) ， 若能唯一确定初始状态 x (t 0 ) ，则称系统是完全可观测的，简称可观测。3.8.2 泛函分析泛函分析（Functional Analysis）是现代数学的一个分支，隶属于分析学，其研究 的主要对象是函数构成的空间。泛函分析是由对变换（如傅立叶变换等）的性质的研究 和对微分方程以及积分方程的研究发展而来的。使用泛函作为表述源自变分法，代表作 用于函数的函数。巴拿赫（Stefan Banach）是泛函分析理论的主要奠基人之一，而数 学家兼物理学家伏尔泰拉（Vito Volterra）对泛函分析的广泛应用有重要贡献。 泛函分析是 20 世纪 30 年代形成的数学分科。是从变分问题，积分方程和理论物理 的研究中发展起来的。它综合运用函数论，几何学，代数学的观点来研究无限维向量空 间上的函数，算子和极限理论。它可以看作无限维向量空间的解析几何及数学分析。主 要内容有拓扑线性空间等。泛函分析在数学物理方程，概率论，计算数学等分科中都有 应用，也是研究具有无限个自由度的物理系统的数学工具。泛函分析是研究拓扑线性空 间到拓扑线性空间之间满足各种拓扑和代数条件的映射的分支学科。 泛函分析的主要定理包括： 1. 一致有界定理（亦称共鸣定理），该定理描述一族有界算子的性质。 2. 谱定理包括一系列结果，其中最常用的结果给出了希尔伯特空间上正规算子的一 个积分表达，该结果在量子力学的数学描述中起到了核心作用。 3. 罕-巴拿赫定理（Hahn-Banach Theorem）研究了如何将一个算子保范数地从一 个子空间延拓到整个空间。另一个相关结果是对偶空间的非平凡性。 4. 开映射定理和闭图像定理。欢迎您转载分享： 更多精彩：