资源描述
单击此处编辑母版标题样式,单击此处编辑母版文本样式,第二级,第三级,第四级,第五级,*,本资料仅供参考，不能作为科学依据。谢谢。本资料仅供参考，不能作为科学依据。谢谢您,常微分方程数值解,杨瑞琰,第1页,0、导言,在许多实际问题中，比如物理中速率问题，人口增加问题，放射性衰变问题，经济学中边际问题等，经常包括到两个变量之间改变规律。微分方程是研究上述问题一个机理分析方法，它在科技、工程、生态、环境、人口以及经济管理等领域中有着十分广泛应用。,在应用微分方程处理实际问题时，必须经过两个阶段。一是微分方程建立，建立一个微分方程实质就是构建函数、自变量以及函数对自变量导数之间一个平衡关系。而正确地构建这种平衡关系，需要对实际问题深入浅出刻画，依据物理和非物理原理、定律或定理，作出合理假设和简化并将它升华成数学问题。,第2页,另一个是方程求解和结果分析。对一些常系数或特殊函数形式微分方程，往往能得到解析解，这对实际问题分析和应用都是有利,不过大多数变系数、非线性函数形式微分方程都是求不出解析解,此时就需要应用求解微分方程另一个主要方法数值解法。,本章简明介绍相关微分方程模型概念，微分方程数值解法和图解法，主要介绍若干建模实例，经过它们展示微分方程模型建模步骤及处理实际问题全过程。,第3页,1、实例及其数学模型,例1 海上缉私,问题,海防某部缉私艇上雷达发觉正东方向,c,海里处有一艘走私船正以速度,a,向正北方向行驶，缉私艇马上以最大速度,b,前往拦截。用雷达进行跟踪时，可保持缉私艇速度方向一直指向走私船。建立任意时刻缉私艇位置和缉私艇航线数学模型，讨论缉私艇能够追上走私船条件，求出追上时间。,第4页,建立直角坐标系如图，设在,t,=0时刻缉私艇发觉走私船，此时缉私艇位置在(0,0)，走私船位置在(,c,0)。走私船以速度,a,平行于,y,轴正向行驶，缉私艇以速度,b,按指向走私船方向行驶。在任意时刻,t,缉私艇位于,P,(,x,y,)点，而走私船抵达,Q,(,c,at,)点，直线,PQ,与缉私艇航线（图中曲线）相切，切线与,x,轴正向夹角为,。,Q,(,c,at,),P,(,x,y,),R,(,c,y,),0,y,x,c,第5页,缉私艇在,x,y,方向速度分别为,，由直角三角形,PQR,写出sin,和cos,表示式，得到微分方程：,（1）,初始条件为,（2）,这就是缉私艇位置(,x,(,t,),y,(,t,)数学模型。不过由方程（1）无法得到,x,(,t,),y,(,t,)解析解，需要用数值算法求解。我们将在后面继续讨论这个问题。,第6页,例2 弱肉强食,问题,自然界中在同一环境下两个种群之间存在着几个不一样生存方式，比如相互竞争，即争夺一样食物资源，造成一个种群趋于灭绝，而另一个趋向环境资源允许最大容量；或者相互依存，即彼此提供部分食物资源，二者和平共处，趋于一个平衡状态；再有一个关系可称之为弱肉强食，即某个种群甲靠丰富自然资源生存，而另一个群乙靠捕食种群甲为生，种群甲称为食饵（Prey），种群乙为捕食者（Predator），二者组成食饵-捕食者系统。海洋中食用鱼和软骨鱼（鲨鱼等）、美洲兔和山猫、落叶松和蚜虫等都是这种生存方式经典。这么两个种群数量是怎样演变呢？近百年来许多数学家和生态学家对这一系统进行了深入研究，建立了一系列数学模型，本节介绍是最初、最简单一个模型，它是意大利数学家Volterra在上个世纪20年代建立。,第7页,模型,用,x,(,t,)表示时刻,t,食饵（如食用鱼）密度，即一定区域内数量，,y,(,t,)表示捕食者（如鲨鱼）密度。假设食饵独立生存时（相对）增加率为常数,r,0，即 ，而捕食者存在使食饵增加率减小，设减小量与捕食者密度成正比，百分比系数为,a,0，则 。,捕食者离开食饵无法生存，设它独自存在时死亡率为常数,d,0，即 ，而食饵存在为捕食者提供了食物，使捕食者死亡率减小，设减小量与食饵密度成正比，百分比系数为,b,0，则 ，实际上，当,bx,d,时捕食者密度将增加。,给定食饵和捕食者密度初始值,x,0,y,0，由上可知,x,(,t,),y,(,t,)满足以下方程：,（3）,（3）解,x,(,t,),y,(,t,)描述了食饵和捕食者密度随时间演变过程。不过我们一样得不到,x,(,t,),y,(,t,)解析解，需要用数值算法求解。我们将在3继续讨论这个问题,第8页,2 欧拉方法和龙格库塔方法,一阶常微分方程初值问题普通形式为,y,=(x,y),axb,(4),y,(a)=,其中(x,y)是已知函数,为给定初值.,假如函数(x,y)在区域axb,-y0为,Lipschitz常数,则初值问题(4)有唯一解.,第9页,所谓数值解法,就是设法将常微分方程离散化,建立差分方程,给出解在一些离散点上近似值.,a=x,0,x,1,x,2,x,n,x,N,=b,其中剖分节点x,n,=a+nh,n=0,1,N,h称为剖分步长.数值解法就是求准确解y(x)在剖分节点x,n,上近似值y,n,y(x,n,),n=1,2,n.,假设初值问题(4)解y=y(x)唯一存在且足够光滑.对求解区域a,b做剖分,我们采取数值积分方法来建立差分公式.,2.1 结构数值解法基本思想,在区间x,n,x,n+1,上对方程(4)做积分,则有,第10页,对右边积分应用左矩形公式，则有,梯形公式,o,x,y,a,b,左矩形公式,y=,(x),右矩形公式,中矩形公式,第11页,对右边积分应用左矩形公式，则有,所以，建立节点处近似值y,n,满足差分公式,称之为,Euler公式,.,称为,梯形公式,.,若对(6)式右边积分应用梯形求积公式，则可导出差分公式,第12页,利用Euler方法求初值问题,解,此时Euler公式为,称为,Euler中点公式,或称,双步Euler公式.,若在区间x,n-1,x,n+1,上对方程(4)做积分,则有,对右边积分应用中矩形求积公式，则得差分公式,例3,数值解.此问题准确解是y(x)=x/(1+x,2,).,第13页,分别取步长h=0.2,0.1,0.05,计算结果以下,第14页,h,x,n,y,n,y(x,n,),y(x,n,)-y,n,h=0.2,0.00,0.40,0.80,1.20,1.60,2.00,0.00000,0.37631,0.54228,0.52709,0.46632,0.40682,0.00000,0.34483,0.48780,0.49180,0.44944,0.40000,0.00000,-0.03148,-0.05448,-0.03529,-0.01689,-0.00682,h=0.1,0.00,0.40,0.80,1.20,1.60,2.00,0.00000,0.36085,0.51371,0.50961,0.45872,0.40419,0.00000,0.34483,0.48780,0.49180,0.44944,0.40000,0.00000,-0.01603,-0.02590,-0.01781,-0.00928,-0.00419,h=0.05,0.00,0.40,0.80,1.20,1.60,2.00,0.00000,0.35287,0.50049,0.50073,0.45425,0.40227,0.00000,0.34483,0.48780,0.49180,0.44944,0.40000,0.00000,-0.00804,-0.01268,-0.00892,-0.00481,-0.00227,第15页,Euler中点公式则不然,计算y,n+1,时需用到前两步值y,n,y,n-1,称其为,两步方法,两步以上方法统称为,多步法,.,在Euler公式和梯形公式中,为求得y,n+1,只需用到前一步值y,n,这种差分方法称为,单步法,这是一个自开始方法.,隐式公式中,每次计算y,n+1,都需解方程,要比显式公式需要更多计算量,但其计算稳定性很好.,在Euler公式和Euler中点公式中,需要计算y,n+1,已被显式表示出来,称这类差分公式为,显式公式,而梯形公式中,需要计算y,n+1,隐含在等式两侧,称其为,隐式公式,.,第16页,从数值积分角度来看,梯形公式,计算数值解精度要比Euler公式好,但它属于隐式公式,不便于计算.,实际上,常将Euler公式与梯形公式结合使用:,2.2,改进Euler方法,第17页,由迭代法收敛角度看，当,(,是给定精度要求)时,取,就能够确保迭代公式收敛,而当h很小时,收敛是很快.,而且,只要,通常采取只迭代一次算法:,称之为,改进Euler方法.,这是一个单步显式方法.,改进Euler方法也能够写成,第18页,y,=y-2x/y ,0 x1,数值解,取步长h=0.1.准确解为y(x)=(1+2x),1/2,.,例4,求初值问题,y,(0)=1,解,(1)利用Euler方法,第19页,计算结果以下:,(2)利用改进Euler方法,第20页,n,x,n,Euler方法y,n,改进Euler法y,n,准确解y(xn),0,1,2,3,4,5,6,7,8,9,10,0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1,1,1.1,1.191818,1.277438,1.358213,1.435133,1.508966,1.580338,1.649783,1.717779,1.784770,1,1.095909,1.184096,1.266201,1.343360,1.416402,1.485956,1.552515,1.616476,1.678168,1.737869,1,1.095445,1.183216,1.264991,1.341641,1.414214,1.483240,1.549193,1.612452,1.673320,1.732051,第21页,在节点x,n+1,误差y(x,n+1,)-y,n+1,不但与y,n+1,这一步计算相关,而且与前n步计算值y,n,y,n-1,y,1,都相关.,为了简化误差分析,着重研究进行一步计算时产生误差.即假设y,n,=y(x,n,),求误差y(x,n+1,)-y,n+1,这时误差称为,局部截断误差,它能够反应出差分公式精度.,2.3,差分公式误差分析,假如单步差分公式局部截断误差为O(h,p+1,),则称该公式为,p阶方法,.这里p为非负整数.显然,阶数越高,方法精度越高.,研究差分公式阶主要伎俩是Taylor展开式,一元函数和二元函数Taylor展开式为:,第22页,另外,在y,n,=y(x,n,)条件下,考虑到y(x)=(x,y(x),则有,y(x,n,)=(x,n,y(x,n,)=(x,n,y,n,)=,n,y(x,n,)=(x,n,y(x,n,)=,x,(x,n,y,n,)+,y,(x,n,y,n,)(x,n,y,n,),第23页,y,n+1,=y,n,+h(x,n,y,n,),对Euler方法，有,=y,n,+(x,n,y,n,)h+O(h,2,),从而有:y(x,n+1,)-y,n+1,=O(h,2,),所以Euler方法是一阶方法.,再看改进Euler方法,因为,可得,第24页,所以,改进Euler方法是二阶方法.,而,从而有:y(x,n+1,)-y,n+1,=O(h,3,),2.4 Taylor,展开方法,设y(x)是初值问题（4)准确解,利用Taylor展开式可得,第25页,称之为p阶Taylor展开方法.,所以，可建立节点处近似值y,n,满足差分公式,其中,第26页,所以,此差分公式是p阶方法.,因为Taylor展开方法包括很多复合函数(x,y(x)导数计算,比较繁琐,因而极少直接使用,经惯用它为多步方法提供初始值.然而,Taylor展开方法给出了一个结构单步显式高阶方法路径.,Euler方法可写为,可见,公式局部截断误差为:y(x,n+1,)-y,n+1,=O(h,p+1,).,2.5 Runge-Kutta,方法,第27页,结构差分公式,改进Euler方法可写为,其中,i,i,ij,为待定参数.,若此公式局部截断误差为,第28页,因为,y,n+1,=y,n,+h,1,n,+h,2,(,n,+h,xn,+h,n,yn,)+O(h,3,),O(,h,3,),称此公式为,p阶Runge-kutta方法,简称,p阶R-K方法,.,对于p=2情形,应有,=y,n,+h(,1,+,2,),n,+h,2,2,(,xn,+,n,yn,)+O(h,3,),所以,只要令,1,+,2,=1,2,=1/2,2,=1/2 (8),第29页,普通地,参数由(8)确定一族差分公式(7)统称为,二阶R-K方法.,称之为,中点公式,或可写为,若取,=1,则得,1,=,2,=1/2,=1,此时公式(7)就是改进Euler公式;,若取,1,=0,则得,2,=1,=1/2,公式(7)为,高阶R-K公式可类似推导.,下面列出惯用三阶、四阶R-K公式.,第30页,四阶标准R-K公式,三阶R-K公式,第31页,解,四阶标准R-K公式为,例3,用四阶标准R-K方法求初值问题,y,=y-2x/y ,0 x1,y,(0)=1,数值解,取步长h=0.2.,计算结果以下:,第32页,n,x,n,y,n,y(x,n,),n,x,n,y,n,y(x,n,),0,1,2,0.0,0.2,0.4,1.00,1.1832,1.3417,1.00,1.1832,1.3416,3,4,5,0.6,0.8,1.0,1.4833,1.6125,1.7321,1.4832,1.6125,1.7321,也能够结构隐式R-K方法,其普通形式为,称之为,p级隐式R-K方法,同显式R-K方法一样确定参数.如,第33页,是二级二阶隐式R-K方法,也就是梯形公式.不过p级隐式R-K方法阶能够大于p,比如,一级隐式中点公式为,或写为,它是二阶方法.,2.6,变步长,Runge-Kutta,方法,以p阶R-K方法为例讨论.设从x,n,以步长h计算y(x,n+1,)近似值为,局部截断误差为,其中,C是与h无关常数.,第34页,假如将步长减半,取h/2为步长,从x,n,经两步计算得到y(x,n+1,)近似值记为,其局部截断误差为,于是有,从而,得到事后误差预计,可见,当,成立时,可取,.不然,应将步长再次减半进行计算.,第35页,求解初值问题,单步显式方法可一统一写为以下形式,y,n+1,=y,n,+h(x,n,y,n,h)（9),对于Euler方法,有,2.7,单步方法收敛性,y,=(x,y),axb,y,(a)=,其中(x,y,h)称为,增量函数,.,(x,y,h)=(x,y),对于改进Euler方法,有,(x,y,h)=1/2(x,y)+(x+h,y+h(x,y),第36页,设y(x)是初值问题（4)解,y,n,是单步法（9)产生近似解.假如对任意固定点x,n,都有,y(x,n,),则称单步法（9)是收敛.,可见,若方法（9)是收敛,则当h0时,整体截断误差e,n,=y(x,n,)-y,n,将趋于零.,定理,设单步方法（9)是p1阶方法,增量函数(x,y,h)在区域axb,-yn)改变均不超出,则称此差分方法是,绝对稳定,.,讨论数值方法稳定性,通常仅限于经典试验方程,y=y,其中是复数且Re()0.,在复平面上,当方法稳定时要求变量h取值范围称为方法,绝对稳定域,它与实轴交集称为,绝对稳定区间.,第40页,将Euler方法应用于方程y=y,得到,设在计算y,n,时产生误差,n,计算值y,n,=y,n,+,n,则,n,将对以后各节点值计算产生影响.记y,m,=y,m,+,m,mn,由上式可知误差,m,满足方程,m,=(1+h),m-1,=(1+h),m-n,n,mn,对隐式单步方法也可类似讨论.如将梯形公式用于方程y=y,则有,y,n+1,=y,n,+h/2(y,n,+y,n+1,),y,n+1,=(1+h)y,n,可见,若要|,m,|,n,|,必须且只须|1+h|1,所以Euler法绝对稳定域为|1+h|1,绝对稳定区间是-2Re()h0.,解出y,n+1,得,第41页,类似前面分析,可知绝对稳定区域为,因为Re()0,所以此不等式对任意步长h恒成立,这是隐式公式优点.,一些惯用方法绝对稳定区间为,方 法,方法阶数,稳 定 区 间,Euler方法,梯形方法,改进Euler方法,二阶R-K方法,三阶R-K方法,四阶R-K方法,1,2,2,2,3,4,(-2 ,0),(-,0),(-2 ,0),(-2 ,0),(-2.51 ,0),(-2.78 ,0),第42页,解,因y,0,=1,计算得y,10,=1024,而y(1)=9.35762310,-14,.,例4,考虑初值问题,y,=-30y ,0 x1,y,(0)=1,取步长h=0.1,利用Euler方法计算y,10,y(1).y(x)=e,-30 x,这是因为,h=-3不属于Euler方法绝对稳定区间.,若取h=0.01,计算得y,100,=3.23447710,-16,.,若取h=0.001,计算得y,1000,=5.91199810,-14,.,若取h=0.0001,计算得y,10000,=8.94505710,-14,.,若取h=0.00001,计算得y,100000,=9.315610,-14,.,第43页,单步显式方法稳定性与步长亲密相关,在一个步长下是稳定差分公式,取大一点步长就可能是不稳定.,收敛性是反应差分公式本身截断误差对数值解影响;稳定性是反应计算过程中舍入误差对数值解影响.只有即收敛又稳定差分公式才有实用价值.,2.9,线性多步方法,因为在计算y,n+1,时,已经知道y,n,y,n-1,及(x,n,y,n,),(x,n-1,y,n-1,),利用这些值结构出精度高、计算量小差分公式就是线性多步法.,2.9.1,利用待定参数法结构线性多步方法,r+1步线性多步方法普通形式为,第44页,当,-1,0时,公式为隐式公式,反之为显式公式.参数,i,i,选择标准是使方法局部截断误差为,y(x,n+1,)-y,n+1,=O(h),r+2,选取参数,0,1,2,使三步方法,y,n+1,=y,n,+h(,0,n,+,1,n-1,+,2,n-2,),这里,局部截断误差是指,在y,n-i,=y(x,n-i,),i=0,1,r前提下,误差y(x,n+1,)-y,n+1,.,为三阶方法.,例5,解,设y,n,=y(x,n,),y,n-1,=y(x,n-1,),y,n-2,=y(x,n-2,),则有,第45页,n,=(x,n,y(x,n,)=y(x,n,),y(x,n+1,)=y(x,n,)+hy(x,n,)+1/2h,2,y(x,n,)+1/6h,3,y(x,n,),于是有,若使:y(x,n+1,)-y,n+1,=O(h,4,),只要,0,1,2,满足:,n-1,=(x,n-1,y(x,n-1,)=y(x,n-1,)=y(x,n,-h),=y(x,n,)-hy(x,n,)+1/2h,2,y(x,n,)-1/6h,3,y,(4),(x,n,)+O(h,4,),n-2,=y(x,n,)-2hy(x,n,)+2h,2,y(x,n,)-4/3h,3,y,(4),(x,n,)+O(h,4,),y,n+1,=y(x,n,)+h(,0,+,1,+,2,)y(x,n,)-h,2,(,1,+2,2,)y(x,n,),+h,3,(1/2,1,+2,2,)y(x,n,)-h,4,/6(,1,+8,2,)y,(4),(x,n,)+O(h,5,),+1/24h,4,y,(4),(x,n,)+O(h,5,),第46页,=1,0,+,1,+,2,=1,1,+2,2,=-1/2,1,+4,2,=1/3,于是有三步三阶显式差分公式,设p,r,(x)是函数(x,y(x)某个r次插值多项式,则有,解之得:,y,n+1,=y,n,+h/12(23,n,-16,n-1,+5,n-2,),因为,2.9.2,利用数值积分结构线性多步方法,其中,第47页,选取不一样插值多项式p,r,(x),就可导出不一样差分公式.下面介绍惯用,Adams公式.,设已求得准确解y(x)在步长为h等距节点x,n-r,x,n,上近似值y,n-r,y,n,记,k,=(x,k,y,k,),利用r+1个数据(x,n-r,n-r,),(x,n,n,)结构r次Lagrange插值多项式,由此,可建立差分公式,1.Adams显式公式,其中,第48页,由此,可建立差分公式,因为,h,rj,则有,称之为,r+1步Adams显式公式.,第49页,下面列出几个带有局部截断误差主项Adams显式公式,r=0 y,n+1,=,y,n,+h,n,+(1/2)h,2,y(x,n,),2.Adams隐式公式,r=1 y,n+1,=,y,n,+(h/2)(3,n,-,n-1,)+(5/12)h,3,y(x,n,),r=2 y,n+1,=,y,n,+(h/12)(23,n,-16,n-1,+5,n-2,),+(3/8)h,4,y,(4),(x,n,),r=3 y,n+1,=,y,n,+(h/24)(55,n,-59,n-1,+37,n-2,-9,n-3,),+(251/720)h,5,y,(5),(x,n,),假如利用r+1个数据(x,n-r+1,n-r+1,),(x,n+1,n+1,)结构r次Lagrange插值多项式p,r,(x),则可导出数值稳定性好隐式公式,称为,Adams隐式公式,其普通形式为,第50页,其中系数为,下面列出几个带有局部截断误差主项Adams隐式公式,r=0 y,n+1,=,y,n,+h,n+1,-(1/2)h,2,y,(x,n,),r=1 y,n+1,=,y,n,+(h/2)(,n,+,n+1,)-(1/12)h,3,y,(x,n,),r=2 y,n+1,=,y,n,+(h/12)(5,n+1,+8,n,-,n-1,),-(1/24)h,4,y,(4),(x,n,),r=3 y,n+1,=,y,n,+(h/24)(9,n+1,+19,n,-5,n-1,+,n-2,),-(19/720)h,5,y,(5),(x,n,),第51页,3.Adams预估-校正公式,由显式公式提供一个预估值，再用隐式公式校正一次，求得数值解，称为,预估校正方法,。,校正 y,n+1,=,y,n,+(h/24)(9,n+1,+19,n,-5,n-1,+,n-2,),普通预估公式和校正公式都采取同阶公式。比如：,预估 y,n+1,=,y,n,+(h/24)(55,n,-59,n-1,+37,n-2,-9,n-3,),n+1,=(x,n+1,y,n+1,),n=3,4,称为四阶Adams预估校正公式.实际计算时通惯用四阶单步方法(如四阶R-K公式)为它提供起始值y,1,y,2,y,3,.,例6,用四阶Adams预估校正公式求解初值问题,第52页,y,=y-2x/y,0 x1,y(0)=1,取步长h=0.1.,解 用四阶R-K公式提供起始值,计算结果以下,x,n,R-k法y,n,预估值y,n,校正值y,n,准确值y(xn),0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1,1,1.095446,1.183217,1.264912,1.341551,1.414045,1.483017,1.548917,1.612114,1.672914,1.731566,1.341641,1.414213,1.483239,1.549192,1.612450,1.673318,1.732048,1,1.095445,1.183216,1.264991,1.341641,1.414214,1.483240,1.549193,1.612452,1.673320,1.732051,第53页,3 RK方法MATLAB实现,对于微分方程（组）初值问题,龙格库塔方法可用以下MATLAB命令实现其计算：,t，x=ode23(f，ts，x0，options),t，x=ode45(f，ts，x0，options),其中ode23用是3级2阶龙格库塔公式，ode45用是以Runge-Kutta-Fehberg命名 5级4阶公式。,第54页,命令输入f是待解方程写成函数m文件：,function dx=f(t,x),dx=f1;f2;,;fn;,若ts=t0,t1,t2,tf，则输出在指定时刻t0,t1,t2,tf函数值；等分点时用ts=t0:k:tf，输出在t0,tf内等分点处函数值。x0为函数初值（n维向量）。options可用于设定误差限（options缺省时设定相对误差10-3，绝对误差10-6），命令为：,options=odeset(reltol,rt,abstol,at),其中rt，at分别为设定相对误差和绝对误差。,命令输出t为指定ts，x为对应函数值（n维向量）。,注意，计算步长是依据误差限自动调整，并不是输入中指定ts分点。,第55页,下面用MATLAB软件处理,1提出两个问题,例1 海上缉私（续）,模型数值解,1.,设,a,=20(海里/小时)，,b,=40(海里/小时)，,c,=15(海里)，由模型（1），（2）求任意时刻缉私艇位置及缉私艇航线。,对于给出,a,b,c,用MATLAB求数值解时，记,x,(1)=,x,x,(2)=,y,x,=(,x,(1),x,(2)T。编写以下m文件：,function dx=jisi(t,x)%建立名为jisi函数m文件,a=20;b=40;c=15;,s=sqrt(c-x(1)2+(a*t-x(2)2);,dx=b*(c-x(1)/s;b*(a*t-x(2)/s;%以向量形式表示方程（1）,第56页,然后运行以下程序：,ts=0:0.05:0.5;%设定t起终点及中间等分点，终点可先作试探，再按照x(t),c=15调整到0.5,x0=0,0;%输入x，y初始值（2）,t,x=ode45(jisi,ts,x0);%调用ode45计算,t,x%输出t,x(t),y(t),plot(t,x),grid,%按照数值输出作x(t),y(t)图形,gtext(x(t),gtext(y(t),pause,plot(x(:,1),x(:,2),grid,%作y(x)图形,gtext(x),gtext(y),第57页,得到数值结果,x,(,t,),y,(,t,)为缉私艇位置，列入表1。走私船位置记作,x,1(,t,),y,1(,t,)，显然,x,1(,t,)=,c,=15，,y,1(,t,)=,at,=20,t,，将,y,1(,t,)列入表1最终一列。可知当,t,=0.5（小时），,x,y,与,x,1,y,1几乎一致，认为缉私艇追上走私船。,x,(,t,),y,(,t,)及,y,(,x,)图形见图2，,y,(,x,)为缉私艇航线。,t,x,(,t,),y,(,t,),y,1,(,t,),0,0,0,0,0.05,1.9984,0.0698,1.0,0.10,3.9854,0.2924,2.0,0.15,5.9445,0.6906,3.0,0.20,7.8515,1.2899,4.0,0.25,9.6705,2.1178,5.0,0.30,11.3496,3.,6.0,0.35,12.8170,4.5552,7.0,0.40,13.9806,6.1773,8.0,0.45,14.7451,8.0273,9.0,0.50,15.0046,9.9979,10.0,a=20，b=40，c=15数值解,x(t),y(t)和y1(t),第58页,a,=20，,b,=40，,c,=15 时,x,(,t,),y,(,t,)和,y,(,x,),图形,第59页,2.设,b,，,c,不变，而,a,变大为30，35，,靠近40(海里/小时)，观察解改变,修改,a,输入，并对应地延长,t,终点。设,a,=35，,t,终点经试探，调整为1.6适当。表2是计算结果，其中,x,(,t,),y,(,t,)有两列数字，左边是用“缺省”精度（即相对误差10-3，绝对误差10-6）计算，中间,y,1(,t,)=,at,=35,t,是走私船抵达位置。可知,t,=1.3,1.4,1.5时缉私艇位置,x,15,但,y,与,y,1(,t,)相差甚远，,t,=1.6时,x,y,与,x,1,y,1也有差距，这是累积误差造成。可利用ode45控制参数options提升精度（上面“调用ode45计算”用以下程序代替），如设,第60页,opt=odeset(reltol,1e-6,abstol,1e-9);,t,x=ode45(jisi,ts,x0,opt);,得到,x,(,t,),y,(,t,)，与走私船抵达位置,x,1(,t,)，,y,1(,t,)相对照，可知,t,=1.6时,x,y,与,x,1,y,1几乎一致，认为缉私艇追上走私船。,x,(,t,),y,(,t,)及,y,(,x,)图形见图，,y,(,x,)为缉私艇航线，当,x,靠近15时航线几乎是正北方向，形成沿走私船逃向追赶态势。,读者可尝试寻求缉私艇追上走私船判断标准，设计程序并计算。,第61页,a,=35，,b,=40，,c,=15 时,x,(,t,),y,(,t,)和,y,(,x,),图形,第62页,模型解析解,要想得到缉私艇拦截到走私船准确时间和位置，必须对模型（1）做深入分析，设法得到某种形式解析解。,将（1）两式相除，消去,dt,得到,（27）,（27）式对,x,求导得,（28）,为消去式中,dt,，利用曲线弧长微分，缉私艇速度，以及微分关系，得,第63页,（29）,这是关于,y,(,x,)二阶微分方程，若令 ，可化为,p,(,x,)一阶微分方程,（30）,依题意其初始条件为,（31）,（30）为可分离变量方程，轻易求得在（31）下解为,（32）,第64页,设，即走私船速度,a,小于缉私艇速度,b,，对（32）积分并注意到得,（35）,这就是缉私艇航线解析表示式。,由（35）式知，当,x=c,时，,这也是走私船,y,坐标。因为走私船速度是,a,，所以缉私艇拦截到走私船时间为,第65页,例2 弱肉强食（续）,模型数值解,为了考查模型（3）描述食饵和捕食者密度随时间演变过程，选取模型参数,r,=1，,d,=0.5，,a,=0.1，,b,=0.02，,x,0=25，,y,0=2，用MATLAB求（3）数值解，记,x,(1)=,x,x,(2)=,y,x,=(,x,(1),x,(2),T,。编写以下m文件：,function xdot=shier(t,x),r=1;d=0.5;a=0.1;b=0.02;,xdot=diag(r-a*x(2),-d+b*x(1)*x;%以向量形式表示方程（3）,第66页,然后运行以下程序：,ts=0:0.1:15;%t终值经试验后定为15。,x0=25,2;,t,x=ode45(shier,ts,x0);,t,x,plot(t,x),grid,gtext(fontsize12x(t),gtext(fontsize12y(t),%将标识x(t),y(t)字体放大,pause,plot(x(:,1),x(:,2),grid,xlabel(x),ylabel(y),第67页,得到数值结果从略，,x,(t),y,(t)及,y,(,x,)图形见图4，曲线,y,(,x,)称,相轨线,，图形称,相图,。能够猜测，,x,(,t,),y,(,t,)是周期函数，,y,(,x,)是封闭曲线，而且从数值结果能够确定周期约为10.7，还能够用数值积分算出,x,(,t,),y,(,t,)在一个周期平均值，得到，它们代表食饵和捕食者平均密度。,第68页,
展开阅读全文
浙ICP备2021020529号-1 | 浙B2-20240490 
