二阶一阶非线性微分方程程可不可以没有一阶导数?！！！！很迷！
点击联系发帖人
时间：2018-03-15 01:26
第四章 微分方程与差分方程方法_甜梦文库
第四章 微分方程与差分方程方法
第四章 微分方程与差分方程方法第一节 微分方程模型我们在数学分析中所研究的函数,是反映客观现实世界运动过程中量与量之 间的一种关系,但我们在构造数学模型时,遇到的大量实际问题往往不能直接写 出量与量之间的关系,却能比较容易地建立这些变量和它们的导数(或微分)间的 关系式,这种联系着自变量、 未知函数及其导数(或微分)的关系式称为微分方程.§4.1.1 微分方程简介这一节,我们将介绍关于微分方程的一些基本概念. 一、微分方程的阶数 首先我们具体的来看一个微分方程的例子. 例 4-1 物体冷却过程的数学模型 将某物体放置于空气中,在时刻 t ? 0 ,测量得它的温度为 u 0 ? 150 0 C ,10 分 钟后测量得温度为 u 1 ? 100 0 C .我们要求决定此物体的温度 u 和时间 t 的关系,并 计算 20 分钟后物体的温度.这里我们假定空气的温度保持为 u ? ? 24 0 C . 解:根据物理学中的牛顿冷却定律可知,热量总是从温度高的物体向温度低的物 体传导;一个物体的温度变化速度与这一物体的温度与其所在介质温度的差值成 正比.设物体在时刻 t 的温度为 u ? u (t ) ,则温度的变化速度可以用 们得到描述物体温度变化的微分方程du dt ? ? k (u ? u ? ) du dt来表示.我(4.1.1)du dt其中 k ? 0 是比例常数. 方程(4.1.1)中含有未知函数 u 及它的一阶导数 ,这样的方程,我们称为一阶微分方程.微分方程中出现的未知函数最高阶导数的阶数称为微分方程的阶 数. 方程d y dt3 3?bdy dt? cy ? f ( t )(4.1.2)中未知函数最高阶导数的阶数是三阶,则方程(4.1.2)称为三阶微分方程. 二、常微分方程与偏微分方程 如果在微分方程中,自变量的个数只有一个,我们称这种微分方程为常微分 方程;自变量的个数为两个或两个以上的微分方程称为偏微分方程. 方程 ? T2?x2?? T2?y2?? T2?z2? 0(4.1.3)就是偏微分方程的例子,其中 T 是未知函数, x 、 y 、 z 都是自变量.而方程 (4.1.1)(4.1.2)都是常微分方程的例子. 三、线性与非线性微分方程 如果 n 阶常微分方程F ( x, y, dy dx ,? , d y dxn n) ? 0(4.1.4)的左端为关于未知函数 y 及其各阶导数的线性组合,则称该方程为线性微分方程, 否则称为非线性方程.一般的 n 阶线性微分方程具有形式d y dxn n? a1 ( x )dn ?1ydxn ?1? ? ? a n ?1 ( x )dy dx? an (x) y ? f (x)(4.1.5)其中 a i ( x ), f ( x ) ( i ? 1 ? n ) 是关于 x 的已知函数. 当 f ( x ) ? 0 时,称(4.1.5)为 n 阶齐次线性微分方程;当 f ( x ) ? 0 时,称(4.1.5)为n阶非齐次线性微分方程. 例如,方程(4.1.2)是三阶线性微分方程.而方程dy ? dy ? ? y ?0 ? ? ?t dx ? dx ?2(4.1.6)是一阶齐次非线性方程. 四、微分方程初边值问题 我们把含有 n 个独立的任意常数 c1 , c 2 , ? c n 的解y ? ? ( x , c1 , c 2 , ? c n )称为 n 阶微分方程(4.1.4)的通解.为了确定微分方程一个特定的解,我们通常给出 这个解所必须满足的条件,这就是定解条件.常见的定解条件是初始条件. n 阶微 分方程(4.1.4)的初始条件是指如下的 n 个条件:y ( x0 ) ? y0 , d y ( x0 ) dx ? y0 ,?(1)dn ?1y ( x0 )n ?1dx? y0( n ? 1)求解微分方程满足初始条件的解,称为初值问题.求解初值问题的过程,就是 通过初始条件确定通解中的常数,从而求得满足初始条件的特解的过程. 若方程所给出的定解条件,既有自变量初始时刻的值,也有自变量取终值时 的值,则称该问题为边值问题. §4.1.2 微分方程的求解及 Matlab 实现线性微分方程和低阶特殊微分方程往往可以通过解析解的方法求解,但一般 的非线性微分方程是没有解析解的,即使可以求得解析解,参数也很难确定,需要 用数值解的方式求解.具体的求解方法很多,本节我们主要介绍如何利用 matlab 来求解微分方程的解析解和数值解.§4.1.2.1 微分方程的解析解一、基本理论 有些简单的常微分方程可用一些技巧,如分离变量法、积分因子法、常数变 异法、降阶法等,化为可直接积分的方程从而求得显示解.例如,一阶常系数线性 常微分方程dy dt ? ay ? b (a ? 0)可化为dy ay ? b ? dt两边积分可得 y ( t ) 的通解为 y ? C e a t ? 一般的常系数线性微分方程d y dxn nb a.? a1dn ?1ydxn ?1? ? ? a n ?1dy dx? an y ? f ( x)(4.1.7)的解满足叠加原理,即方程(4.1.7)的通解是对应齐次方程d y dxn n? a1dn ?1ydxn ?1? ? ? a n ?1dy dx? an y ? 0(4.1.8)的通解与非齐次方程(4.1.7)的一个特解的和. 一阶常系数线性常微分方程总可用这一思路求得显示解. 高阶线性常系数微分方程可用特征根法求对应齐次微分方程(4.1.8)的通 解.(4.1.8)对应的代数特征方程s ? a1 sn n ?1? a2 sn?2? ? ? a n ?1 s ? a n ? 0若方程的特征根 s i 均可以求出,且两两相异,则方程(4.1.7)的解可以表示为y ( x ) ? C 1es1 x? C 2es2 x? ? ? C nesn x? ? (x)其中, C i 为待定系数, ? ( x ) 是满足方程(4.1.8)的一个特解,这个特解可以通过常数 变异法,比较系数法得到. s i 有重根的情况也可以写出相应的解析解形式,这里我 们不做过多介绍,因为我们下面要介绍一种更简单的计算机求解法. 二、用 matlab 求解线性常系数微分方程解析解 dsolve 函数调用格式y= d so lve('f 1 ','f 2 ', ? , ' f m ') y= d so lve('f 1 ','f 2 ', ? , ' f m ','x ') 指明自变量其中, f i 既可以描述微分方程,也可以描述初边值条件.在描述微分方程时,用字 母 D 代表微分,后面加数字表示微分的阶数,例如 D4y 表示 y ( 4 ) ( t ) ,也可以用 D2y(0)=3 这类记号来表示初值条件 y ''(0) ? 3 .一般而言任何 D 后面所跟的字母 为因变量,自变量默认为 t ,也可以在函数调用时指明. 例 4-2 求方程 解:输入命令 &&syms t,x; x=dsolve(‘D2x-2*Dx-3*x=exp(-t)’) latex(x) 运行结果 x =exp(-t)*C2+exp(3*t)*C1-1/4*t*exp(-t) ans ={e^{3t}}{\it C2}+{e^{-t}}{\it C1}-1/4\,t{e^{-t}} 原方程的通解为x ? C 1e3td x dt22?2dx dt? 3x ? e?t的通解.%定义符号变量 %求解方程 %用 LATEX 语句显示结果? C 2e?t?1 4te?t其中, C i 为任意常数.若给出初始或边界条件,则可以通过这些条件建立方程,求 出 C i 的值. 仍考虑上面的微分方程,假设已知 x (0 ) ? 3, 令求满足该方程的特解. && x=dsolve(‘D2x-2*Dx-3*x=exp(-t)’,’x(0)=3’,’Dx(0)=2’) 运行结果 x =27/16*exp(-t)+21/16*exp(3*t)-1/4*t*exp(-t)dx dt (0 ) ? 2,则可以通过下面的命 所以，满足这一初值问题的特解是x ? 21 16 e3t?27 16e?t?1 4te?t对于线性方程组我们可以通过定义多维的因变量同样求解. 例 4-3 试求解下面的线性微分方程? dx ? 2x ? 3 y ? 3z ? dt ? ? dy ? 4x ? 5 y ? 3z ? dt ? ? dz ? 4x ? 4y ? 2z ? ? dt解:输入命令 &&[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')； x=collect(simple(x)) %化简 x 的表达式 y=collect(simple(y)) z=collect(simple(z)) 结果为x ? C 2e y ? C 2e z ? C 3e?t ?t? C 3e ? C 3e ? C 1e2t ?2t2t? C 1e2t?2 t部分特殊的非线性微分方程也可以用 dsolve()函数来求解析解.这样的方程 描述方式和前面介绍的线性微分方程是一致的.下面我们通过例子来演示非线性 方程的解析解求解问题,同时还将演示不能求解的例子. 例 4-4 试求出一阶非线性微分方程 解: 输入命令 &&x=dsolve('Dx=x*(1-x^2)') 运行结果 x= 1/(1+exp(-2*t)*C1)^(1/2) -1/(1+exp(-2*t)*C1)^(1/2) 即该方程的解析解为 x ( t ) ? ?1 1 ? C 1e?2 tdx dt? x ( t )(1 ? x ( t )) 的解析解.2.但是如果稍微改变原方程,例如将等号右侧加上 1,则可以用下面语句试解 方程.我们会发现该方程是没有解析解的. &&x=dsolve('Dx=x*(1-x^2)+1') 运行结果 Warning: Explicit solutio implicit solution returned. & In dsolve at 312 x= t+Int(-1/(_a-_a^3+1),_a = .. x)+C1 = 0 从上例我们可以看出,可以求得解析解的方程是很有限的,对于一般的非线 性方程只能用数值解法去求解.下面我们介绍微分方程数值解的求解方法.§4.1.2.2 微分方程的数值解一、基本理论 (一) 数值解的定义 考虑一阶常微分方程组初值问题? y ' ? f (t , y ) ? ? y (t0 ) ? y 0 t0 ? t ? t f(4.1.9)1 其中 y ? ( y 1 , y 2 , ? , y m ) T , f ? ( f 1 , f 2 , ? , f m ) T , y 0 ? ( y 0 , y 02 , ? , y 0m ) T ,这里 T 表示转置.微分方程的数值解,是指不求出解 y ( t ) 的解析表达式,而是寻求 y ( t ) 在一系 列离散结点 t 0 ? t1 ? ? ? t n ? t f 上的近似值 y k ( k ? 0,1,? , n ) .称 hk ? t k ? 1 ? t k 为步长. 通常,取为常数 h ? ( t n ? t 0 ) / n ,此时结点变为等距结点,有t k ? t k ?1 ? h ? t 0 ? kh.(二) 求解数值解的方法 1.欧拉法 基本思想:在结点处用差商近似替代导数y (t k ? h ) ? y (t k ) h ? y '( t k ) ? f ( t k , y ( t k )) ? f ( t k , y k )这样导出欧拉格式的计算公式? y k ?1 ? y k ? h f ( t k , y k ) ? ? y 0 ? y (t0 ) k ? 0,1, ? n ? 1(4.1.10)欧拉法可以求解各种形式的微分方程,但是它只有一阶精度,所以实际应用 效率比较差. 2.改进的欧拉法 基本思想: 使用数值积分离散化.对方程 y ' ? f ( t , y ) ,两边由 t k 到 t k ? 1 积分, 得到y ( t k ?1 ) ? y ( t k ) ?t k ?1 tk?f ( t , y ( t )) d t对右端的定积分用数值积分方法做离散化,可得计算公式,如用矩形公式可 得欧拉公式,若用梯形公式可得改进的欧拉公式,下面我们用梯度公式来做数值 积分y ( t k ?1 ) ? y ( t k ) ??t k ?1 tkf ( t , y ( t )) d t ?t k ?1 ? t k 2[ f ( t k , y ( t k )) ? f ( t k ? 1 , y ( t k ? 1 ))]故有如下迭代计算公式:h ? ? y k ? 1 ? y k ? [ f ( t k , y k ) ? f ( t k ? 1 , y k ? 1 )] 2 ? ? y ? y (t ) 0 ? 0(4.1.11)实际应用时,与欧拉公式结合使用:? y k( 0 )1 ? y k ? h f ( t k , y k ) ? ? ? ( l ? 1) h (l ) ? y k ? 1 ? y k ? [ f ( t k , y k ) ? f ( t k ? 1 , y k ? 1 )] ? 2l ? 0,1, 2, ?(4.1.12)对于已给的精度 ? ,当满足 y k( l??11) ? y k( l?)1 ? ? 时,取 y k ? 1 ? y k( l??11) ,然后继续下一 步 y i ? 2 的计算. 3.Runge-Kutta 方法 基本思想:利用 Taylor 展式进行离散化.假设(4.1.9)中的 f ( t , y ) 充分光滑,将y ( t k ?1 ) 在 t k点作 Taylor 展开:y ( t k ? 1 ) ? y ( t k ) ? h y '( t k ) ? h22!y ''( t k ) ? ? ?hpy( p)p!(t k ) ? ?其中y '( t ) ? f ( t , y ( t )) y ''( t ) ? [ f ( t , y ( t ))] 't ? f t ? f ? f y ? y( p)( t ) ? [ f ( t , y ( t ))] t( p)对照标准形式 y k ? 1 ? y k ? h ? ( t k , h ) .若取 ? ( t , h ) ? y '( t ) ?h 2!y ''( t ) ? ? ?hp ?1y p!( p)(t )并用 y k 代替 y ( t k ) ,则得到一个 p 阶近似公式? y k ?1 ? y k ? h ? ( t k , h ) ? ? y 0 ? y (t0 ) k ? 0,1, ? n ? 1(4.1.13)显然 p=1 时,式(4.1.13)就是(4.1.10),即为欧拉方法.但当 p ? 2 时,如果直接 利用公式(4.1.13)求解数值解问题,就需要计算 f ( t , y ) 的高阶微商.这个计算量是 很大的.因此,直接利用式(4.1.13)构造高阶公式是不实用的. 通常,在求解高阶精度的数值解时,我们并不是直接使用 Taylor 级数,而是利 用它的思想,即计算 f ( t , y ) 在不同结点的函数值,然后作这些函数值的线性组合, 构造近似公式,式中有一些可供选择的参数.将近似公式与 Taylor 展开式相比较, 使前面的若干项密合,从而使近似公式达到一定的精度. 若线性组合选取的函数值个数为 n ,计算达到的精度阶数为 p ,这样的 Runge-Kutta 方法就是我们通常所说的 p 阶 n 级 Runge-Kutta 方法. 后来,德国学者 Felhberg 对传统的 Runge-Kutta 方法进行了改进,使原来的定 步长算法变为自动变换步长的 Runge-Kutta-Felhberg 方法,保证了更高的精度和数 值稳定性,它已经成为求解数值解问题中最常见的方法.具体的算法我们在这里不 做介绍,下面我们主要来讨论如何通过 matlab 软件来运用这一方法求解微分方程 的数值解问题. 二、用 matlab 求解微分方程数值解 (一) 一阶微分方程组初值问题 Matlab 下求解一阶微分方程组初值问题(4.1.9)的数值解的最常见的方法 ode45()函数,该函数实现了前面介绍的变步长四阶五级 Runge-Kutta-Felhberg 算 法,可以采用变步长的算法求解微分方程. ode45()函数的调用格式[t, y]= o d e4 5 (o d efu n ,[t 0 ,t f ],y 0 ) [t, y]= o d e4 5 (o d efu n ,[t 0 ,t f ],y 0 ,o p tio n ,p 1 ,p 2 , ? p m ) 常用格式 完整格式参数说明 [t, y] 中 t 表示自变量,是标量.y 表示函数,可以是标量也可以是向量. y 0 表示 初值.当 y 为 n 维向量时,表明求解的问题是有 n 个未知函数的方程组,则 y 0 也应 为 n 维向量. [ t 0 , t f ] 表示微分方程的求解区间,即自变量的取值范围.如果只给出 一个值tf则表示初始时刻为 t 0? 0 .另外,该函数还允许 t 0 ? t f,即可以认为 t 0 为终值时刻, t f 为初始时刻,则 y 0 表示状态变量的终值,而该函数可以直接求解这样 一个终值问题. odefun 表示待解的微分方程(4.1.9)中 f ( t , y ) 所写成的 m-文件的名称,m-文件 的编写格式为 function ydot=odefun(t,y) ydot=f(t,y); 也可以不把函数单独用 m-文件保存,而是直接用 inline()函数输入,从而获得函数 名,编写格式为 odefun=inline(‘f(t,y)’,’t’,’y’); 其中,t 是自变量,即使所求方程是自治系统,也需要写出 t,否则 matlab 在变量传递 过程中将出现问题.y 为未知函数,ydot 表示未知函数的导数. 如果求解的问题是微分方程组问题,即我们之前讨论过的 y 为向量的情况,则 在书写上述 m-文件或 inline()函数时,待解方程 f ( t , y ) 应以 ydot 的分量形式写成, 或者仍用单个方程的编写格式,多个方程用[ ]括起来,并以;隔开,在后面的例题 中我们将具体的看到. options 表示控制选项.在微分方程求解中有时需要对求解算法及控制条件 进行进一步设置,这可以通过求解过程中的 options 变量来进行修改.初始 options 变量可以通过 odeset()函数来获取,该变量是一个结构体变量,其中有众多成员变 量,表 4-1 列出了常用的一些成员变量.表 4-1 微分方程求解函数的控制参数表参数名 RelTol参数说明为相对误差容许上限,默认值为 10 较高的精度,可以适当减小该值. 为绝对误差容许上限,默认值为 10 分量允许的绝对误差.?6?3,在一些特殊的微分方程求解时,为了保证AbsTol MaxStep,它是一个向量,每个分量表示各未知函数为求解方程允许的最大步长 Mass Jacobian微分代数方程中的质量矩阵,用于描述微分代数方程 为描述 Jacobi 矩阵函数的函数名,如果已知 Jacobi 矩阵,可以加速仿真过程修改 option 变量有两种方式: ①用 odeset()函数设置,其格式为 options=odeset(‘RelTol’,1e-7); 达到的目的是,将相对误差上限设置为较小的 10 ? 7 ②直接修改 options 的成员变量,也可以达到这样的目的,其格式为 options= options.RelTol=1e-7;p1 , p 2 , ? , p m 表示附加参数.即在函数表达式中,除了 t , y还有其他的可变的参数,引入这样的附加参数,使得微分方程的某些参数可以选择不同的值,对不同 值求解时,考虑附加参数可以避免每次修改模型文件.但是值得注意的是,这种情 况下的函数定义格式有特殊要求,格式为 function ydot=odefun(t,y,flag,p1,p2,?,pm)ydot=f(t,y,p1,p2,?,pm); 或 odefun=inline(‘f(t,y,p1,p2, ?,pm)’,’t’,’y’,’flag’,’p1’,’p2’,?,’pm’); 注意,这里的 flag 变量不能省略,用于占位. 例 4-5 解微分方程组? x ' ? ? x3 ? y ? 3 ?y'? x? y ? x (0) ? 1, y (0) ? 0.5 ?0 ? t ? 30解: 该方程是非线性微分方程,不能用我们之前讨论的求解解析解的方法来处理, 只能用数值解法求解.首先,我们将待解方程写成 m-文件,用文件名 eg1fun.m 保 存.我们将 2 个未知变量 x,y,以分量的形式写成一个向量变量 x. %函数 eg1fun.m function f=eg1fun(t,x)f(1)=-x(1)^3-x(2); f(2)=x(1)-x(2)^3; f=f(:); %保证 f 为列向量 这个 m-文件还有另一种等价的书写方式. %函数 eg1fun.m 的另一种书写方式 function f=eg1fun(t,x)f=[-x(1)^3-x(2);x(1)-x(2)^3]; 这样我们就完成了对待解方程的描述,然后输入命令 &&[t,x]=ode45('eg1fun',[0 30],[1;0.5]); subplot(1,2,1); plot(t,x(:,1),t,x(:,2),':'); title('函数曲线图'); xlabel('t'); legend('x(t)','y(t)'); axis square subplot(1,2,2); plot(x(:,1),x(:,2)); title('相平面图'); xlabel('x'); ylabel('y'); axis square 运行结果: %求解方程 %分 2 部分绘图,编辑第一部分 %在第一个图中绘制 x,y 的函数图象 %图象名称 %横轴说明 %说明图例 %显示正方形的坐标系 %编辑第二部分图象 %在第二个图中绘制相平面图 %图象名称 %横轴说明 %纵轴说明 %显示正方形的坐标系图 4-1 例 4-5 函数曲线和相平面图 运行程序,得到函数图象和相平面图.根据相平面上的相轨线,可以看出解从初值 开始的变化趋势,这对于研究微分方程平衡点的稳定性是很有意义的,在后面的章 节中我们将具体讨论. 例 4-6 解微分方程组y1 ' ? y 2 y 3 ? ? y 2 ' ? ? y1 y 3 ? ? y 3 ' ? ? 0 . 51 y 1 y 2 ? ? y 1 ( 0 ) ? 0 , y 2 ( 0 ) ? 1, y 3 ( 0 ) ? 1 ?解:此问题仍为非线性微分方程组,我们求它的数值解,不妨将求解区间取为t ? [0, 2 0 ] .下面我们用 inline()函数来描述方程组,这样就不用单独编写 m-文件来描述方程. 输入命令: &&eg2fun=inline('[y(2)*y(3);-y(1)*y(3);-0.51*y(1)*y(2)]','t','y'); tf=20; y0=[0;1;1]; [t,y]=ode45(eg2fun,[0,tf],y0); plot(t,y(:,1),'-',t,y(:,2),'*',t,y(:,3),'+'); %绘出解的图象 legend('y1','y2','y3'); %说明图例 运行得到函数曲线图图4-2 例4-6函数曲线图例4-7 编写带有附加参数的matlab函数,来描述食饵-捕食者模型? x '1 ( t ) ? x1 ( a 1 ? b1 x 2 ) ? ? x ' 2 ( t ) ? ? x 2 ( a 2 ? b 2 x1 ) ? x (0 ) ? x (0 ) ? 1 2 ? 1 并对 a1 ? 3, b1 ? 2, a 2 ? 2.5, b 2 ? 1 和 a1 ? 2, b1 ? 1, a 2 ? 1, b 2 ? 0 .5 分别求解,画出解的曲 线图和相轨线图. 解:选定附加参数为 a1 , b1 , a 2 , b 2 ,编写如下 m-文件来描述题述方程 %函数 eg3fun.m function xdot=eg3fun(t,x,flag,a1,b1,a2,b2) xdot=[x(1)*(a1-b1*x(2));-x(2)*(a2-b2*x(1))]; 注意,这里的 flag 变量不能省略.这样我们就完成了对方程的描述.我们不妨 取自变量区间为[0,10],对方程求数值解. 输入命令 &&tf=10;x0=[1;1]; a1=3;b1=2;a2=2.5;b2=1; [t,x]=ode45(‘eg3fun’,[0,tf],x0,[],a1,b1,a2,b2); subplot(1,2,1);plot(t,x(:,1),t,x(:,2),':');title('函数曲线图'); xlabel('t');legend('x1','x2'); subplot(1,2,2);plot(x(:,1),x(:,2));title('相平面图'); xlabel('x1');ylabel('x2'); 运行结果图 4-3a1 ? 3, b1 ? 2, a 2 ? 2.5, b 2 ? 1当 a1 ? 2, b1 ? 1, a 2 ? 1, b 2 ? 0 .5 时,只需要修改上述命令中的给 a1,b1,a2,b2 的赋 值语句,就可以得到需要的结果. 图 4-4a1 ? 2, b1 ? 1, a 2 ? 1, b 2 ? 0.5(二) 一阶微分方程组边值问题 对于一阶常微分方程组的边值问题? y ' ? f (t , y ) ? ? g ( y ( a ), y ( b )) ? 0 a ?t ?b(4.1.14)这里 y,f,g 都可以是向量. 用 matlab 求(4.1.14)的数值解的基本格式 sinit=bvpinit(tint,yinit); 由在粗略结点 tinit 的预估解 yinit 生成粗略解网络 sinit. sol=bvp4c(odefun,bcfun,sinit) 其中 odefun 是微分方程组函数,bcfun 为边值条件函数,sinit 是由 bvpinit 得到的粗 略解网络.求得的边值问题的解 sol 是一个结构体变量,sol.x 为求解结点,sol.y 为数 值解. sx=deval(sol,ti) 计算由 bvp4c 得到的解在 ti 的值. 例 4-8 求解边值问题? y '1 ? y 2 ? ? y '2 ? ? y1 ? ? y 1 ( a ) ? 0, y 1 ( b ) ? 2 ? 0a ?t ?b解:输入命令 &&sinit=bvpinit(0:4,[1;0]); odefun=inline('[y(2);-abs(y(1))]','t','y'); bcfun=inline('[ya(1);yb(1)+2]','ya','yb'); sol=bvp4c(odefun,bcfun,sinit); t=linspace(0,4,101); y=deval(sol,t); plot(t,y(1,:),sol.x,sol.y(1,:),'o',sinit.x,sinit.y(1,:),'x'); legend('解曲线','求解点','粗略解'); 运行结果图 4-5 例 4-8 函数曲线图(三) 高阶微分方程 由于 matlab 只能处理一阶微分方程组问题,对于高阶微分方程的初边值问题, 我们必须通过变量代换,使它转换为一阶方程组,才能用 matlab 求解数值解.我们 考虑高阶常微分方程的初值问题y(n)? f ( t , y , y ', ? , y( n ? 1))(4.1.15)且初始值 y (0), y '(0), ? y ( n ?1) (0) 为已知. 令 x1 ? y , x 2 ? y ', ? , x n ? y ( n ?1) ,通过这样的变量代换,高阶微分方程(4.1.15)就 转化为一阶方程组的形式 ? x '1 ? ? x '2 ? ? ?x ' ? n? x2 ? x3 ? ? f ( t , x1 , x 2 , ? , x n )(4.1.16)且初值条件为 x1 (0) ? y (0), x 2 (0) ? y '(0), ? , x n (0) ? y ( n ?1) (0) .这样,我们就可以通过 之前讨论过的一阶微分方程组的数值解法来研究这一问题. 对于高阶微分方程组,只要依次将每个未知变量的各阶导数均单独定义变量, 就可以同样的转化为一阶微分方程组. 例 4-9 用数值解法求解微分方程? y '' ? ( y 2 ? 1) y ' ? y ? 0 ? ? y (0 ) ? ? 0 .2, y '(0 ) ? ? 0 .7解:首先我们需要将这个 2 阶微分方程转化为一阶微分方程组,我们令x1 ? y, x 2 ? y ' ,则原方程转换为? x '1 ? x 2 ? 2 ? x ' 2 ? ? ( x1 ? 1) x 2 ? x1 ? x (0 ) ? ? 0 .2 , x (0 ) ? ? 0 .7 2 ? 1对[0,20]区间求数值解,输入命令 &&tf=20;x0=[-0.2;-0.7]; eg5fun=inline('[x(2);-((x(1))^2-1)*x(2)-x(1)]','t','x'); [t,y]=ode45(eg5fun,[0,tf],x0); plot(t,y); legend('y','y'''); 运行得到函数图象 图 4-6 例 4-9 函数曲线图例 4-10 求解微分方程组(竖直加热板的自然对流问题)2 2 ?d3 f ? df ? d f ?3f ? 2? ? ? ?T ? 0 3 2 d? ? d? ? ? d? ? 2 dT ?d T ? 2 .1 f ?0 ? 2 d? ? d? 2 ? df d f dT (0 ) ? 0, (0 ) ? 0 .6 8, T (0 ) ? 1, (0 ) ? ? 0 .5 ? f (0 ) ? 0, 2 d? d? d? ? ?解:这是一个高阶方程组问题,我们需要通过变量代换转化为一阶微分方程组,从 而寻求数值解.我们令y1 ? f , y 2 ? df d? , y3 ? d f d?2 2, y4 ? T , y5 ?dT d?原方程转化为如下一阶方程组dy ? d y1 ? y2 , 2 ? y3 ? d? d? ? ? dy3 2 ? ? 3 y1 y 3 ? 2 y 2 ? y 4 ? d? ? ? dy dy ? 4 ? y 5 , 5 ? ? 2 .1 y 1 y 5 d? ? d? ? y (0 ) ? 0, y (0 ) ? 0, y (0 ) ? 0 .6 8, y (0 ) ? 1, y (0 ) ? ? 0 .5 2 3 4 5 ? 1输入命令 &&y0=[0;0;0.68;1;-0.5];tf=5; eg6fun=inline('[y(2);y(3);-3*y(1)*y(3)+2*y(2)^2-y(4);y(5);-2.1*y(1)*y(5)]','t','y'); [t,y]=ode45(eg6fun,[0,tf],y0); plot(t,y(:,1),t,y(:,4),':'); legend('f','T'); 运行得到函数图象图 4-7 例 4-10 函数曲线图§4.1.3 导弹追踪问题一、导弹追击问题 一艘敌舰在某海域内沿正北方向航行时,我方战舰恰位于敌舰的正西方向 1km 处.我舰向敌舰发射导弹,导弹头始终对准敌舰,敌舰速度为 v 0 =1km/min,导弹 速度为敌舰速度的 5 倍.问需要多长时间,在何处导弹击中敌舰? 以我舰位置为坐标原点,正北方向为 y 轴方向建立坐标系,设 t 时刻导弹所处 的位置为 P ( x ( t ), y ( t )) ,敌舰所处的位置为 Q (1, v 0 t ) .图 4-8 导弹追击问题由于导弹头始终对准敌舰,因此直线 PQ 是导弹运行轨迹 OP 在 P 点的切线, 即dy dx ? v0t ? y 1? x dy dx从而v 0 t ? (1 ? x ) ? y(4.1.17)又因为导弹速度为敌舰速度的 5 倍,所以 OP 弧的长度为 AQ 长度的 5 倍. 即?联立(4.1.17)和(4.1.18)得到x 0? dy ? 1? ? ? d x ? 5v0t ? dx ?2(4.1.18)?x 0dy ? dy ? 1? ? ? 5y ? d x ? 5(1 ? x ) dx ? dx ?2两边对 x 求导,得到微分方程? d y? 1? ? 1 ? ? 5 (? x ? d x?2d y ) 2 d x2(4.1.19)该问题的初始条件为y (0 ) ? 0, y '(0 ) ? 0首先,我们将二阶方程化为一阶方程组,令 y1 ? y , y 2 ?? d y1 ? dx ? y2 ? 1 / 2 ? dy ?1 ? y 22 ? ? 2 ? ? dx 5(1 ? x ) ? ? y (0 ) ? 0, y (0 ) ? 0 1 2 ? ? ?dy dx,(4.1.19)转化为(4.1.20)对于这一问题,我们可以求得方程的解析解,从而确定导弹击中敌舰位置和时 间,也可以通过求数值解的方法得到同样的结论. 1) 解析解. 关于 y 2 的方程1/ 2 ? ?1 ? y 22 ? dy2 ? ? ? dx 5(1 ? x ) ? ? y 2 (0 ) ? 0分离变量,有 dy2?1 ? y ?2 21/ 2?dx 5(1 ? x )两边积分,并代入初始条件,得ln ( y 2 ? 1 ? y2 ) ? ?21 5ln (1 ? x )即y2 ? 1 ? y 2 ? (1 ? x )2 ? 1 5(4.1.21)1又因为y2 ? 1 ? y2 ?2?1 y2 ? 1 ? y22? ? (1 ? x ) 5(4.1.22)(4.1.21)和(4.1.22)两式相加,得到y2 ? 1 2 ((1 ? x )? 1 5 1? (1 ? x ) 5 )下面,我们求解关于 y1 的方程1 1 ? ? d y1 1 ? y 2 ? ((1 ? x ) 5 ? (1 ? x ) 5 ) ? 2 ? dx ? y (0 ) ? 0 ? 1直接积分,并代入初始条件得到y ? y1 ? ? 5 84(1 ? x ) 5 ?5 126(1 ? x ) 5 ?5 24当 x ? 1 时, y ?t ? y v0 ? 1 2 .5 s5 24,即当敌舰航行到点(1,5 24)处时,被导弹击中.被击中的时间为.2) 数值解. 首先,编写出描述方程(4.1.20)的 m-文件 function dy=equ20(x,y) dy=[y(2);1/5*sqrt(1+y(2)^2)/(1-x)]; 由于方程在 x ? 1 处没有定义,我们取求解区间为[0,1-1e-8],输入命令 xf=1-1e-8; [x,y]=ode45('equ20',[0,xf],[0;0]) plot(x,y(:,1)); hold on y=0:0.01:2; plot(1,y,'*'); 运行结果 y= 0 0.0000 ?? 0.3 0 0.0001 ?? 16.8图 4-9 导弹击中曲线说明当 x ? 1 时, y ? 0.2083 ,导弹击中敌舰的位置大致在(1,0.2083)处,击中的 时间 t ?y v0=0.2083min=12.498s.这与求解析解得到的结果是一致的.比较上面 2 种求解方法,我们发现,解析解可以得到解的真实值,但是方法比较 复杂,技巧性强,只能处理比较简单的方程问题;而数值解虽然得到的解与真实值 之间会有一些允许范围内的误差,但是方法比较容易掌握,更具一般性. 二、导弹系统的改进 现根据情报,这种敌舰能在我舰发射导弹后 T 分钟做出反应并摧毁导弹.因 此,我们要求改进电子导弹系统,使其根据敌舰与我舰的距离,行使方向和速度, 能自动判断出敌舰是否在有效打击范围之内. 对于更一般的情况,我们做出如下假设,设敌舰在我舰正东方向 d km 处,行驶 速度为 v 0 km/min,行驶方向与正东方向的夹角为 ? ,导弹的飞行速度为 v km/min. 问题的关键是计算出导弹击中敌舰所需要的时间 t * ,并将 t * 与 T 比较,若 t * &T,则敌 舰在打击范围内. 我们仍以我舰位置为坐标原点,以正北方向为 y 轴建立坐标系,设 t 时刻导弹 所处的位置为 P ( x ( t ), y ( t )) ,敌舰所处位置为 Q ( d ? v 0 t co s ? , v 0 t sin ? ) .图 4-10 导弹击中问题由于导弹头始终对准敌舰,因此直线 PQ 是导弹运行轨迹 OP 在 P 点的切线, 即dy dx v v0x 0?v 0 t sin ? ? y d ? v 0 t co s ? ? x v v0(4.1.23)又因为导弹速度为敌舰速度的倍,所以 OP 弧的长度为 AQ 长度的2倍.即?? d y? 1? ? ? d x? ? d x?v0v0v t(4.1.24)联立方程(4.1.23)和(4.1.24),消去 v 0 t ,再对方程两边对 x 求导,得到微分方程d y dx2 2d y 2 1 / 2 d y 2 ( 1? ( ) ) (? ? n s i ? co s ) v dx dx ? dy dy (d ? x ) ( s ? n i ? ? o?s ) ? c o s ? (x ? (y c d dx dxv0(4.1.25)) )同样,我们令 y1 ? y , y 2 ?dy dx,(4.1.25)转化为一阶微分方程组 ? d y1 ? dx ? y2 ? v0 2 1 / 2 2 ? (1 ? ( y 2 ) ) (sin ? ? y 2 co s ? ) ? dy2 v ? ? dx ( d ? x )(sin ? ? y 2 co s ? ) ? co s ? ( y 2 ( d ? x ) ? y 1 ) ? ? y (0 ) ? 0, y (0 ) ? 0 1 2 ? ? ?(4.1.26)从图象上来判断,我们知道,当 P 和 Q 两点的运动曲线相遇时,导弹击中敌舰. 因此我们可以认为,若 x * , y * 满足方程(4.1.25)且y ? ( x ? d ) tan ?* *(4.1.27)则 点 ( x* , y* ) 为 导 弹 击 中 敌 舰 的 击 中 点 . 再 根 据 Q 点 的 表 达 式Q ( d ? v 0 t co s ? , v 0 t sin ? ) ,可以计算出击中时间 t ?*y*v 0 sin ?.若 t * &T,则敌舰在打击范围内,可以发射. 从上述分析,我们知道,当 v 0 、v、d、 ? 、T 已知时,根据(4.1.25)和(4.1.27)可 以计算出导弹击中敌舰所需要的时间 t * ,且当 t * &T 时,敌舰在打击范围内.这种对 于每组 v 0 、v、d、? 、T 分别计算从而作出判断的方法,称为在线算法.如果在敌舰 和我舰的情况已知,即 v 0 、v、T 已知,计算出所有在打击范围内的 d 和 ? ,之后要判 断敌舰是否在打击范围,只需要在计算结果中查询,这种算法称为离线算法.从原 理上来看,我们发现,在线算法灵活,容易调整参数和模型,但是速度比较慢;离线算 法事先计算好,实时使用查询方式,不需要计算,速度快. (1)在线算法 现已知 v 0 =90km/h,v=450km/h,d=50km, ? ??4,T=0.1h.首先,我们编写带参数的 m-文件描述方程(4.1.26).为了使方程始终有意义,我们通常在分母上加上一个 无穷小量 1e-8. %方程(4.1.26) function dy=equ26(x,y,v0,v,d,theta) dy(1)=y(2); dy(2)=((v0/v)*(1+(y(2))^2)^(1/2)*(sin(theta)-y(2)*cos(theta))^2) /((d-x)*(sin(theta)-y(2)*cos(theta))+cos(theta)*(y(2)*(d-x)+y(1))+1e-8); dy=dy(:); 输入命令 v0=90;v=450;d=50;theta=pi/4;T=0.1; %初始化参数 [x,y]=ode45(@equ26,[0,60],[0;0],[],v0,v,d,theta); %计算导弹运行轨迹方程的数值解 z=(x-d)*tan(theta); %计算敌舰运行的轨迹 n=length(x); for i=1:n if abs(z(i)-y(i,1))&1e-6 xk=x(i); yk=y(i,1); end end xk yk %求出击中点xk,yk坐标 for i=1:n if z(i)&0 z1(i)=z(i); else z1(i)=0; end end plot(x,y(:,1),x,z1(:),xk,yk,'*'); legend('导弹运行轨迹','敌舰运行轨迹','击中点'); tk=yk/(v0*sin(theta)+1e-8) %计算击中时间 if tk&T disp(['敌舰在打击范围内，击中地点在', num2str(xk),num2str(yk),'击中时间为', num2str(tk)]); else disp(['敌舰不在打击范围内']); end 得到结果 击中点 xk=58.4056,yk=8.4056,击中时间 tk =0.1321,所以敌舰不在打击范围内,应 等接近一些再发射. 图4-11 在线算法(2)离线算法 在介绍离线算法之前,我们先来看方程(4.1.26)的另一种表达.因为导弹的线 速度v ? ? dx ? ? dy ? ? ? ?? ? ? dt ? ? dt ?2 2(4.1.28)方程(4.1.28)与(4.1.23)联立,可以建立以 t 为参数的关于 x,y 的参数方程? dx ? ? dt ? ? ? ? ? dy ? ? dt ? ? ? v ? dy ? 1? ? ? ? dx ? v ? dx ? 1? ? ? ? dy ?2 2?v ? v 0 t sin ? ? y ? 1? ? ? ? d ? v 0 t co s ? ? x ?2?v ? d ? v 0 t co s ? ? x ? 1? ? ? ? v 0 t sin ? ? y ?2(4.1.29)当 x ( t ) 满足 x ( t ) ? d ? v0 t c o s? ,则导弹击中了敌舰.由于我们所要求的打击范围是 在 t ? T 中讨论的,所以考虑以 t 为参数的方程(4.1.29)能使求解过程大大简化. 现已知 v 0 =90km/h,v=450km/h,T=0.1h,要计算出所有在打击范围内的 d 和 ? . 依题意, 0 ? ? ? ? .据此,我们可以确定 d 的取值范围.当 ? ? 0 时,敌舰正好背向行驶, 导弹直线运行,击中时间t ? d /( v ? v 0 ) ? T求得 d m in ? T ( v ? v 0 ) ? 36 km .当 ? ? ? 时,敌舰迎面驶来,导弹直线运行,击中时间 t ? d /( v ? v 0 ) ? T则 d m ax ? T ( v ? v 0 ) ? 54 km .所以 3 6 ? d ? 5 4 .这样,我们对于所有可能的 d 和 ? 的 取值,计算击中所需时间,从而对不同的 ? ,得到 d 的临界值.具体应用时直接查 询判断. 编写如下 m-文件描述方程(4.1.29) %方程(4.1.29) function dy=equ29(t,y,v0,v,d,theta) dydx=(v0*t*sin(theta)-y(2))/(d+v0*t*cos(theta)-y(1)+1e-8); dy(1)=v/(1+dydx^2)^(1/2); dy(2)=v/(1+dydx^(-2))^(1/2); dy=dy(:); 输入命令 v0=90;v=450;d=50;theta=pi/4;T=0.1; i=1; for d=54:-1:36 for theta=0:0.1:pi [t,y]=ode45(@equ29,[0,T],[0;0],[],v0,v,d,theta); if max(y(:,1)-d-v0*t*cos(theta))&0 range(i,:)=[d,theta]; i=i+1; plot(range(:,1),range(:,2)); xlabel('d');ylabel('theta'); 运行得到临界曲线,即在该曲线上方的 d 和 ? 值,所对应的敌舰位置在打击范 围内,曲线下方不在打击范围内. 图4-12 离线算法因此 d=50km, ? ??4,不在打击范围内.§4.1.4 微分方程稳定性理论简介在处理实际问题时,对于有些微分方程模型我们不仅要得到问题的解,有时 还需要研究解的稳定性,即解对初始值的连续依赖性,如果解在一定范围内是稳 定的,那么初始条件发生一些小的扰动(如实验测量误差等),对问题的解不会造 成影响.另外,还有一些问题我们并不需要求解,而通过解的变化趋势的研究,并 分析一些特殊解的稳定性就可以解决问题.本节仅介绍几类特殊,但常用的常微 分方程稳定性分析的基本方法. 一、单个常微分方程的平衡点及稳定性 若微分方程dx dt ? f (x)(4.1.30)方程右端不显含自变量 t ,称为自治方程.代数方程的实根 x ? x 0 称为方程(4.1.30) 的平衡点(或奇点).注意到,平衡点也是方程的解(奇解). 如果从一定范围内的初始条件出发,方程(4.1.30)的解 x (t ) 都满足t ? ??lim x ( t ) ? x 0 ,则称平衡点 x 0 是稳定的. 对于一些不易求解的问题,我们可以不求方程(4.1.30)的解,不用定义来判断 平衡点 x 0 的稳定性,下面我们来介绍这种方法. 将 f ( x ) 在 x 0 点作泰勒(Taylor)展开,只取一次项,得到方程(4.1.30)的近似线 性方程dx dtx0? f ' ( x 0 )( x ? x 0 ),0(4.1.31))t|也是方程(4.1.31)的平衡点,方程(4.1.31)的通解为 x ( t ) ? ce | f '( x? x 0 .关于 x 0点稳定性有如下结论: ① 若 f ' ( x 0 ) ? 0 ,则 x 0 对于方程(4.1.31)和(4.1.30)都是稳定的; ② 若 f ' ( x 0 ) ? 0 ,则 x 0 对于方程(4.1.31)和(4.1.30)都是不稳定的. 二、二阶常微分方程组的平衡点及稳定性 现在讨论二阶微分方程组 ? dx ? f ( x , y ), ? dt ? ? ? d y ? g ( x , y ). ? dt ?(4.1.32)它的解x ? x ( t ), y ? y ( t )在以 t , x , y 为坐标的欧氏空间中决定了一条曲线,如果把时间 t 看作参数,仅考虑x, y为坐标的欧氏空间,此空间称为方程(4.1.32)的相平面(若方程组是高阶,则称为相空间).对于右端函数不显含时间 t 的自治系统? dx ? dt ? f ( x , y ), ? dy ? ? g ( x , y ). ? dt(4.1.33)代数方程组? f ( x, y ) ? 0, ? ? g ( x, y ) ? 0.的 实 根 x ? x 0 , y ? y 0 称 为 方 程 (4.1.33) 的 平 衡 点 ,记 作 P0 ( x 0 , y 0 ) . 它 也 是 方 程 (4.1.33)的解. 如果从一定的范围内的初始条件出发,方程(4.1.33)的解 x ( t ), y ( t ) 都满足t ? ??lim x ( t ) ? x 0 , lim y ( t ) ? y 0 ,t ? ??则称平衡点 P0 是稳定的,否则称 P0 是不稳定的. 与单个方程的讨论类似,对于不易求解的微分方程组,我们可以通过研究与 其对应的近似线性方程组,从而得到平衡点和稳定性.在这里,我们省略证明过程, 只给出判别平衡点 P0 是否稳定的判别准则.令? f ( P0 ) ? ? f ( P0 ) ? g ( P0 ) ? ?x p ? ?? ? ?, q ? ?g ( P ) ?y ? 0 ? ?x ?x ? f ( P0 ) ?y , ? g ( P0 ) ?y关于 P0 点稳定性有如下结论: ① 当 p ? 0 且 q ? 0 时,平衡点 P0 是稳定的; ② 当 p ? 0 或 q ? 0 时,平衡点 P0 是不稳定的. §4.1.5 动物种群模型分析一、单一种群模型 (一) Malthus 增长模型 英国经济学家和人口统计学家 Malthus()根据百余年的统计资料, 在 1798 年提出了闻名于世的人口指数增长模型(Malthus 人口模型).此模型也适用 于单种族问题.此模型有两个基本假设前提： ① 假设人口数量 x (t )是时间 t 的连续可微函数, 且 x (0) = x0. 虽然人口数量 总是离散变化的, 但实际上增加一个人所引起的变化与庞大的人口数量相比是 微不足道的. 这里实际上是将离散变量作连续化处理, 在差分方程中则恰好相 反. ② 人口数量的增长速度与现有人口数量成正比, 比例系数为 r. 这个假设在 当时是比较合理的. 在此基础上, Malthus 提出了如下的人口模型(即指数增长模型)：dx dt= rx,x (0) = x0 .不难求得其解为 x (t ) = x0ert, 这个模型比较准确地反映了在
年期间的 人口总数. (二)Logistic 增长模型 因为 r＞0, x0＞0, x (t ) = x0ert 表明人口数量 x (t ) 将随时间 t 按指数规律无限 地增长, 对于有限的地球资源来说, 这当然是不合理的. 即假设②在自然环境条 件较好, 人口密度较小时是合理的, 而对于人口密度较大时它是不合理的, 这是 因为人类生存空间及可利用资源(食物、水、空气)等环境影响对人口的增长起着 阻滞作用. 从这里我们可以看出, 在进行数学建模时, 我们必须作合理的假设. 现在我们对指数增长模型加以改进.设人类生存空间及可利用资源(食物、 水、 空气)等环境因素所能容纳的最大人口容量为 K (称为饱和系数). 人口数量 x 的增 长速率不仅与现有人口数量成正比, 而且还与人口尚未实现的部分(相对最大容 量 K 而言)所占比例(K ??x)/K 成正比, 比例系数为固有增长率 r. 修改后的模型为? dx ?K ? x? ? rx ? ?, ? ? K ? ? dt ? x (0) ? x . 0 ?Kx不难求得其解为 x ( t ) ?0x 0 ? ( K ? x 0 ) e ? rt.这就是非常著名的 Logistic 增长模型(即阻滞增长模型), 它是荷兰数学家、生物学家 Verhulst 在 1839 年首次提出的.因子(K ??x)/K 的生物学含义是“剩余空间”或尚未利用的增长机会. 图 4-13 描绘了dx dt与 x 之间的关系, 它表明人口变化率dx dt随着人口数量 x 的 增加而先增后减, 在 x = K / 2 处到达最大值, 在 x = K 处, Logistic 曲线 x (t ), 且有 lim 义看这是明显的结果.t? ?dx dt= 0. 图 4-14 描绘了x (t ) ? K, 即 x = K 是稳定的平衡点, 从模型本身的意图 4 -13图 4 -14Logistic 模型的曲线Logistic 模型用途十分广泛, 除了用于预测人口增长外, 也可完全类似地用 于虫口增长、疾病的传播、谣言的传播、技术革新的推广、销售预测等.但这个 模型的最大缺点是饱和系数 K 不易确定. 二、两种物种相互竞争与相互依存模型当某个自然环境中只有一种生物的群体(生态学上称为种群)生存时, 人们常用 Logistic 模型来描述这个种群数量的演变过程.即dx x ? ? ? rx ? 1 ? ? dt K ? ?x(t ) 是种群在时刻 t 的数量, r 是固有增长率, K 是环境资源容许的种群最大数量, 且 x = K 是 稳定的平衡点. 如果一个自然环境中有两个或两个以上种群生存, 那么它们之间就要存在着或是相互 竞争, 或是相互依存, 或是弱肉强食(捕食--被捕食模型)的关系. ⑴ 种群的相互竞争 当两个种群为了争夺有限的同一种食物来源和生活空间而进行生存竞争时, 最常见的 结局是竞争力较弱的种群灭绝, 竞争力较强的种群达到环境容许的最大数量. 设有甲乙两个种群, 当他们独自在一个自然环境中生存时, 数量的演变均遵从 Logistic 规律, x1 (t ), x2 (t )是两个种群的数量, r1, r2 是它们固有增长率, N1, N2 是它们的最大容量. 于 是对种群甲有d x1 ? x1 ? r1 x 1 ? 1 ? ? dt N1 ? ? ? ? ?, ? ?x1 N1其中因子 ? 1 ? ?x1 ? ? 反映由于甲对有限资源的消耗导致对它本身增长的阻滞作用, ? N1 ?可解释为相对于 N1 而言单位数量的甲消耗的供养甲的食物量(设食物总量为 1). 当两个种群在同一自然环境中生存时, 考察由于乙消耗同一种有限资源对甲的增长产 生的影响, 可以合理地在因子 ? 1 ? 1 ? 中再减去一项, 该项与种群乙的数量 x2 (相对于 N2 而 ? ? N1 ? ? 言)成正比, 得到种群甲的增长方程d x1 ? x1 x2 ? r1 x 1 ? 1 ? ??1 ? dt N1 N2 ? ? ?, ? ??x?这里的?1 意义是, 单位数量乙(相对 N2 而言)消耗的供给甲的食物量为单位数量甲(相对 N1) 消耗的供给甲的食物量的?1 倍. 类似地, 得到种群乙的增长方程dx 2 ? ? r2 x 2 ? 1 ? ? ? dt ? x12?N1x2 ? ?, ? N2 ?对?2 可作相应地解释. 关于种群的相互竞争模型? d x1 ? x1 x2 ? ?, ? r1 x 1 ? 1 ? ??1 ? ? ? N1 N2 ? ? dt ? ? ? x1 x2 ? ? dx 2 ?. ? r2 x 2 ? 1 ? ? 2 ? ? ? ? dt N1 N2 ? ? ?的稳定性分析见表 4-1. 表 4-2 种群相互竞争模型的平衡点及其稳定性 平衡点 P1 (N1, 0) P2 (0, N2 ) P3 ( N1 (1???1) N2 (1???2) , ) 1???1?2 1???1?2 P4 (0, 0) p r1 ??r2 (1???2) r2 ??r1 (1???1) r1 (1???1) + r2 (1???2) 1???1??2 ??(r1 + r2) q ??r1 r2 (1???2) ??r1 r2 (1???1) r1 r2 (1???1) (1???2) 1???1??2 r1 r2 稳定条件?1＜1,????2＞1 ?1＞1,????2＜1 ?1＜1,????2＜1不稳定根据表 4-2 和?1、??2 的含义, 稳定的平衡点 P1, P2, P3 在生态学上的意义解释如下： ① 当?1＜1,????2＞1 时, ?1＜1 意味着在对供养甲的资源的竞争中乙弱于甲, ??2＞1 意味 着在对供养乙的资源的竞争中甲强于乙, 于是种群乙终将灭绝, 种群甲趋向最大容量, 即 x1(t)、x2 (t)趋向平衡点 P1 (N1, 0). ② 当?1＞1,????2＜1 时, 与①正好相反. ③ 当?1＜1,????2＜1 时, 因为在竞争甲的资源中乙较弱, 在竞争乙的资源中甲较弱, 于 是可以达到一个双方共存的稳定的平衡状态 P3. ④ 当?1＞1,????2＞1 时, P1, P2, P3 都是稳定的平衡点, 需根据初始条件确定, 这是种群竞 争中很少出现的情况, 因为在竞争甲的资源中乙较强, 在竞争乙的资源中甲较强. ⑵ 种群的相互依存模型 自然界中处于同一环境下两个种群相互依存而共生的现象是很普遍的, 其类型可分为 三种：①其中一种能独立生存, 另其中一种则不能独立生存；②两个种群都能独立生存；③ 两个种群都不能独立生存.下面介绍类型①, 另两种留作习题 7. 设种群甲能独立存在, 按 Logistic 规律增长, 种群乙为甲提供食物, 有助于甲的增长.例 如, 植物可以独立生存, 昆虫的授粉作用又可以提高植物的增长率, 而以花粉为食物的昆虫 却不能离开植物单独存活.种群甲的数量演变规律为d x1 ? x1 x2 ? r1 x 1 ? 1 ? ??1 ? dt N1 N2 ? ? ?, ? ?这里 x1 (t), x2 (t), r1, N1, N2 的意义同前, ?1 表示单位数量乙(相对 N2 而言)提供的供养甲的食物 量为单位数量甲(相对 N1)消耗的供给甲的食物量的?1 倍. 种群乙没有甲的存在会灭亡, 设其死亡率为 r2, 甲为乙提供食物对乙的增长有促进作用, 同时乙的增长又会受到自身的阻滞作用, 因此种群乙的数量演变规律为dx 2 ? ? r2 x 2 ? ? 1 ? ? ? dt ? x12?N1x2 ? ?, ? N2 ?这里的?2 表示单位数量甲(相对 N1 而言)提供的供养乙的食物量为单位数量乙(相对 N2)消耗 的供给乙的食物量的?2 倍. 对于种群的相互依存模型? d x1 ? x1 x2 ? ?, ? r1 x 1 ? 1 ? ??1 ? ? ? N1 N2 ? ? dt ? ? ? x1 x2 ? dx 2 ? r2 x 2 ? ? 1 ? ? 2 ? ? ? dt N1 N2 ? ?? ?. ? ?的稳定性分析见表 4-3. 表 4-3 种群相互依存模型的平衡点及其稳定性 平衡点 P1 (N1, 0) P3 ( N1 (1???1) N2 (?2 ??1) , ) 1???1?2 1???1?2 P4 (0, 0) p r1 ??r2 (?2 ??1) r1 (1???1) + r2 (?2 ??1) 1???1??2 ??r1 + r2 q ??r1 r2(?2 ??1) r1 r2 (1???1) (?2 ??1) 1???1??2 r1 r2 稳定条件?2＜1, ?1?2＞1 ?2＞1,???1?2＜1不稳定根据表 4-3 和?1、?2 的含义, 稳定的平衡点 P1, P2 在生态学上的意义解释如下： ① 当?2＜1,???1?2＜1 时, ?2＜1 意味着甲供养乙的资源不足以维持乙的生存, 于是种群 乙终将灭绝, 种群甲趋向最大容量, 即 x1 (t)、x2 (t)趋向平衡点 P1 (N1, 0). ② 当?2＞1,???1?2＜1 时, ?2＞1 意味着甲供养乙的资源能足以维持乙的生存, ?1?2＜1 即 有?1＜1, 意味着对乙供养甲的资源加以限制, 防止甲过分地增长, 于是可以达到一个双方 共存的稳定的平衡状态 P2. ③ 当?2＜1,???1?2＞1 时, 基本上同①. ④ 当?2＞1,???1?2＞1 时, 需要适当对它们加以限制. ⑶ 种群的弱肉强食模型 处于同一自然环境中两个种群之间的关系出了相互竞争和相互依存之外, 还有一种更 为有趣的生存方式： 种群甲靠丰富的天然资源生长, 而种群乙则靠掠食甲为生.本世纪二十年 代, 意大利著名生物学家 D’Ancona 在研究时无意中发现, 第一次世界大战期间从地中海捕 获的鱼中, 鲨鱼等食肉鱼所占的比重十分明显地上升了. D’Ancona 认为, 这一现象决非偶然, 它是由战争期间捕鱼量减少造成的.可是, 为什么捕鱼量减少, 食用鱼反而减少呢？他无法 解释这个现象, 于是就找了著名数学家 Vloterra, 希望他能建立一个数学模型来解释之. Vloterra 根据建模目的将鱼分为两大类：食用鱼 (被捕食者) 和食肉鱼(捕食者), 用 x (t) 和 y (t )分别表示 t 时刻两者的数量.他分析, 假设如果没有食肉鱼, 食用鱼以净自然增长率 r ? (r＞0)的速度无限增长, 即 x /x = r. 但是, 食肉鱼以食用鱼为食物, 致使鱼的增长率降低, 设 ? 降低的程度与食用鱼数量 y (t )成正比. 于是, x / x = r ????y. 常数? 反映食肉鱼掠取食物的 ? 能力. 如果没有食用鱼, 食肉鱼无法生存. 设食肉鱼的自然死亡率为 d (d＞0), 则 y / y= ??d. 食用鱼为食肉鱼提供了食物, 致使食肉鱼的死亡率降低, 或者说为食肉鱼提供了增长的条件. ? 设增长率与食用鱼的数量 x (t )成正比, 于是, y / y = ?d + ? x, 常数 ??反映食用鱼对食肉鱼 的供养能力. 根据上述分析可得该模型为? dx ? rx ? ? xy , ? ? dt ? ? d y ? ? dy ? ? xy . ? dt ?4-4这就是著名的捕食--被捕食模型. 容易看出它的两个平衡点是 A (0, 0) 和 B (x0, y0), 其中 x0 = d /??, y0 = r/ ? . 对于平衡点 A(0, 0), p = r ??d, q = ??rd＜0,所以 A (0, 0)不稳定；对于平衡点 B (x0, y0), p = 0, q = rd＞0, 处 于临界状态, 不能用平衡点稳定性判别准则判别. 下面通过求出该模型的解来进行分析. 从微分方程组(4-1)中消去 d t 后, 用分离变量法得到其通解为 xde???x yre???y = c 4-5 式中的 c＞0, 为任意常数, 可由初始条件确定. 下面再来研究方程(4-5)的图形(通常称之为 微分方程组(4-4)的轨线或相轨线). 记 f (x) = xde???x, g (y) = yre?? y, 利用微积分中函数作图的方法作出它们的图形 (图 4-3, 图 4-4), 函数的极值点分别为 x = x0, y = y0, 函数的极大值分别记为 f m 和 gm.图 4 -3图 4 -4当 c＞f m gm 时, 方程(4-5)无解； 当 c = f m gm 时, 方程(4-5)有唯一解：x = x0, y = y0,轨线退化为一点 B (x0, y0) (平衡点)； 当 c＜f m gm 时, 设 c = ??gm (0＜? ＜ f m ＝, 令 y = y0, 注意到 g ( y0) = gm, 由方程(4-5) 可得代数方程 f (x) = ??, 从图 4-3 知道这个方程有两个根 x1 和 x2, 且 x1＜x0＜x2, 于是在图 4-5 上得到这条轨线上的两个点 S1 (x1, y0), S2 (x1, y0). 对于区间 [x1, x2 ] 之外的任意一点 x, 从图 4-3 知 f (x)＜?, 再由 f (x) g ( y) = ? gm 得 g ( y) ＞gm, 与 gm 的定义矛盾. 对于区间(x1, x2)内的任意一点 x, 因为 f (x)＞?,图 4 -5 f (x) g ( y) = ? gm, 所以 g ( y)＜gm . 从图 4-4 知道方程 g ( y) = ???( ? ＜gm )有两个根 y1 和 y2, 且 y1＜ y0＜y2, 于是在图 4-5 上得到同一条轨线上的另外两个点 S3 和 S4. 因为 x 是任意的, 所 以可知这条轨线是封闭曲线. 通过上面的分析可知, 当 c 从 f m gm 变小而趋于零时, 轨线是一族以平衡点 B (x0, y0)为中 ? ? 心, 越来越扩展的封闭曲线(图 4-6). 轨线的方向可由图上四个区域中 x 和 y 的符号决定.图 4 -6 方程(1)的轨线族图 4 -7封闭轨线对应着方程(4-4)的周期解(图 4-7) x = x (t), y = y (t), 设周期为 T, 利用方程(4-4) 可以求出 x (t ) 和 y (t )在一周期内的平均值 x , y .x ? 1 T? x (t )d t ? T ?0T1T(0d??1 dy y dt) d t = d / ??+ [ ln y (T ) ??ln y (0)]/T = d / ??,类似地有 y = r /??. 即 x (t ) 和 y (t )在一周期内的平均值 x , y 为它们在平衡点 B (x0, y0)的值. 对于周期性变化的 x (t ) 和 y (t ), 我们用它们的平均值 x , y 表示其大小.综上所述, 食肉 鱼的数量取决于两个参数 r 和??. 当食用鱼的自然增长率 r 下降时, 食肉鱼的数量将减少, 这 就是说, 在弱肉强食情况下, 如果人们想使强者减少, 应当去降低弱者的繁殖率, 而当食肉 鱼掠取鱼的能力? 提高时, 对食用鱼没有影响, 只会使食肉鱼减少.另一方面, 食肉鱼死亡率 d 上升将导致食用鱼的增多, 而食用鱼对食肉鱼供养能力 ??的提高会使食用鱼减少. 下面用这个模型解释开始提出的为什么战争时捕捞量下降食用鱼反而减少的问题. 上面我们讨论的是在自然环境中食用鱼和鲨鱼数量的变化规律, 在人工捕鱼的情况下, 设表示捕捞能力的系数为 e, 这相当于食用鱼的自然增长率由 r 降为 r ??e, 鲨鱼的死亡率由 d 升为 d + e.在这种改动下上述模型仍然成立, 于是食用鱼和鲨鱼的(平均)数量 x , y 变为： x = (d + e)/?, y = (r ??e)/??. 因此捕捞能力 e 下降, 鲨鱼的数量( y )增加, 而食用鱼的数量( x )减少. 上述分析用于解释使用杀虫剂的影响时也是出乎人们意料的.害虫靠吃农作物生存, 为 被捕食者；益虫是害虫的天敌, 为捕食者.如果某种杀虫剂不仅杀死害虫, 也能杀死益虫, 那 么增强杀虫能力相当于提高上述分析中的捕捞系数 e, (从长期效果上看)会使害虫( x )增多, 益虫( y )减少, 事与原违！这个事实说明, 捕食--被捕食模型需要改进, 改进后的模型如下：? d x1 ? x1 x2 ? ?, ? r1 x 1 ? 1 ? ??1 ? ? ? N1 N2 ? ? dt ? ? ? x1 x2 ? dx 2 ? r2 x 2 ? ? 1 ? ? 2 ? ? ? dt N1 N2 ? ?? ?. ? ?各符号的含义及其稳定性分析留作习题 8.第二节 差分方程模型 在上一节中,我们利用微分方程方法研究了一些连续变化的变量,但现实世 界中常常会遇到一些问题,其变量只在一定的离散集合中取值,例如,在对经济进 行动态分析时,一般总根据一计划周期期末的指标值来判断某经济计划执行得如 何.有些变量虽然可以认为是连续变化或近似连续变化得,但在研究时,并不时时 刻刻去统计它,而只在某些特定时刻才去统计它.这实际上就是差分的基本思想, 利用这种方法,我们可以将微分方程的变量离散化,得到相应的差分方程模型.§4.2.1 差分方程模型及求解这一节中,我们先介绍差分及差分方程的一些基本概念,然后就几类特殊的 差分方程的解法进行讨论.§4.2.1.1 基本概念以 n 表示时间,规定 n 只取非负整数,记 x n 为变量 x 在 n 时刻的取值,则称? x n ? x n? 1 ? x n 为 x n的一阶差分,称 ? 2 x n ? ? ( ? x n ) ? ? x n ?1 ? ? x n ? x n ? 2 ? 2 x n ?1 ? x n 为 阶差分 ? k x n .x n 的二阶差分.类似地,可以定义 x n 的 k由 n , x n 以及 x n 的差分给出的方程称为 x n 的差分方程,其中含有 x n 的最高阶 差分的阶数称为差分方程的阶.差分方程也可以写成不显含差分的形式.例如,二 阶差分方程 ? 2 x n ? ? x n ? x n ? 0 ,可以改写成 x n ? 2 ? x n ? 1 ? x n ? 0 .因此,k 阶差分方程 也可以写成F ( n , x n , x n ?1 , ? , x n ? k ) ? 0(4.2.1)满足差分方程(4.2.1)的序列 x n ? x ( n ) ,称为差分方程(4.2.1)的解.包含 k 个任 意常数的解称为(4.2.1)的通解 , x 0 , x1 , ? , x k ?1 为已知时,称为(4.2.1)的初始条件 , 通 解中的任意常数都由初始条件确定后的解称为(4.2.1)的特解. 形如x n ? k ? a1 ( n ) x n ? k ? 1 ? ? ? a k ( n ) x n ? f ( n )(4.2.2)的差分方程,称为 k 阶线性差分方程.其中 a i ( n ) 为已知系数,且 a k ( n ) ? 0 . 若差分方程(4.2.2)中 f ( n ) ? 0 ,则称差分方程(4.2.2)为 k 阶齐次线性差分方程, 否则称为 k 阶非齐次线性差分方程.§4.2.1.2 差分方程的求解 一、常系数线性齐次差分方程的解析解法 形如x n ? k ? a1 x n ? k ?1 ? ? ? a k x n ? 0(4.2.3)的差分方程,称为 k 阶常系数线性齐次差分方程,其中 a i 为常数,且 a k ? 0 .代数方程? ? a1 ?k k ?1? ? ? ak ? 0(4.2.4)称为差分方程(4.2.3)的特征方程,其根称为特征根. 根据特征根的情况,我们可以得到方程(4.2.3)的通解. 1) 特征值均为单根. 若差分方程(4.2.3)的特征方程(4.2.4),有 k 个相异的特征根 ?1 , ? 2 , ? , ? k ,则x n ? C 1 ?1 ? C 2 ? 2 ? ? ? C k ? kn n n是方程(4.2.3)的一个通解,其中 C i 为任意常数. 如果有一组初始条件x 0 ? u 0 , x1 ? u 1 , ? , x k ? 1 ? u k ? 1则可以确定出一个满足初始条件的特解. 2) 特征值中有重根 若差分方程(4.2.3)的特征方程(4.2.4),有 l 个相异特征根 ?1 , ? 2 , ? , ? l ,重数依次 为 m 1 , m 2 , ? m l .显然 ? m i ? k .则差分方程(4.2.3)的通解为i ?1 lm1 m2 mixn ??C1 j nj ?1?1 ?nj ?1?C2 jnj ?1?2 ? ? ?nj ?1?Cj ?1ijnj ?1?in3) 特征值中有复根 若 差 分 方 程 (4.2.3) 的 特 征 方 程 (4.2.4) 的 特 征 根 中 出 现 一 对 共 轭 复 根?1, 2 ? ? (co s ? ? i sin ? ) ,和相异的 k ? 2 个根 ? 3 , ? 4 , ? , ? k ,则差分方程(4.2.3)的通解为x n ? C 1 ? co s n? ? C 2 ? sin n? ? C 3 ? 3 ? ? ? C k ? kn n n n.注 1:实际上,如果我们把复根的形式用指数表示,即 ?1 ? e ? ? i? , ? 2 ? e ? ? i? .则可以直 接运用情况 1)全为单根的相应结论得到同样的通解,即x n ? C 1 ?1 ? C 2 ? 2 ? ? ? C k ? k ? C 1 en n n n n n ? ? in?? C 2enn ? ? in?? C 3?3 ? ? ? C k ? kn nn? C '1 ? sin n? ? C ' 2 ? co s n? ? C ' 3 ? 3 ? ? ? C ' k ? k其中 C 'i 在复数集中取值. 同理可得,如果出现特征值为复根且为重根的情形,我们可以用同样的方法,将 复根化成指数表达的形式,然后运用情况 2)的结论得到通解. 例 4-11 求 Fibonacci 数,通项为 Fn ? 2 ? Fn ? 1 ? Fn 的通解,并确定满足初始条件F1 ? F 2 ? 1 的特解.1? 2 5 1? 2 5解:该差分方程的特征方程为 ? 2 ? ? ? 1 ? 0 .特征根为 ?1 ?, ?2 ?,是互异的,根据情形 1)的结论,我们可以得到,该差分方程的通解为Fn ? C 1 ( 1? 2 5 ) ? C2(n1? 25)n将初始条件 F1 ? F2 ? 1 代入通解,得C1 ( 1? 2 5 ) ? C2( 1? 2 5 ) ? 1, C 1 ( 1? 2 5 ) ? C2(21? 25) ?12所以 C 1 ?1 5,C2 ? ?1 5,则满足初始条件得特解为1 ? 1? 5 n 1? 5 n? ) ?( ) ? ?( 2 2 5 ? ?Fn ?例 4-12 求差分方程 x n ? 4 ? x n ? 3 ? 3 x n ? 2 ? 5 x n ? 1 ? 2 x n ? 0 的通解,并求满足初始条件x 0 ? 1, x1 ? 0, x 2 ? 1, x 3 ? 2的特解.解:该方程对应的特征方程为? ? ? ? 3? ? 5 ? ? 2 ? 04 3 2特征根为-1,-1,-1,2,出现重根的情形,利用 2)的结论,得方程的通解为x n ? ( C 1 ? C 2 n ? C 3 n )( ? 1) ? C 4 22 n n将初始条件代入通解,得到 C 1 ? 的特解为xn ? 1 5242 52,C2 ??29 52, C3 ?7 52,C4 ?10 52,所以满足初始条件( 4 2 ? 2 9 n ? 7 n )( ? 1) ?2 n10 522n例 4-13 求差分方程 x n ? 2 ? x n ? 0 的通解,并求满足初始条件 x 0 ? 1, x1 ? 2 时的特解. 解:该方程是齐次线性方程,其对应的特征方程为 ? ?1 ? 02对应的特征根为一对共轭复根 ?1 ? i ? co sn? 2?2? i sin?2, ? 2 ? ? i ? co s?2? i sin?2.所以根据情形 1)和 3)的结论,我们都可以得到,该差分方程的通解为x n ? C 1 co s ? C 2 sin n? 2当 x 0 ? 1, x1 ? 2 时,代入 x n 的表达式,求得 C 1 ? 1, C 2 ? 2 ,所以满足初始条件的特解为x n ? co s n? 2 ? 2 sin n? 2二、系数线性非齐次差分方程的解析解法 形如x n ? k ? a1 x n ? k ? 1 ? ? ? a k x n ? f ( n )(4.2.5)的差分方程,称为 k 阶常系数线性齐次差分方程,其中 a i 为常数,且 a k ? 0 , f ( n ) ? 0 . 差分方程(4.2.5)对应的齐次差分方程为(4.2.3). 与微分方程类似,非齐次差分方程的解也满足叠加原理.即非齐次差分方程 (4.2.5)的通解等于对应的齐次差分方程(4.2.3)的通解加上非齐次方程的特解. 非齐次差分方程(4.2.5)特解的求解方法,与非齐次方程微分方程类似,可以通 过 常 数 变 异 法 , 对 一 些 特 殊 的 f (n) 函 数 , 还 可 以 用 比 较 系 数 法 . 例 如 , 当f (n ) ? b pl (n )n时,其中 p l ( n ) 为 n 的 l 次多项式.1)若 b 不是特征根,则非齐次方程(4.2.5)有形如 b n q l ( n ) 的特解, 其中 q l ( n ) 为 n 的 l 次多项式; 2)若 b 是 r 重特征根,则非齐次方程(4.2.5)有形如 b n n r g l ( n ) 的特解, 其中 g l ( n ) 为 n 的 l 次多项式; 进而,通过待定系数法,将带有待定未知常数的特解代回原方程,可以求出q l ( n ) , g l ( n ) ,从而求得非齐次方程的特解.另外,由于特解的随意性,只要能找到一个满足方程的特殊解就可以解决问题, 所以有时对于一些比较简单的差分方程,我们可以直接通过观察得到特解. 例 4-14 求解差分方程 x n ? 1 ? 2 x n ? 1 ,满足初始条件 x 0 ? 1 的解. 解:该方程是一个简单的非齐次问题,其对应的齐次方程为x n ?1 ? 2 x n ? 0对应的特征方程为 ? ? 2 ? 0 ,特征根为 ? ? 2 ,所以齐次方程的通解为xn ? C1 2n通过观察我们发现,如果我们取 x n ? a 始终为常数,代入原差分方程可以找到 a ? 2 a ? 1 ,即 x n ? a ? ? 1 ,为满足非齐次方程的一个特解.根据叠加原理,我们得到非齐次方程的特解为xn ? C1 2 ? 1n代入初始条件 x 0 ? 1 ,得到 C 1 ? 2 .所以满足初始条件的原方程的解为xn ? 2n ?1?1例 4-15 求差分方程 x n ? 2 ? 4 x n ? 1 ? 4 x n ? 2 n 的通解. 解: 该方程为非齐次差分方程,其对应的齐次方程的特征方程为? ? 4? ? 4 ? 02求得特征根为二重根 ?1 ? ? 2 ? 2 ,所以齐次方程的通解为xn ? C1 2 ? C 2 n 2n n根据题述非齐次差分方程的右端 f ( n ) ? 2 n ,又因为 2 是特征方程的二重根,可得非 齐次方程有形如xn ? C 3n 22 n的特解.将 x n 代入原方程,解得 C 3 ?n1 8n.所以非齐次方程的通解为n n 2 n?3xn ? C1 2 ? C 2 n 2 ? xn ? C1 2 ? C 2 n 2 ? n 2三、计算机求解 若 x0, x1, ? , xk?1 已知,则形如 xn+k = g (n; xn, xn+1, ? , xn+k?1 )的差分方程的解可 以在计算机上直接用循环语句实现.这种方法虽然不能求得差分方程的解的具体 表达式,但是可以知道每个时刻的数值解,而且适用于一切有显示表达的差分方 程. 例 4-16 计算差分方程xn ? 4 ? x2 n?3? 3 x n ? 2 ? 5 x n ?1 ? 2 x n ? 2n满足初值条件 x 0 ? 1, x1 ? 0, x 2 ? 1, x 3 ? 2 的解的前 10 项. 解:该差分方程含有非线性项,根据我们之前讨论的求解方法,不能处理这类问题, 我们只能通过通项求解差分方程的有限项的数值解. 在 matlab 中输入命令 &&x(1)=1;x(2)=0;x(3)=1;x(4)=2; %赋初值,注意,matlab 中的数组中元素的下标从 1 开始 k=4; n=10; %差分方程阶数 %求解项数 for i=k+1:n x(i)=-(x(i-1))^2+3*x(i-2)+5*x(i-3)+2*x(i-4)+2^(i-4); %通项公式 end x 运行结果: x= 1 0 1 2 3 6 -7 4 31 -908 即为所求差分方程解的前 10 项. 对于我们之前求解的常系数差分方程,也可以用此方法求解数值解.例如,例 4-14 中的差分方程,我们将上述命令稍做改动,便可以得到它的解的前 20 项. 编写命令 &&x(1)=1; k=1; n=20; for i=k+1:n x(i)=2*x(i-1)+1; end x 运行结果: x = 1 3 7 15 31 63 127 255 511 95 8191 %输出结果 %通项公式 %差分方程阶数 %求解项数 %输出结果
287 1048575 运用我们在例 4-14 中得到的解的通项公式 x n ? 2 n ? 1 ? 1 ,可以计算出第 20 项x1 9 ? 220? 1 ? 1 0 4 8 5 7 5 ,与我们的数值解是一致的.§4.2.2 差分方程稳定性理论简介与微分方程模型类似,差分方程平衡点的稳定性研究也是有重大意义的.在 这一节中,我们先给出差分方程平衡点及其稳定性的概念,再对几类特殊的差分 方程的平衡点的稳定性进行具体研究. 若有常数 a 是差分方程(4.2.1)的解,即 F ( n , a , a , ? , a ) ? 0 ,则称 a 是差分方程 (4.2.1)的平衡点.若对差分方程(4.2.1)由任意初始条件确定的解 xn 都有xn ? a (n ? ? )则称这个平衡点 a 是稳定的. 一、一阶线性差分方程的平衡点和稳定性 一阶线性常系数差分方程x n ?1 ? a x n ? b(4.2.6) .的平衡点可以由 x ? ax ? b 解得,为 x * ?b 1? a注意到,通过变量代换可以把方程(4.2.6)的平衡点的稳定性问题转换为x n ?1 ? a x n ? 0(4.2.7)的平衡点 x * ? 0 的稳定性问题.根据上一节介绍的齐次线性差分方程的解法,可知, 方程(4.2.7)的解可以表示为xn ? ( ? a ) x0n根据稳定性的定义,可知,当且仅当 a ? 1 时,方程(4.2.7)的平衡点才是稳定的, 从而方程(4.2.6)的平衡点也是稳定的. 对于 m 维向量 x ( n ) 和 m ? m 常数矩阵 A 构成的方程组x ( n ? 1) ? A x ( n ) ? 0(4.2.8)其平衡点稳定的条件是 A 的特征值 ? i , i ? 1, 2, ? ,均有?i ? 1即均在复平面上的单位圆内. 二、二阶线性差分方程的平衡点和稳定性 二阶齐次线性差分方程x n ? 2 ? a 1x n ? 1 a x2n ? 0 ?(4.2.9)的平衡点为 x * ? 0 .差分方程(4.2.9)的特征方程为 ? 2 ? a1 ? ? a 2 ? 0 ,记特征根为? 1 , ? 2 ,则方程(4.2.9)的通解可以表示为x n ? C 1? 1 ? C ? 2n n 2(4.2.10)其中常数 C 1 , C 2 可以由初始条件 x 0 , x1 确定.由(4.2.10)可得,当且仅当?i ? 1( i ? 1, 2 )时方程(4.2.9)的平衡点才是稳定的. 与一阶线性方程类似,通过变量代换可得,非齐次线性方程x n ? 2 ? a 1x n ? 1 a x2n ? b ?(4.2.11) 的平衡点的稳定性和方程(4.2.9)相同. 二阶方程的上述结果可以推广到 k 阶线性方程,即 k 阶方程x n ? k ? a1 x n ? k ?1 ? ? ? a k x n ? b的平衡点稳定的充要条件是其对应的特征方程的特征根 ? i ,均有?i ? 1i ?( 1, 2 , ? n. , )(4.2.12)注 2：高阶方程与一阶方程组通过一定的变量代换是可以相互转化的,我们注意 到,(4.2.12)这一条件与一阶方程组(4.2.8)的平衡点稳定性条件是一致的. 三、一阶非线性差分方程的平衡点和稳定性 一阶非线性差分方程x n ?1 ? f ( x n )(4.2.13)的平衡点 x * 可以由代数方程 x ? f ( x ) 解出.为分析 x * 的稳定性,我们采用与微分 方程类似的线性近似的方法,将方程(4.2.13)的右端在 x * 处作 Taylor 展开,只取到 一次项,得到方程(4.2.13)的线性近似方程x n ? 1 ? f ( x ) ? f '( x )( x n ? x )* * *(4.2.14)显然, x * 也是(4.2.14)的平衡点.线性方程(4.2.14)的平衡点的讨论与(4.2.6)相 同,而当 f '( x * ) ? 1 时,方程(4.2.14)与(4.2.13)的平衡点的稳定性相同.于是我们得到 ①当 f '( x * ) ? 1 时,方程(4.2.13)的平衡点 x * 是稳定的; ②当 f '( x * ) ? 1 时,方程(4.2.13)的平衡点 x * 是不稳定的.§4.2.3 市场经济的蛛网模型在自由贸易的集市上你注意过这样的现象吗？一个时期由于猪肉的上市量远大于需求, 销售不畅导致价格下降, 农民觉得养猪赔钱, 于是转而经营其它农副业.过一段时间后猪肉 上市量大减, 供不应求导致价格上涨, 原来的饲养户看到有利可图, 又重操旧业.这样下一个 时期会重现供大于求、价格下降的局面, 在没有外界干扰的情况下, 这种现象将如此循环下 去. 在完全自由竞争市场经济中上述现象通常是不可避免的, 因为商品的价格是由消费者 的需求关系决定的, 商品数量越多价格越低.而下一时段商品的数量由生产者的供应关系决 定, 商品价格越低生产的数量就越小.这样的需求和供应关系决定了市场经济中商品的价格 和数量必然是振荡的.这种振荡越小越好, 如果振荡太大就会影响人民群众的正常生活. 现在我们用差分方程建模, 讨论市场经济趋于稳定的条件, 再用图形方法建立“蛛网模 型”对上述现象进行分析, 对结果进行解释, 然后作适当推广. 记第 n 时段商品数量为 xn, 价格为 yn, n = 1, 2, ?. 这里我们把时间离散化为时段, 1 个 时段相当于商品的 1 个生产周期, 如蔬菜、水果可以是一年, 肉类则是 1 个饲养周期. 在 n 时段商品的价格 yn 取决于数量 xn, 设 yn = f (xn ). 它反映消费者对这种商品的需求关 系, 称需求函数.因为商品的数量越多价格越低, 所以在图 4-8 中用一条下降曲线 f 表示它, f 为需求曲线. 在 n + 1 时段商品的数量 xn+1 由上一时段价格 yn 决定, 用 xn+1 = g ( yn )表示, 它反映生产者的供应关系, 称供应函 数.因为价格越高生产量才越大, 所以在图 4-8 中供应曲线 g 是一条上升曲线. 设图 4-8 中两条曲线 f 和 g 相交于 P0(x0, y0)点, 在 P0 点附近取函数 f 和 g 的线性近似, 即 需求函数( f )：yn ? y0 = ????(xn ??x0), ??＞0 4-13 图 4 -8 供应函数(g)：xn+1 ??x0 = ???( yn ?? y0 ), ??＞0. 4-14 由式(4-13), (4-14)消去 yn 得一阶线性差分方程 xn+1 = ??????xn + (1 + ????) x0, n = 1,2,?. 4-5 因此 x0 是其平衡点, 即 P0 是平衡点, 对 n 递推不难得到 xn+1 = (???????) n x1 +[1? (??????) n ]x0, n = 1,2,?. 由此可得, 平衡点 P0 稳定的条件是：??? ＜1; 不稳定的条件是：??? ＞1. 下面用图形解释此模型. 若对某一个 k 有 xk = x0, 则由式(4-5)递推可得, 当 n≥k 时 xn = x0, 从而 yn = y0, 即商品的 数量和价格将永远保持在 P0(x0, y0)点.但是实际生活中的种种干扰使得 xn, yn 不可能停止在 P0 点上.不妨设 x1 偏离 x0 (如图 4-9, 图 4-10), 我们分析随着 n 的增加 xn, yn 的变化.图 4 -9 P0 点是稳定的图 4 -10 P0 点是不稳定的数量 x1 给定后, 价格 y1 由曲线 f 上的 P1 点决定, 下一时段的数量 x2 由曲线 g 上的 P2 点 决定, 这样得到一系列的点 P1(x1, y1), P2(x2, y1), P3(x2, y2), P4(x3, y2),?,在图 4-9 上这些点将接 箭头所示方向趋向 P0(x0, y0), 表明 P0 是稳定平衡点, 意味着市场经济（商品的数量和价格） 将趋向稳定.但是如果需求函数和供应函数由图 4-10 的曲线所示, 则类似的分析发现, 市场 经济将按照 P1, P2, P3, P4,?, 的规律变化而远离 P0, 即 P0 是不稳定平衡点, 市场经济趋向不 稳定. 由此可见, 需求曲线越平, 供应曲线越陡, 越有利于经济稳定. 图 4-9 和图 4-10 中折线 P1P2P3P4?形似蛛网, 于是这种用需求曲线和供应曲线分析市场 经济稳定性的图示法在经济学中称蛛网模型.实际上, 需求曲线 f 和供应曲线 g 的具体形式通 常是根据各个时段商品的数量和价格的一系列统计资料得到的.一般地说, f 取决于消费者对 这种商品的需要程度和他们的消费水平, g 则与生产者的生产能力, 经营水平等因素有关. 下面再来解释此模型的实际意义. 首先考察参数?, ??的含义, 需求函数 f 的斜率??(取绝对值)表示商品供应量减少 1 个单位 时价格的上涨幅度；供应函数 g 的斜率?表示价格上涨 1 个单位时(下一时期)商品供应增加 量.所以? 的数值反映消费者对商品需求的敏感程度. 如果这种商品是生活必需品, 消费者 处于持币待购状态, 商品数量稍缺, 人们立即蜂拥抢购, 那么??会比较大； 反之, 若这种商品 非必需品, 消费者购物心理稳定, 或者消费水平低下, 则??较小. ? 的数值反映生产经营者 的对商品价格的敏感程度, 如果他们目光短浅, 热衷于追逐一时的高利润, 价格稍有上涨立 即大量增加生产, 那么??会比较大；反之, 若他们素质较高, 有长远的计划, 则???较小. 根据?, ??的意义很容易对市场经济稳定与否的条件作出解释.当供应函数 g, 即???固定 时, ? 越小, 需求曲线越平, 表明消费者对商品需求的敏感程度越小, 越利于经济稳定.当需 求函数 f, 即? 固定时, ???越小, 供应曲线越陡, 表明生产者对价格的敏感程度越小, 越利于 经济稳定.反之, 当??? 较大, 表明消费者对商品的需求和生产者对商品的价格都很敏感, 则 会导致经济不稳定. 当市场经济趋向不稳定时, 政府有两种干预办法：一种办法是控制价格, 无论商品数量 多少, 命令价格不得改变, 于是??= 0, 不管曲线 g 如何, 总是稳定的；另一种办法是控制市 场上的商品数量, 当上市量小于需求时, 政府从外地收购或调拔, 投入市场, 当上市量多于 需求时, 政府收购过剩部分, 于是??= 0, 不管曲线 f 如何, 也总是稳定的. 最后我们将模型稍加推广. 如果生产者的管理水平更高一些, 他们在决定商品生产数量时, 不是仅根据前一时期的 价格, 而是根据前两个时期的价格, 为简单起见不妨设根据二者的平均值(yn + yn?1)/2, 于是 供应函数为 xn+1 = g[(yn + yn?1)/2].在 P0 点附近取线性近似时 供应函数(g)：xn+1 ??x0 = ? [(yn + yn?1)/2? y0], ??＞0. 4-16 又设需求函数仍由(1)式表示, 则由(4-13),(4-16)式得到 2xn+2 + ?? xn+1 + ?? xn = (1 + ?? )x0, n =1,2,?. 4-17 (4-17)式是二阶线性差分方程. P0 点稳定的条件可由其特征方程 2??2 +???? + ??? = 0 的根??1, 2 =? ?? ? (?? ) 2 ? 8?? 4确定.当?? ＞8 时, 显然有??2 =? ?? ?(?? ) 2 ? 8?? 4? ???4＜?2,从而|???2 |＞2, P0 是不稳定的. 当???＜8 时, 可以算出|???1, 2 | =??2, 由|???1, 2 |＜1 得到 P0 点稳定的条件为??? ＜2.这与原有模型中 P0 点稳定的条件??? ＜1 相比, 保持经济稳定的参数??,????的范围放大了(?, ???的含义未变). 可以想到, 这是生产经营者的生产管理水平提高, 对市场经济稳定起着有利 影响的必然结果.参考文献[1]胡良剑,孙晓军.Matlab 数学实验.北京:高等教育出版社,2006. [2]姜启源,谢金星,叶俊.数学模型(第三版).北京:高等教育出版社,2003. [3]刘承平.数学建模方法.北京:高等教育出版社,2002. [4]薛定宇,陈阳泉.高等应用数学问题的 Matlab 求解.北京:清华大学出版社,2006 [5]王高雄,周之铭,朱思铭,王寿松.常微分方程(第二版).北京:高等教育出版社,1983.习题 41. 一角度为 60 度圆锥形漏斗里装着 10cm 高的水, 其下端小孔的面积为 0.6cm2,求这此水流完需用多少 时间？ 2. 某容器温度为 60℃, 将其内的温度计移入另一容器中, 十分钟后读数为 70℃, 又过十分钟后读数为 76℃, 问另一容器的温度为多少度？（提示：温度变化率与温度差成正比） 3. 一大学教师小李从 31 岁开始建立自己的养老基金, 他把已有的积蓄 1 万元也一次性地存入, 已知月 利率为 0.01（以复利计）, 每月存入 300 元, 试问当小李 60 岁退休时, 他的退休基金有多少？又若, 他退休 后每月要从银行提取 1000 元, 试问多少年后他的退休基金将用完？你能否根据你了解的实际情况建立一个 较好的养老基金的数学模型及相应的算法和程序(软件). 4. 某人每天由饭食获取 2500 卡热量, 其中 1200 卡用于新陈代谢, 此外每公斤体重需支付 16 卡热量作 为运动消耗, 其余热量则转化为脂肪, 已知以脂肪形成贮存的热量利用率为 100%, 每公斤脂肪含热量 10000 卡, 问此人的体重如何随时间而变化？ 5. 1650 年世界人口为 5 亿, 当时的年增长率为 3%.用指数增长模型计算什么时候世界人口达到 10 亿. (实际上 1850 年前已超过 10 亿), 1970 年世界人口为 36 亿, 年增长率为 2.1%, 用指数增长模型预测什么时 候世界人口会翻一番(这个结果可信吗？)你对同样的模型得到的两个结果有什么看法. 6. 设某城市共有 n + 1 人, 其中一人出于某种目的编造了一个谣言.该城市具有初中以上文化程度的人 占总人数一半, 这些人只有 1/4 相信这一谣言, 而其他人约有 1/3 会相信.又设凡相信此谣言的人在每单位时 间内传播的平均人数正比于时当未听说此谣言的人数, 而不相信此谣言的人不传播谣言.试建立一个反映谣 言传播情况的微分方程模型. 7. 如果两个种群都能独立生存, 共处时又能相互提供食物, 试建立种群相互依存模型并讨论平衡点的 稳定性, 解释稳定的意义.又如果两个种群都不能独立生存, 试建模以讨论共处的可能性. 8. 请你捕食--被捕食模型加以改进, 并讨论平衡点的稳定性, 解释稳定的意义. 9. 如果在捕食--被捕食模型中, 捕食者掠食的对象只是成年的被捕食者, 而未成年的因体积太小免遭 捕获, 在适当的假设下建立这三者之间关系的模型, 求平衡点. 10. 大陆上物种数目是常数, 各物种独立地从大陆向附近一岛迁移, 岛上物种数量的增加与尚未迁移 的物种数目有关, 而随着迁移物种数的增加又导致岛上物种的减少, 在适当假设下建立岛上物种数的模型, 并讨论稳定状况. 11. 对于技术革新的推广, 在下列几种情况下分别建立模型. ① 推广工作通过已经采用新技术的人进行, 推广速度与已采用新技术的人数成正比, 推广是无限的. ② 总人数有限, 因而推广速度随着尚未采用新技术人数的减少而降底. ③ 在②的前提下还要考虑广告等媒介的传播作用. 12. 在市场经济的蛛网模型中, 继续讨论下列问题： ① 因为一个时段上市的商品不能立即售完, 其数量也会影响到下一段的价格, 所以 yn+1 由 xn+1 和 xn 决 定, 再设 xn+1 仍只取决于 yn, 给出稳定平衡的条件. ② 若除了 yn+1 由 xn+1 和 xn 决定之外, xn+1 也由 yn 和 yn?1 确定, 试分析稳定平衡的条件是否还会放宽. 13. 在按年龄分组的的种群增长模型中, 设一群的动物最高年龄为 15 岁, 每 5 岁一组分成 3 个年龄组, 各组的繁殖率为 b1 = 0, b2 = 4, b3 = 3, 存活率为 s1 = 1/2, s2 = 1/4 , 开始时 3 组各有 1000 只.求 15 年后各组分别 有多少只？以及时间充分长以后种群的增长率和按年龄组的分布. 14. 最优捕鱼策略(1996 年全国大学生数学建模竞赛题 A 题)为了保护人类赖以生存的自然环境, 可再 生资源(如渔业, 林业资源)的开发必须适度, 一种合理, 简化的策略是, 在实际可持续收获的前提下, 追求 最大产量或最佳效益. 考虑对某种鱼(鱼)的最优捕捞策略：假设这种鱼分 4 个年龄组, 称 1 龄鱼, …, 4 龄鱼, 各年龄组每条 鱼的平均重量分别为 5.07, 11.55, 17.86, 22.99(g), 各年龄组自然死亡率均为 0.8(1/年), 这种鱼为季节性集中 产卵繁殖, 平均每条 4 龄鱼的产卵量为 1.109×105(个), 3 龄鱼的卵量为这个数的一半, 2 龄鱼和 1 龄鱼不产 卵, 产卵和孵卵期为每年的最后 4 个月, 卵孵化并成活为 1 龄鱼, 成活率(1 龄鱼条数与产卵总量 n 之比)为 1.22××1011+n). 渔业管理部门, 每年只允许在产卵孵化期前的 8 个月内进行捕捞作业, 如果每年投入的捕捞能力(如渔 船数, 下网次数等)固定不变, 这时单位时间捕捞量将与各年龄组鱼数成正比, 比例系数不妨称为捕捞强度 系数, 通常使用 13mm 网眼的拉网, 这种网只能捕捞 3 龄鱼和 4 龄鱼, 其中两个捕捞强度系数之比为 0.42:1. 渔业上称为这种方式为固定努力量捕捞. ① 建立数学模型分析如何实现可持续捕捞(即每年开始捕捞渔场中各年龄组鱼群条不变)并且在此前 提下得到最高的年收获量(捕捞总重量). ② 某渔业公司承包这种鱼的捕捞业务 5 年, 合同要求 5 年后鱼群的生产能力不受到太大破坏, 已知承 包时各年龄组鱼群的数量分别为：122, 29.7, 10.1, 3.29(×109 条), 如果要用固定努力量的捕捞方式, 该公司 应采取怎样的策略才能使总收获量最高.
更多相关文档}
第四章 微分方程与差分方程方法
第四章 微分方程与差分方程方法第一节 微分方程模型我们在数学分析中所研究的函数,是反映客观现实世界运动过程中量与量之 间的一种关系,但我们在构造数学模型时,遇到的大量实际问题往往不能直接写 出量与量之间的关系,却能比较容易地建立这些变量和它们的导数(或微分)间的 关系式,这种联系着自变量、 未知函数及其导数(或微分)的关系式称为微分方程.§4.1.1 微分方程简介这一节,我们将介绍关于微分方程的一些基本概念. 一、微分方程的阶数 首先我们具体的来看一个微分方程的例子. 例 4-1 物体冷却过程的数学模型 将某物体放置于空气中,在时刻 t ? 0 ,测量得它的温度为 u 0 ? 150 0 C ,10 分 钟后测量得温度为 u 1 ? 100 0 C .我们要求决定此物体的温度 u 和时间 t 的关系,并 计算 20 分钟后物体的温度.这里我们假定空气的温度保持为 u ? ? 24 0 C . 解:根据物理学中的牛顿冷却定律可知,热量总是从温度高的物体向温度低的物 体传导;一个物体的温度变化速度与这一物体的温度与其所在介质温度的差值成 正比.设物体在时刻 t 的温度为 u ? u (t ) ,则温度的变化速度可以用 们得到描述物体温度变化的微分方程du dt ? ? k (u ? u ? ) du dt来表示.我(4.1.1)du dt其中 k ? 0 是比例常数. 方程(4.1.1)中含有未知函数 u 及它的一阶导数 ,这样的方程,我们称为一阶微分方程.微分方程中出现的未知函数最高阶导数的阶数称为微分方程的阶 数. 方程d y dt3 3?bdy dt? cy ? f ( t )(4.1.2)中未知函数最高阶导数的阶数是三阶,则方程(4.1.2)称为三阶微分方程. 二、常微分方程与偏微分方程 如果在微分方程中,自变量的个数只有一个,我们称这种微分方程为常微分 方程;自变量的个数为两个或两个以上的微分方程称为偏微分方程. 方程 ? T2?x2?? T2?y2?? T2?z2? 0(4.1.3)就是偏微分方程的例子,其中 T 是未知函数, x 、 y 、 z 都是自变量.而方程 (4.1.1)(4.1.2)都是常微分方程的例子. 三、线性与非线性微分方程 如果 n 阶常微分方程F ( x, y, dy dx ,? , d y dxn n) ? 0(4.1.4)的左端为关于未知函数 y 及其各阶导数的线性组合,则称该方程为线性微分方程, 否则称为非线性方程.一般的 n 阶线性微分方程具有形式d y dxn n? a1 ( x )dn ?1ydxn ?1? ? ? a n ?1 ( x )dy dx? an (x) y ? f (x)(4.1.5)其中 a i ( x ), f ( x ) ( i ? 1 ? n ) 是关于 x 的已知函数. 当 f ( x ) ? 0 时,称(4.1.5)为 n 阶齐次线性微分方程;当 f ( x ) ? 0 时,称(4.1.5)为n阶非齐次线性微分方程. 例如,方程(4.1.2)是三阶线性微分方程.而方程dy ? dy ? ? y ?0 ? ? ?t dx ? dx ?2(4.1.6)是一阶齐次非线性方程. 四、微分方程初边值问题 我们把含有 n 个独立的任意常数 c1 , c 2 , ? c n 的解y ? ? ( x , c1 , c 2 , ? c n )称为 n 阶微分方程(4.1.4)的通解.为了确定微分方程一个特定的解,我们通常给出 这个解所必须满足的条件,这就是定解条件.常见的定解条件是初始条件. n 阶微 分方程(4.1.4)的初始条件是指如下的 n 个条件:y ( x0 ) ? y0 , d y ( x0 ) dx ? y0 ,?(1)dn ?1y ( x0 )n ?1dx? y0( n ? 1)求解微分方程满足初始条件的解,称为初值问题.求解初值问题的过程,就是 通过初始条件确定通解中的常数,从而求得满足初始条件的特解的过程. 若方程所给出的定解条件,既有自变量初始时刻的值,也有自变量取终值时 的值,则称该问题为边值问题. §4.1.2 微分方程的求解及 Matlab 实现线性微分方程和低阶特殊微分方程往往可以通过解析解的方法求解,但一般 的非线性微分方程是没有解析解的,即使可以求得解析解,参数也很难确定,需要 用数值解的方式求解.具体的求解方法很多,本节我们主要介绍如何利用 matlab 来求解微分方程的解析解和数值解.§4.1.2.1 微分方程的解析解一、基本理论 有些简单的常微分方程可用一些技巧,如分离变量法、积分因子法、常数变 异法、降阶法等,化为可直接积分的方程从而求得显示解.例如,一阶常系数线性 常微分方程dy dt ? ay ? b (a ? 0)可化为dy ay ? b ? dt两边积分可得 y ( t ) 的通解为 y ? C e a t ? 一般的常系数线性微分方程d y dxn nb a.? a1dn ?1ydxn ?1? ? ? a n ?1dy dx? an y ? f ( x)(4.1.7)的解满足叠加原理,即方程(4.1.7)的通解是对应齐次方程d y dxn n? a1dn ?1ydxn ?1? ? ? a n ?1dy dx? an y ? 0(4.1.8)的通解与非齐次方程(4.1.7)的一个特解的和. 一阶常系数线性常微分方程总可用这一思路求得显示解. 高阶线性常系数微分方程可用特征根法求对应齐次微分方程(4.1.8)的通 解.(4.1.8)对应的代数特征方程s ? a1 sn n ?1? a2 sn?2? ? ? a n ?1 s ? a n ? 0若方程的特征根 s i 均可以求出,且两两相异,则方程(4.1.7)的解可以表示为y ( x ) ? C 1es1 x? C 2es2 x? ? ? C nesn x? ? (x)其中, C i 为待定系数, ? ( x ) 是满足方程(4.1.8)的一个特解,这个特解可以通过常数 变异法,比较系数法得到. s i 有重根的情况也可以写出相应的解析解形式,这里我 们不做过多介绍,因为我们下面要介绍一种更简单的计算机求解法. 二、用 matlab 求解线性常系数微分方程解析解 dsolve 函数调用格式y= d so lve('f 1 ','f 2 ', ? , ' f m ') y= d so lve('f 1 ','f 2 ', ? , ' f m ','x ') 指明自变量其中, f i 既可以描述微分方程,也可以描述初边值条件.在描述微分方程时,用字 母 D 代表微分,后面加数字表示微分的阶数,例如 D4y 表示 y ( 4 ) ( t ) ,也可以用 D2y(0)=3 这类记号来表示初值条件 y ''(0) ? 3 .一般而言任何 D 后面所跟的字母 为因变量,自变量默认为 t ,也可以在函数调用时指明. 例 4-2 求方程 解:输入命令 &&syms t,x; x=dsolve(‘D2x-2*Dx-3*x=exp(-t)’) latex(x) 运行结果 x =exp(-t)*C2+exp(3*t)*C1-1/4*t*exp(-t) ans ={e^{3t}}{\it C2}+{e^{-t}}{\it C1}-1/4\,t{e^{-t}} 原方程的通解为x ? C 1e3td x dt22?2dx dt? 3x ? e?t的通解.%定义符号变量 %求解方程 %用 LATEX 语句显示结果? C 2e?t?1 4te?t其中, C i 为任意常数.若给出初始或边界条件,则可以通过这些条件建立方程,求 出 C i 的值. 仍考虑上面的微分方程,假设已知 x (0 ) ? 3, 令求满足该方程的特解. && x=dsolve(‘D2x-2*Dx-3*x=exp(-t)’,’x(0)=3’,’Dx(0)=2’) 运行结果 x =27/16*exp(-t)+21/16*exp(3*t)-1/4*t*exp(-t)dx dt (0 ) ? 2,则可以通过下面的命 所以，满足这一初值问题的特解是x ? 21 16 e3t?27 16e?t?1 4te?t对于线性方程组我们可以通过定义多维的因变量同样求解. 例 4-3 试求解下面的线性微分方程? dx ? 2x ? 3 y ? 3z ? dt ? ? dy ? 4x ? 5 y ? 3z ? dt ? ? dz ? 4x ? 4y ? 2z ? ? dt解:输入命令 &&[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')； x=collect(simple(x)) %化简 x 的表达式 y=collect(simple(y)) z=collect(simple(z)) 结果为x ? C 2e y ? C 2e z ? C 3e?t ?t? C 3e ? C 3e ? C 1e2t ?2t2t? C 1e2t?2 t部分特殊的非线性微分方程也可以用 dsolve()函数来求解析解.这样的方程 描述方式和前面介绍的线性微分方程是一致的.下面我们通过例子来演示非线性 方程的解析解求解问题,同时还将演示不能求解的例子. 例 4-4 试求出一阶非线性微分方程 解: 输入命令 &&x=dsolve('Dx=x*(1-x^2)') 运行结果 x= 1/(1+exp(-2*t)*C1)^(1/2) -1/(1+exp(-2*t)*C1)^(1/2) 即该方程的解析解为 x ( t ) ? ?1 1 ? C 1e?2 tdx dt? x ( t )(1 ? x ( t )) 的解析解.2.但是如果稍微改变原方程,例如将等号右侧加上 1,则可以用下面语句试解 方程.我们会发现该方程是没有解析解的. &&x=dsolve('Dx=x*(1-x^2)+1') 运行结果 Warning: Explicit solutio implicit solution returned. & In dsolve at 312 x= t+Int(-1/(_a-_a^3+1),_a = .. x)+C1 = 0 从上例我们可以看出,可以求得解析解的方程是很有限的,对于一般的非线 性方程只能用数值解法去求解.下面我们介绍微分方程数值解的求解方法.§4.1.2.2 微分方程的数值解一、基本理论 (一) 数值解的定义 考虑一阶常微分方程组初值问题? y ' ? f (t , y ) ? ? y (t0 ) ? y 0 t0 ? t ? t f(4.1.9)1 其中 y ? ( y 1 , y 2 , ? , y m ) T , f ? ( f 1 , f 2 , ? , f m ) T , y 0 ? ( y 0 , y 02 , ? , y 0m ) T ,这里 T 表示转置.微分方程的数值解,是指不求出解 y ( t ) 的解析表达式,而是寻求 y ( t ) 在一系 列离散结点 t 0 ? t1 ? ? ? t n ? t f 上的近似值 y k ( k ? 0,1,? , n ) .称 hk ? t k ? 1 ? t k 为步长. 通常,取为常数 h ? ( t n ? t 0 ) / n ,此时结点变为等距结点,有t k ? t k ?1 ? h ? t 0 ? kh.(二) 求解数值解的方法 1.欧拉法 基本思想:在结点处用差商近似替代导数y (t k ? h ) ? y (t k ) h ? y '( t k ) ? f ( t k , y ( t k )) ? f ( t k , y k )这样导出欧拉格式的计算公式? y k ?1 ? y k ? h f ( t k , y k ) ? ? y 0 ? y (t0 ) k ? 0,1, ? n ? 1(4.1.10)欧拉法可以求解各种形式的微分方程,但是它只有一阶精度,所以实际应用 效率比较差. 2.改进的欧拉法 基本思想: 使用数值积分离散化.对方程 y ' ? f ( t , y ) ,两边由 t k 到 t k ? 1 积分, 得到y ( t k ?1 ) ? y ( t k ) ?t k ?1 tk?f ( t , y ( t )) d t对右端的定积分用数值积分方法做离散化,可得计算公式,如用矩形公式可 得欧拉公式,若用梯形公式可得改进的欧拉公式,下面我们用梯度公式来做数值 积分y ( t k ?1 ) ? y ( t k ) ??t k ?1 tkf ( t , y ( t )) d t ?t k ?1 ? t k 2[ f ( t k , y ( t k )) ? f ( t k ? 1 , y ( t k ? 1 ))]故有如下迭代计算公式:h ? ? y k ? 1 ? y k ? [ f ( t k , y k ) ? f ( t k ? 1 , y k ? 1 )] 2 ? ? y ? y (t ) 0 ? 0(4.1.11)实际应用时,与欧拉公式结合使用:? y k( 0 )1 ? y k ? h f ( t k , y k ) ? ? ? ( l ? 1) h (l ) ? y k ? 1 ? y k ? [ f ( t k , y k ) ? f ( t k ? 1 , y k ? 1 )] ? 2l ? 0,1, 2, ?(4.1.12)对于已给的精度 ? ,当满足 y k( l??11) ? y k( l?)1 ? ? 时,取 y k ? 1 ? y k( l??11) ,然后继续下一 步 y i ? 2 的计算. 3.Runge-Kutta 方法 基本思想:利用 Taylor 展式进行离散化.假设(4.1.9)中的 f ( t , y ) 充分光滑,将y ( t k ?1 ) 在 t k点作 Taylor 展开:y ( t k ? 1 ) ? y ( t k ) ? h y '( t k ) ? h22!y ''( t k ) ? ? ?hpy( p)p!(t k ) ? ?其中y '( t ) ? f ( t , y ( t )) y ''( t ) ? [ f ( t , y ( t ))] 't ? f t ? f ? f y ? y( p)( t ) ? [ f ( t , y ( t ))] t( p)对照标准形式 y k ? 1 ? y k ? h ? ( t k , h ) .若取 ? ( t , h ) ? y '( t ) ?h 2!y ''( t ) ? ? ?hp ?1y p!( p)(t )并用 y k 代替 y ( t k ) ,则得到一个 p 阶近似公式? y k ?1 ? y k ? h ? ( t k , h ) ? ? y 0 ? y (t0 ) k ? 0,1, ? n ? 1(4.1.13)显然 p=1 时,式(4.1.13)就是(4.1.10),即为欧拉方法.但当 p ? 2 时,如果直接 利用公式(4.1.13)求解数值解问题,就需要计算 f ( t , y ) 的高阶微商.这个计算量是 很大的.因此,直接利用式(4.1.13)构造高阶公式是不实用的. 通常,在求解高阶精度的数值解时,我们并不是直接使用 Taylor 级数,而是利 用它的思想,即计算 f ( t , y ) 在不同结点的函数值,然后作这些函数值的线性组合, 构造近似公式,式中有一些可供选择的参数.将近似公式与 Taylor 展开式相比较, 使前面的若干项密合,从而使近似公式达到一定的精度. 若线性组合选取的函数值个数为 n ,计算达到的精度阶数为 p ,这样的 Runge-Kutta 方法就是我们通常所说的 p 阶 n 级 Runge-Kutta 方法. 后来,德国学者 Felhberg 对传统的 Runge-Kutta 方法进行了改进,使原来的定 步长算法变为自动变换步长的 Runge-Kutta-Felhberg 方法,保证了更高的精度和数 值稳定性,它已经成为求解数值解问题中最常见的方法.具体的算法我们在这里不 做介绍,下面我们主要来讨论如何通过 matlab 软件来运用这一方法求解微分方程 的数值解问题. 二、用 matlab 求解微分方程数值解 (一) 一阶微分方程组初值问题 Matlab 下求解一阶微分方程组初值问题(4.1.9)的数值解的最常见的方法 ode45()函数,该函数实现了前面介绍的变步长四阶五级 Runge-Kutta-Felhberg 算 法,可以采用变步长的算法求解微分方程. ode45()函数的调用格式[t, y]= o d e4 5 (o d efu n ,[t 0 ,t f ],y 0 ) [t, y]= o d e4 5 (o d efu n ,[t 0 ,t f ],y 0 ,o p tio n ,p 1 ,p 2 , ? p m ) 常用格式 完整格式参数说明 [t, y] 中 t 表示自变量,是标量.y 表示函数,可以是标量也可以是向量. y 0 表示 初值.当 y 为 n 维向量时,表明求解的问题是有 n 个未知函数的方程组,则 y 0 也应 为 n 维向量. [ t 0 , t f ] 表示微分方程的求解区间,即自变量的取值范围.如果只给出 一个值tf则表示初始时刻为 t 0? 0 .另外,该函数还允许 t 0 ? t f,即可以认为 t 0 为终值时刻, t f 为初始时刻,则 y 0 表示状态变量的终值,而该函数可以直接求解这样 一个终值问题. odefun 表示待解的微分方程(4.1.9)中 f ( t , y ) 所写成的 m-文件的名称,m-文件 的编写格式为 function ydot=odefun(t,y) ydot=f(t,y); 也可以不把函数单独用 m-文件保存,而是直接用 inline()函数输入,从而获得函数 名,编写格式为 odefun=inline(‘f(t,y)’,’t’,’y’); 其中,t 是自变量,即使所求方程是自治系统,也需要写出 t,否则 matlab 在变量传递 过程中将出现问题.y 为未知函数,ydot 表示未知函数的导数. 如果求解的问题是微分方程组问题,即我们之前讨论过的 y 为向量的情况,则 在书写上述 m-文件或 inline()函数时,待解方程 f ( t , y ) 应以 ydot 的分量形式写成, 或者仍用单个方程的编写格式,多个方程用[ ]括起来,并以;隔开,在后面的例题 中我们将具体的看到. options 表示控制选项.在微分方程求解中有时需要对求解算法及控制条件 进行进一步设置,这可以通过求解过程中的 options 变量来进行修改.初始 options 变量可以通过 odeset()函数来获取,该变量是一个结构体变量,其中有众多成员变 量,表 4-1 列出了常用的一些成员变量.表 4-1 微分方程求解函数的控制参数表参数名 RelTol参数说明为相对误差容许上限,默认值为 10 较高的精度,可以适当减小该值. 为绝对误差容许上限,默认值为 10 分量允许的绝对误差.?6?3,在一些特殊的微分方程求解时,为了保证AbsTol MaxStep,它是一个向量,每个分量表示各未知函数为求解方程允许的最大步长 Mass Jacobian微分代数方程中的质量矩阵,用于描述微分代数方程 为描述 Jacobi 矩阵函数的函数名,如果已知 Jacobi 矩阵,可以加速仿真过程修改 option 变量有两种方式: ①用 odeset()函数设置,其格式为 options=odeset(‘RelTol’,1e-7); 达到的目的是,将相对误差上限设置为较小的 10 ? 7 ②直接修改 options 的成员变量,也可以达到这样的目的,其格式为 options= options.RelTol=1e-7;p1 , p 2 , ? , p m 表示附加参数.即在函数表达式中,除了 t , y还有其他的可变的参数,引入这样的附加参数,使得微分方程的某些参数可以选择不同的值,对不同 值求解时,考虑附加参数可以避免每次修改模型文件.但是值得注意的是,这种情 况下的函数定义格式有特殊要求,格式为 function ydot=odefun(t,y,flag,p1,p2,?,pm)ydot=f(t,y,p1,p2,?,pm); 或 odefun=inline(‘f(t,y,p1,p2, ?,pm)’,’t’,’y’,’flag’,’p1’,’p2’,?,’pm’); 注意,这里的 flag 变量不能省略,用于占位. 例 4-5 解微分方程组? x ' ? ? x3 ? y ? 3 ?y'? x? y ? x (0) ? 1, y (0) ? 0.5 ?0 ? t ? 30解: 该方程是非线性微分方程,不能用我们之前讨论的求解解析解的方法来处理, 只能用数值解法求解.首先,我们将待解方程写成 m-文件,用文件名 eg1fun.m 保 存.我们将 2 个未知变量 x,y,以分量的形式写成一个向量变量 x. %函数 eg1fun.m function f=eg1fun(t,x)f(1)=-x(1)^3-x(2); f(2)=x(1)-x(2)^3; f=f(:); %保证 f 为列向量 这个 m-文件还有另一种等价的书写方式. %函数 eg1fun.m 的另一种书写方式 function f=eg1fun(t,x)f=[-x(1)^3-x(2);x(1)-x(2)^3]; 这样我们就完成了对待解方程的描述,然后输入命令 &&[t,x]=ode45('eg1fun',[0 30],[1;0.5]); subplot(1,2,1); plot(t,x(:,1),t,x(:,2),':'); title('函数曲线图'); xlabel('t'); legend('x(t)','y(t)'); axis square subplot(1,2,2); plot(x(:,1),x(:,2)); title('相平面图'); xlabel('x'); ylabel('y'); axis square 运行结果: %求解方程 %分 2 部分绘图,编辑第一部分 %在第一个图中绘制 x,y 的函数图象 %图象名称 %横轴说明 %说明图例 %显示正方形的坐标系 %编辑第二部分图象 %在第二个图中绘制相平面图 %图象名称 %横轴说明 %纵轴说明 %显示正方形的坐标系图 4-1 例 4-5 函数曲线和相平面图 运行程序,得到函数图象和相平面图.根据相平面上的相轨线,可以看出解从初值 开始的变化趋势,这对于研究微分方程平衡点的稳定性是很有意义的,在后面的章 节中我们将具体讨论. 例 4-6 解微分方程组y1 ' ? y 2 y 3 ? ? y 2 ' ? ? y1 y 3 ? ? y 3 ' ? ? 0 . 51 y 1 y 2 ? ? y 1 ( 0 ) ? 0 , y 2 ( 0 ) ? 1, y 3 ( 0 ) ? 1 ?解:此问题仍为非线性微分方程组,我们求它的数值解,不妨将求解区间取为t ? [0, 2 0 ] .下面我们用 inline()函数来描述方程组,这样就不用单独编写 m-文件来描述方程. 输入命令: &&eg2fun=inline('[y(2)*y(3);-y(1)*y(3);-0.51*y(1)*y(2)]','t','y'); tf=20; y0=[0;1;1]; [t,y]=ode45(eg2fun,[0,tf],y0); plot(t,y(:,1),'-',t,y(:,2),'*',t,y(:,3),'+'); %绘出解的图象 legend('y1','y2','y3'); %说明图例 运行得到函数曲线图图4-2 例4-6函数曲线图例4-7 编写带有附加参数的matlab函数,来描述食饵-捕食者模型? x '1 ( t ) ? x1 ( a 1 ? b1 x 2 ) ? ? x ' 2 ( t ) ? ? x 2 ( a 2 ? b 2 x1 ) ? x (0 ) ? x (0 ) ? 1 2 ? 1 并对 a1 ? 3, b1 ? 2, a 2 ? 2.5, b 2 ? 1 和 a1 ? 2, b1 ? 1, a 2 ? 1, b 2 ? 0 .5 分别求解,画出解的曲 线图和相轨线图. 解:选定附加参数为 a1 , b1 , a 2 , b 2 ,编写如下 m-文件来描述题述方程 %函数 eg3fun.m function xdot=eg3fun(t,x,flag,a1,b1,a2,b2) xdot=[x(1)*(a1-b1*x(2));-x(2)*(a2-b2*x(1))]; 注意,这里的 flag 变量不能省略.这样我们就完成了对方程的描述.我们不妨 取自变量区间为[0,10],对方程求数值解. 输入命令 &&tf=10;x0=[1;1]; a1=3;b1=2;a2=2.5;b2=1; [t,x]=ode45(‘eg3fun’,[0,tf],x0,[],a1,b1,a2,b2); subplot(1,2,1);plot(t,x(:,1),t,x(:,2),':');title('函数曲线图'); xlabel('t');legend('x1','x2'); subplot(1,2,2);plot(x(:,1),x(:,2));title('相平面图'); xlabel('x1');ylabel('x2'); 运行结果图 4-3a1 ? 3, b1 ? 2, a 2 ? 2.5, b 2 ? 1当 a1 ? 2, b1 ? 1, a 2 ? 1, b 2 ? 0 .5 时,只需要修改上述命令中的给 a1,b1,a2,b2 的赋 值语句,就可以得到需要的结果. 图 4-4a1 ? 2, b1 ? 1, a 2 ? 1, b 2 ? 0.5(二) 一阶微分方程组边值问题 对于一阶常微分方程组的边值问题? y ' ? f (t , y ) ? ? g ( y ( a ), y ( b )) ? 0 a ?t ?b(4.1.14)这里 y,f,g 都可以是向量. 用 matlab 求(4.1.14)的数值解的基本格式 sinit=bvpinit(tint,yinit); 由在粗略结点 tinit 的预估解 yinit 生成粗略解网络 sinit. sol=bvp4c(odefun,bcfun,sinit) 其中 odefun 是微分方程组函数,bcfun 为边值条件函数,sinit 是由 bvpinit 得到的粗 略解网络.求得的边值问题的解 sol 是一个结构体变量,sol.x 为求解结点,sol.y 为数 值解. sx=deval(sol,ti) 计算由 bvp4c 得到的解在 ti 的值. 例 4-8 求解边值问题? y '1 ? y 2 ? ? y '2 ? ? y1 ? ? y 1 ( a ) ? 0, y 1 ( b ) ? 2 ? 0a ?t ?b解:输入命令 &&sinit=bvpinit(0:4,[1;0]); odefun=inline('[y(2);-abs(y(1))]','t','y'); bcfun=inline('[ya(1);yb(1)+2]','ya','yb'); sol=bvp4c(odefun,bcfun,sinit); t=linspace(0,4,101); y=deval(sol,t); plot(t,y(1,:),sol.x,sol.y(1,:),'o',sinit.x,sinit.y(1,:),'x'); legend('解曲线','求解点','粗略解'); 运行结果图 4-5 例 4-8 函数曲线图(三) 高阶微分方程 由于 matlab 只能处理一阶微分方程组问题,对于高阶微分方程的初边值问题, 我们必须通过变量代换,使它转换为一阶方程组,才能用 matlab 求解数值解.我们 考虑高阶常微分方程的初值问题y(n)? f ( t , y , y ', ? , y( n ? 1))(4.1.15)且初始值 y (0), y '(0), ? y ( n ?1) (0) 为已知. 令 x1 ? y , x 2 ? y ', ? , x n ? y ( n ?1) ,通过这样的变量代换,高阶微分方程(4.1.15)就 转化为一阶方程组的形式 ? x '1 ? ? x '2 ? ? ?x ' ? n? x2 ? x3 ? ? f ( t , x1 , x 2 , ? , x n )(4.1.16)且初值条件为 x1 (0) ? y (0), x 2 (0) ? y '(0), ? , x n (0) ? y ( n ?1) (0) .这样,我们就可以通过 之前讨论过的一阶微分方程组的数值解法来研究这一问题. 对于高阶微分方程组,只要依次将每个未知变量的各阶导数均单独定义变量, 就可以同样的转化为一阶微分方程组. 例 4-9 用数值解法求解微分方程? y '' ? ( y 2 ? 1) y ' ? y ? 0 ? ? y (0 ) ? ? 0 .2, y '(0 ) ? ? 0 .7解:首先我们需要将这个 2 阶微分方程转化为一阶微分方程组,我们令x1 ? y, x 2 ? y ' ,则原方程转换为? x '1 ? x 2 ? 2 ? x ' 2 ? ? ( x1 ? 1) x 2 ? x1 ? x (0 ) ? ? 0 .2 , x (0 ) ? ? 0 .7 2 ? 1对[0,20]区间求数值解,输入命令 &&tf=20;x0=[-0.2;-0.7]; eg5fun=inline('[x(2);-((x(1))^2-1)*x(2)-x(1)]','t','x'); [t,y]=ode45(eg5fun,[0,tf],x0); plot(t,y); legend('y','y'''); 运行得到函数图象 图 4-6 例 4-9 函数曲线图例 4-10 求解微分方程组(竖直加热板的自然对流问题)2 2 ?d3 f ? df ? d f ?3f ? 2? ? ? ?T ? 0 3 2 d? ? d? ? ? d? ? 2 dT ?d T ? 2 .1 f ?0 ? 2 d? ? d? 2 ? df d f dT (0 ) ? 0, (0 ) ? 0 .6 8, T (0 ) ? 1, (0 ) ? ? 0 .5 ? f (0 ) ? 0, 2 d? d? d? ? ?解:这是一个高阶方程组问题,我们需要通过变量代换转化为一阶微分方程组,从 而寻求数值解.我们令y1 ? f , y 2 ? df d? , y3 ? d f d?2 2, y4 ? T , y5 ?dT d?原方程转化为如下一阶方程组dy ? d y1 ? y2 , 2 ? y3 ? d? d? ? ? dy3 2 ? ? 3 y1 y 3 ? 2 y 2 ? y 4 ? d? ? ? dy dy ? 4 ? y 5 , 5 ? ? 2 .1 y 1 y 5 d? ? d? ? y (0 ) ? 0, y (0 ) ? 0, y (0 ) ? 0 .6 8, y (0 ) ? 1, y (0 ) ? ? 0 .5 2 3 4 5 ? 1输入命令 &&y0=[0;0;0.68;1;-0.5];tf=5; eg6fun=inline('[y(2);y(3);-3*y(1)*y(3)+2*y(2)^2-y(4);y(5);-2.1*y(1)*y(5)]','t','y'); [t,y]=ode45(eg6fun,[0,tf],y0); plot(t,y(:,1),t,y(:,4),':'); legend('f','T'); 运行得到函数图象图 4-7 例 4-10 函数曲线图§4.1.3 导弹追踪问题一、导弹追击问题 一艘敌舰在某海域内沿正北方向航行时,我方战舰恰位于敌舰的正西方向 1km 处.我舰向敌舰发射导弹,导弹头始终对准敌舰,敌舰速度为 v 0 =1km/min,导弹 速度为敌舰速度的 5 倍.问需要多长时间,在何处导弹击中敌舰? 以我舰位置为坐标原点,正北方向为 y 轴方向建立坐标系,设 t 时刻导弹所处 的位置为 P ( x ( t ), y ( t )) ,敌舰所处的位置为 Q (1, v 0 t ) .图 4-8 导弹追击问题由于导弹头始终对准敌舰,因此直线 PQ 是导弹运行轨迹 OP 在 P 点的切线, 即dy dx ? v0t ? y 1? x dy dx从而v 0 t ? (1 ? x ) ? y(4.1.17)又因为导弹速度为敌舰速度的 5 倍,所以 OP 弧的长度为 AQ 长度的 5 倍. 即?联立(4.1.17)和(4.1.18)得到x 0? dy ? 1? ? ? d x ? 5v0t ? dx ?2(4.1.18)?x 0dy ? dy ? 1? ? ? 5y ? d x ? 5(1 ? x ) dx ? dx ?2两边对 x 求导,得到微分方程? d y? 1? ? 1 ? ? 5 (? x ? d x?2d y ) 2 d x2(4.1.19)该问题的初始条件为y (0 ) ? 0, y '(0 ) ? 0首先,我们将二阶方程化为一阶方程组,令 y1 ? y , y 2 ?? d y1 ? dx ? y2 ? 1 / 2 ? dy ?1 ? y 22 ? ? 2 ? ? dx 5(1 ? x ) ? ? y (0 ) ? 0, y (0 ) ? 0 1 2 ? ? ?dy dx,(4.1.19)转化为(4.1.20)对于这一问题,我们可以求得方程的解析解,从而确定导弹击中敌舰位置和时 间,也可以通过求数值解的方法得到同样的结论. 1) 解析解. 关于 y 2 的方程1/ 2 ? ?1 ? y 22 ? dy2 ? ? ? dx 5(1 ? x ) ? ? y 2 (0 ) ? 0分离变量,有 dy2?1 ? y ?2 21/ 2?dx 5(1 ? x )两边积分,并代入初始条件,得ln ( y 2 ? 1 ? y2 ) ? ?21 5ln (1 ? x )即y2 ? 1 ? y 2 ? (1 ? x )2 ? 1 5(4.1.21)1又因为y2 ? 1 ? y2 ?2?1 y2 ? 1 ? y22? ? (1 ? x ) 5(4.1.22)(4.1.21)和(4.1.22)两式相加,得到y2 ? 1 2 ((1 ? x )? 1 5 1? (1 ? x ) 5 )下面,我们求解关于 y1 的方程1 1 ? ? d y1 1 ? y 2 ? ((1 ? x ) 5 ? (1 ? x ) 5 ) ? 2 ? dx ? y (0 ) ? 0 ? 1直接积分,并代入初始条件得到y ? y1 ? ? 5 84(1 ? x ) 5 ?5 126(1 ? x ) 5 ?5 24当 x ? 1 时, y ?t ? y v0 ? 1 2 .5 s5 24,即当敌舰航行到点(1,5 24)处时,被导弹击中.被击中的时间为.2) 数值解. 首先,编写出描述方程(4.1.20)的 m-文件 function dy=equ20(x,y) dy=[y(2);1/5*sqrt(1+y(2)^2)/(1-x)]; 由于方程在 x ? 1 处没有定义,我们取求解区间为[0,1-1e-8],输入命令 xf=1-1e-8; [x,y]=ode45('equ20',[0,xf],[0;0]) plot(x,y(:,1)); hold on y=0:0.01:2; plot(1,y,'*'); 运行结果 y= 0 0.0000 ?? 0.3 0 0.0001 ?? 16.8图 4-9 导弹击中曲线说明当 x ? 1 时, y ? 0.2083 ,导弹击中敌舰的位置大致在(1,0.2083)处,击中的 时间 t ?y v0=0.2083min=12.498s.这与求解析解得到的结果是一致的.比较上面 2 种求解方法,我们发现,解析解可以得到解的真实值,但是方法比较 复杂,技巧性强,只能处理比较简单的方程问题;而数值解虽然得到的解与真实值 之间会有一些允许范围内的误差,但是方法比较容易掌握,更具一般性. 二、导弹系统的改进 现根据情报,这种敌舰能在我舰发射导弹后 T 分钟做出反应并摧毁导弹.因 此,我们要求改进电子导弹系统,使其根据敌舰与我舰的距离,行使方向和速度, 能自动判断出敌舰是否在有效打击范围之内. 对于更一般的情况,我们做出如下假设,设敌舰在我舰正东方向 d km 处,行驶 速度为 v 0 km/min,行驶方向与正东方向的夹角为 ? ,导弹的飞行速度为 v km/min. 问题的关键是计算出导弹击中敌舰所需要的时间 t * ,并将 t * 与 T 比较,若 t * &T,则敌 舰在打击范围内. 我们仍以我舰位置为坐标原点,以正北方向为 y 轴建立坐标系,设 t 时刻导弹 所处的位置为 P ( x ( t ), y ( t )) ,敌舰所处位置为 Q ( d ? v 0 t co s ? , v 0 t sin ? ) .图 4-10 导弹击中问题由于导弹头始终对准敌舰,因此直线 PQ 是导弹运行轨迹 OP 在 P 点的切线, 即dy dx v v0x 0?v 0 t sin ? ? y d ? v 0 t co s ? ? x v v0(4.1.23)又因为导弹速度为敌舰速度的倍,所以 OP 弧的长度为 AQ 长度的2倍.即?? d y? 1? ? ? d x? ? d x?v0v0v t(4.1.24)联立方程(4.1.23)和(4.1.24),消去 v 0 t ,再对方程两边对 x 求导,得到微分方程d y dx2 2d y 2 1 / 2 d y 2 ( 1? ( ) ) (? ? n s i ? co s ) v dx dx ? dy dy (d ? x ) ( s ? n i ? ? o?s ) ? c o s ? (x ? (y c d dx dxv0(4.1.25)) )同样,我们令 y1 ? y , y 2 ?dy dx,(4.1.25)转化为一阶微分方程组 ? d y1 ? dx ? y2 ? v0 2 1 / 2 2 ? (1 ? ( y 2 ) ) (sin ? ? y 2 co s ? ) ? dy2 v ? ? dx ( d ? x )(sin ? ? y 2 co s ? ) ? co s ? ( y 2 ( d ? x ) ? y 1 ) ? ? y (0 ) ? 0, y (0 ) ? 0 1 2 ? ? ?(4.1.26)从图象上来判断,我们知道,当 P 和 Q 两点的运动曲线相遇时,导弹击中敌舰. 因此我们可以认为,若 x * , y * 满足方程(4.1.25)且y ? ( x ? d ) tan ?* *(4.1.27)则 点 ( x* , y* ) 为 导 弹 击 中 敌 舰 的 击 中 点 . 再 根 据 Q 点 的 表 达 式Q ( d ? v 0 t co s ? , v 0 t sin ? ) ,可以计算出击中时间 t ?*y*v 0 sin ?.若 t * &T,则敌舰在打击范围内,可以发射. 从上述分析,我们知道,当 v 0 、v、d、 ? 、T 已知时,根据(4.1.25)和(4.1.27)可 以计算出导弹击中敌舰所需要的时间 t * ,且当 t * &T 时,敌舰在打击范围内.这种对 于每组 v 0 、v、d、? 、T 分别计算从而作出判断的方法,称为在线算法.如果在敌舰 和我舰的情况已知,即 v 0 、v、T 已知,计算出所有在打击范围内的 d 和 ? ,之后要判 断敌舰是否在打击范围,只需要在计算结果中查询,这种算法称为离线算法.从原 理上来看,我们发现,在线算法灵活,容易调整参数和模型,但是速度比较慢;离线算 法事先计算好,实时使用查询方式,不需要计算,速度快. (1)在线算法 现已知 v 0 =90km/h,v=450km/h,d=50km, ? ??4,T=0.1h.首先,我们编写带参数的 m-文件描述方程(4.1.26).为了使方程始终有意义,我们通常在分母上加上一个 无穷小量 1e-8. %方程(4.1.26) function dy=equ26(x,y,v0,v,d,theta) dy(1)=y(2); dy(2)=((v0/v)*(1+(y(2))^2)^(1/2)*(sin(theta)-y(2)*cos(theta))^2) /((d-x)*(sin(theta)-y(2)*cos(theta))+cos(theta)*(y(2)*(d-x)+y(1))+1e-8); dy=dy(:); 输入命令 v0=90;v=450;d=50;theta=pi/4;T=0.1; %初始化参数 [x,y]=ode45(@equ26,[0,60],[0;0],[],v0,v,d,theta); %计算导弹运行轨迹方程的数值解 z=(x-d)*tan(theta); %计算敌舰运行的轨迹 n=length(x); for i=1:n if abs(z(i)-y(i,1))&1e-6 xk=x(i); yk=y(i,1); end end xk yk %求出击中点xk,yk坐标 for i=1:n if z(i)&0 z1(i)=z(i); else z1(i)=0; end end plot(x,y(:,1),x,z1(:),xk,yk,'*'); legend('导弹运行轨迹','敌舰运行轨迹','击中点'); tk=yk/(v0*sin(theta)+1e-8) %计算击中时间 if tk&T disp(['敌舰在打击范围内，击中地点在', num2str(xk),num2str(yk),'击中时间为', num2str(tk)]); else disp(['敌舰不在打击范围内']); end 得到结果 击中点 xk=58.4056,yk=8.4056,击中时间 tk =0.1321,所以敌舰不在打击范围内,应 等接近一些再发射. 图4-11 在线算法(2)离线算法 在介绍离线算法之前,我们先来看方程(4.1.26)的另一种表达.因为导弹的线 速度v ? ? dx ? ? dy ? ? ? ?? ? ? dt ? ? dt ?2 2(4.1.28)方程(4.1.28)与(4.1.23)联立,可以建立以 t 为参数的关于 x,y 的参数方程? dx ? ? dt ? ? ? ? ? dy ? ? dt ? ? ? v ? dy ? 1? ? ? ? dx ? v ? dx ? 1? ? ? ? dy ?2 2?v ? v 0 t sin ? ? y ? 1? ? ? ? d ? v 0 t co s ? ? x ?2?v ? d ? v 0 t co s ? ? x ? 1? ? ? ? v 0 t sin ? ? y ?2(4.1.29)当 x ( t ) 满足 x ( t ) ? d ? v0 t c o s? ,则导弹击中了敌舰.由于我们所要求的打击范围是 在 t ? T 中讨论的,所以考虑以 t 为参数的方程(4.1.29)能使求解过程大大简化. 现已知 v 0 =90km/h,v=450km/h,T=0.1h,要计算出所有在打击范围内的 d 和 ? . 依题意, 0 ? ? ? ? .据此,我们可以确定 d 的取值范围.当 ? ? 0 时,敌舰正好背向行驶, 导弹直线运行,击中时间t ? d /( v ? v 0 ) ? T求得 d m in ? T ( v ? v 0 ) ? 36 km .当 ? ? ? 时,敌舰迎面驶来,导弹直线运行,击中时间 t ? d /( v ? v 0 ) ? T则 d m ax ? T ( v ? v 0 ) ? 54 km .所以 3 6 ? d ? 5 4 .这样,我们对于所有可能的 d 和 ? 的 取值,计算击中所需时间,从而对不同的 ? ,得到 d 的临界值.具体应用时直接查 询判断. 编写如下 m-文件描述方程(4.1.29) %方程(4.1.29) function dy=equ29(t,y,v0,v,d,theta) dydx=(v0*t*sin(theta)-y(2))/(d+v0*t*cos(theta)-y(1)+1e-8); dy(1)=v/(1+dydx^2)^(1/2); dy(2)=v/(1+dydx^(-2))^(1/2); dy=dy(:); 输入命令 v0=90;v=450;d=50;theta=pi/4;T=0.1; i=1; for d=54:-1:36 for theta=0:0.1:pi [t,y]=ode45(@equ29,[0,T],[0;0],[],v0,v,d,theta); if max(y(:,1)-d-v0*t*cos(theta))&0 range(i,:)=[d,theta]; i=i+1; plot(range(:,1),range(:,2)); xlabel('d');ylabel('theta'); 运行得到临界曲线,即在该曲线上方的 d 和 ? 值,所对应的敌舰位置在打击范 围内,曲线下方不在打击范围内. 图4-12 离线算法因此 d=50km, ? ??4,不在打击范围内.§4.1.4 微分方程稳定性理论简介在处理实际问题时,对于有些微分方程模型我们不仅要得到问题的解,有时 还需要研究解的稳定性,即解对初始值的连续依赖性,如果解在一定范围内是稳 定的,那么初始条件发生一些小的扰动(如实验测量误差等),对问题的解不会造 成影响.另外,还有一些问题我们并不需要求解,而通过解的变化趋势的研究,并 分析一些特殊解的稳定性就可以解决问题.本节仅介绍几类特殊,但常用的常微 分方程稳定性分析的基本方法. 一、单个常微分方程的平衡点及稳定性 若微分方程dx dt ? f (x)(4.1.30)方程右端不显含自变量 t ,称为自治方程.代数方程的实根 x ? x 0 称为方程(4.1.30) 的平衡点(或奇点).注意到,平衡点也是方程的解(奇解). 如果从一定范围内的初始条件出发,方程(4.1.30)的解 x (t ) 都满足t ? ??lim x ( t ) ? x 0 ,则称平衡点 x 0 是稳定的. 对于一些不易求解的问题,我们可以不求方程(4.1.30)的解,不用定义来判断 平衡点 x 0 的稳定性,下面我们来介绍这种方法. 将 f ( x ) 在 x 0 点作泰勒(Taylor)展开,只取一次项,得到方程(4.1.30)的近似线 性方程dx dtx0? f ' ( x 0 )( x ? x 0 ),0(4.1.31))t|也是方程(4.1.31)的平衡点,方程(4.1.31)的通解为 x ( t ) ? ce | f '( x? x 0 .关于 x 0点稳定性有如下结论: ① 若 f ' ( x 0 ) ? 0 ,则 x 0 对于方程(4.1.31)和(4.1.30)都是稳定的; ② 若 f ' ( x 0 ) ? 0 ,则 x 0 对于方程(4.1.31)和(4.1.30)都是不稳定的. 二、二阶常微分方程组的平衡点及稳定性 现在讨论二阶微分方程组 ? dx ? f ( x , y ), ? dt ? ? ? d y ? g ( x , y ). ? dt ?(4.1.32)它的解x ? x ( t ), y ? y ( t )在以 t , x , y 为坐标的欧氏空间中决定了一条曲线,如果把时间 t 看作参数,仅考虑x, y为坐标的欧氏空间,此空间称为方程(4.1.32)的相平面(若方程组是高阶,则称为相空间).对于右端函数不显含时间 t 的自治系统? dx ? dt ? f ( x , y ), ? dy ? ? g ( x , y ). ? dt(4.1.33)代数方程组? f ( x, y ) ? 0, ? ? g ( x, y ) ? 0.的 实 根 x ? x 0 , y ? y 0 称 为 方 程 (4.1.33) 的 平 衡 点 ,记 作 P0 ( x 0 , y 0 ) . 它 也 是 方 程 (4.1.33)的解. 如果从一定的范围内的初始条件出发,方程(4.1.33)的解 x ( t ), y ( t ) 都满足t ? ??lim x ( t ) ? x 0 , lim y ( t ) ? y 0 ,t ? ??则称平衡点 P0 是稳定的,否则称 P0 是不稳定的. 与单个方程的讨论类似,对于不易求解的微分方程组,我们可以通过研究与 其对应的近似线性方程组,从而得到平衡点和稳定性.在这里,我们省略证明过程, 只给出判别平衡点 P0 是否稳定的判别准则.令? f ( P0 ) ? ? f ( P0 ) ? g ( P0 ) ? ?x p ? ?? ? ?, q ? ?g ( P ) ?y ? 0 ? ?x ?x ? f ( P0 ) ?y , ? g ( P0 ) ?y关于 P0 点稳定性有如下结论: ① 当 p ? 0 且 q ? 0 时,平衡点 P0 是稳定的; ② 当 p ? 0 或 q ? 0 时,平衡点 P0 是不稳定的. §4.1.5 动物种群模型分析一、单一种群模型 (一) Malthus 增长模型 英国经济学家和人口统计学家 Malthus()根据百余年的统计资料, 在 1798 年提出了闻名于世的人口指数增长模型(Malthus 人口模型).此模型也适用 于单种族问题.此模型有两个基本假设前提： ① 假设人口数量 x (t )是时间 t 的连续可微函数, 且 x (0) = x0. 虽然人口数量 总是离散变化的, 但实际上增加一个人所引起的变化与庞大的人口数量相比是 微不足道的. 这里实际上是将离散变量作连续化处理, 在差分方程中则恰好相 反. ② 人口数量的增长速度与现有人口数量成正比, 比例系数为 r. 这个假设在 当时是比较合理的. 在此基础上, Malthus 提出了如下的人口模型(即指数增长模型)：dx dt= rx,x (0) = x0 .不难求得其解为 x (t ) = x0ert, 这个模型比较准确地反映了在
年期间的 人口总数. (二)Logistic 增长模型 因为 r＞0, x0＞0, x (t ) = x0ert 表明人口数量 x (t ) 将随时间 t 按指数规律无限 地增长, 对于有限的地球资源来说, 这当然是不合理的. 即假设②在自然环境条 件较好, 人口密度较小时是合理的, 而对于人口密度较大时它是不合理的, 这是 因为人类生存空间及可利用资源(食物、水、空气)等环境影响对人口的增长起着 阻滞作用. 从这里我们可以看出, 在进行数学建模时, 我们必须作合理的假设. 现在我们对指数增长模型加以改进.设人类生存空间及可利用资源(食物、 水、 空气)等环境因素所能容纳的最大人口容量为 K (称为饱和系数). 人口数量 x 的增 长速率不仅与现有人口数量成正比, 而且还与人口尚未实现的部分(相对最大容 量 K 而言)所占比例(K ??x)/K 成正比, 比例系数为固有增长率 r. 修改后的模型为? dx ?K ? x? ? rx ? ?, ? ? K ? ? dt ? x (0) ? x . 0 ?Kx不难求得其解为 x ( t ) ?0x 0 ? ( K ? x 0 ) e ? rt.这就是非常著名的 Logistic 增长模型(即阻滞增长模型), 它是荷兰数学家、生物学家 Verhulst 在 1839 年首次提出的.因子(K ??x)/K 的生物学含义是“剩余空间”或尚未利用的增长机会. 图 4-13 描绘了dx dt与 x 之间的关系, 它表明人口变化率dx dt随着人口数量 x 的 增加而先增后减, 在 x = K / 2 处到达最大值, 在 x = K 处, Logistic 曲线 x (t ), 且有 lim 义看这是明显的结果.t? ?dx dt= 0. 图 4-14 描绘了x (t ) ? K, 即 x = K 是稳定的平衡点, 从模型本身的意图 4 -13图 4 -14Logistic 模型的曲线Logistic 模型用途十分广泛, 除了用于预测人口增长外, 也可完全类似地用 于虫口增长、疾病的传播、谣言的传播、技术革新的推广、销售预测等.但这个 模型的最大缺点是饱和系数 K 不易确定. 二、两种物种相互竞争与相互依存模型当某个自然环境中只有一种生物的群体(生态学上称为种群)生存时, 人们常用 Logistic 模型来描述这个种群数量的演变过程.即dx x ? ? ? rx ? 1 ? ? dt K ? ?x(t ) 是种群在时刻 t 的数量, r 是固有增长率, K 是环境资源容许的种群最大数量, 且 x = K 是 稳定的平衡点. 如果一个自然环境中有两个或两个以上种群生存, 那么它们之间就要存在着或是相互 竞争, 或是相互依存, 或是弱肉强食(捕食--被捕食模型)的关系. ⑴ 种群的相互竞争 当两个种群为了争夺有限的同一种食物来源和生活空间而进行生存竞争时, 最常见的 结局是竞争力较弱的种群灭绝, 竞争力较强的种群达到环境容许的最大数量. 设有甲乙两个种群, 当他们独自在一个自然环境中生存时, 数量的演变均遵从 Logistic 规律, x1 (t ), x2 (t )是两个种群的数量, r1, r2 是它们固有增长率, N1, N2 是它们的最大容量. 于 是对种群甲有d x1 ? x1 ? r1 x 1 ? 1 ? ? dt N1 ? ? ? ? ?, ? ?x1 N1其中因子 ? 1 ? ?x1 ? ? 反映由于甲对有限资源的消耗导致对它本身增长的阻滞作用, ? N1 ?可解释为相对于 N1 而言单位数量的甲消耗的供养甲的食物量(设食物总量为 1). 当两个种群在同一自然环境中生存时, 考察由于乙消耗同一种有限资源对甲的增长产 生的影响, 可以合理地在因子 ? 1 ? 1 ? 中再减去一项, 该项与种群乙的数量 x2 (相对于 N2 而 ? ? N1 ? ? 言)成正比, 得到种群甲的增长方程d x1 ? x1 x2 ? r1 x 1 ? 1 ? ??1 ? dt N1 N2 ? ? ?, ? ??x?这里的?1 意义是, 单位数量乙(相对 N2 而言)消耗的供给甲的食物量为单位数量甲(相对 N1) 消耗的供给甲的食物量的?1 倍. 类似地, 得到种群乙的增长方程dx 2 ? ? r2 x 2 ? 1 ? ? ? dt ? x12?N1x2 ? ?, ? N2 ?对?2 可作相应地解释. 关于种群的相互竞争模型? d x1 ? x1 x2 ? ?, ? r1 x 1 ? 1 ? ??1 ? ? ? N1 N2 ? ? dt ? ? ? x1 x2 ? ? dx 2 ?. ? r2 x 2 ? 1 ? ? 2 ? ? ? ? dt N1 N2 ? ? ?的稳定性分析见表 4-1. 表 4-2 种群相互竞争模型的平衡点及其稳定性 平衡点 P1 (N1, 0) P2 (0, N2 ) P3 ( N1 (1???1) N2 (1???2) , ) 1???1?2 1???1?2 P4 (0, 0) p r1 ??r2 (1???2) r2 ??r1 (1???1) r1 (1???1) + r2 (1???2) 1???1??2 ??(r1 + r2) q ??r1 r2 (1???2) ??r1 r2 (1???1) r1 r2 (1???1) (1???2) 1???1??2 r1 r2 稳定条件?1＜1,????2＞1 ?1＞1,????2＜1 ?1＜1,????2＜1不稳定根据表 4-2 和?1、??2 的含义, 稳定的平衡点 P1, P2, P3 在生态学上的意义解释如下： ① 当?1＜1,????2＞1 时, ?1＜1 意味着在对供养甲的资源的竞争中乙弱于甲, ??2＞1 意味 着在对供养乙的资源的竞争中甲强于乙, 于是种群乙终将灭绝, 种群甲趋向最大容量, 即 x1(t)、x2 (t)趋向平衡点 P1 (N1, 0). ② 当?1＞1,????2＜1 时, 与①正好相反. ③ 当?1＜1,????2＜1 时, 因为在竞争甲的资源中乙较弱, 在竞争乙的资源中甲较弱, 于 是可以达到一个双方共存的稳定的平衡状态 P3. ④ 当?1＞1,????2＞1 时, P1, P2, P3 都是稳定的平衡点, 需根据初始条件确定, 这是种群竞 争中很少出现的情况, 因为在竞争甲的资源中乙较强, 在竞争乙的资源中甲较强. ⑵ 种群的相互依存模型 自然界中处于同一环境下两个种群相互依存而共生的现象是很普遍的, 其类型可分为 三种：①其中一种能独立生存, 另其中一种则不能独立生存；②两个种群都能独立生存；③ 两个种群都不能独立生存.下面介绍类型①, 另两种留作习题 7. 设种群甲能独立存在, 按 Logistic 规律增长, 种群乙为甲提供食物, 有助于甲的增长.例 如, 植物可以独立生存, 昆虫的授粉作用又可以提高植物的增长率, 而以花粉为食物的昆虫 却不能离开植物单独存活.种群甲的数量演变规律为d x1 ? x1 x2 ? r1 x 1 ? 1 ? ??1 ? dt N1 N2 ? ? ?, ? ?这里 x1 (t), x2 (t), r1, N1, N2 的意义同前, ?1 表示单位数量乙(相对 N2 而言)提供的供养甲的食物 量为单位数量甲(相对 N1)消耗的供给甲的食物量的?1 倍. 种群乙没有甲的存在会灭亡, 设其死亡率为 r2, 甲为乙提供食物对乙的增长有促进作用, 同时乙的增长又会受到自身的阻滞作用, 因此种群乙的数量演变规律为dx 2 ? ? r2 x 2 ? ? 1 ? ? ? dt ? x12?N1x2 ? ?, ? N2 ?这里的?2 表示单位数量甲(相对 N1 而言)提供的供养乙的食物量为单位数量乙(相对 N2)消耗 的供给乙的食物量的?2 倍. 对于种群的相互依存模型? d x1 ? x1 x2 ? ?, ? r1 x 1 ? 1 ? ??1 ? ? ? N1 N2 ? ? dt ? ? ? x1 x2 ? dx 2 ? r2 x 2 ? ? 1 ? ? 2 ? ? ? dt N1 N2 ? ?? ?. ? ?的稳定性分析见表 4-3. 表 4-3 种群相互依存模型的平衡点及其稳定性 平衡点 P1 (N1, 0) P3 ( N1 (1???1) N2 (?2 ??1) , ) 1???1?2 1???1?2 P4 (0, 0) p r1 ??r2 (?2 ??1) r1 (1???1) + r2 (?2 ??1) 1???1??2 ??r1 + r2 q ??r1 r2(?2 ??1) r1 r2 (1???1) (?2 ??1) 1???1??2 r1 r2 稳定条件?2＜1, ?1?2＞1 ?2＞1,???1?2＜1不稳定根据表 4-3 和?1、?2 的含义, 稳定的平衡点 P1, P2 在生态学上的意义解释如下： ① 当?2＜1,???1?2＜1 时, ?2＜1 意味着甲供养乙的资源不足以维持乙的生存, 于是种群 乙终将灭绝, 种群甲趋向最大容量, 即 x1 (t)、x2 (t)趋向平衡点 P1 (N1, 0). ② 当?2＞1,???1?2＜1 时, ?2＞1 意味着甲供养乙的资源能足以维持乙的生存, ?1?2＜1 即 有?1＜1, 意味着对乙供养甲的资源加以限制, 防止甲过分地增长, 于是可以达到一个双方 共存的稳定的平衡状态 P2. ③ 当?2＜1,???1?2＞1 时, 基本上同①. ④ 当?2＞1,???1?2＞1 时, 需要适当对它们加以限制. ⑶ 种群的弱肉强食模型 处于同一自然环境中两个种群之间的关系出了相互竞争和相互依存之外, 还有一种更 为有趣的生存方式： 种群甲靠丰富的天然资源生长, 而种群乙则靠掠食甲为生.本世纪二十年 代, 意大利著名生物学家 D’Ancona 在研究时无意中发现, 第一次世界大战期间从地中海捕 获的鱼中, 鲨鱼等食肉鱼所占的比重十分明显地上升了. D’Ancona 认为, 这一现象决非偶然, 它是由战争期间捕鱼量减少造成的.可是, 为什么捕鱼量减少, 食用鱼反而减少呢？他无法 解释这个现象, 于是就找了著名数学家 Vloterra, 希望他能建立一个数学模型来解释之. Vloterra 根据建模目的将鱼分为两大类：食用鱼 (被捕食者) 和食肉鱼(捕食者), 用 x (t) 和 y (t )分别表示 t 时刻两者的数量.他分析, 假设如果没有食肉鱼, 食用鱼以净自然增长率 r ? (r＞0)的速度无限增长, 即 x /x = r. 但是, 食肉鱼以食用鱼为食物, 致使鱼的增长率降低, 设 ? 降低的程度与食用鱼数量 y (t )成正比. 于是, x / x = r ????y. 常数? 反映食肉鱼掠取食物的 ? 能力. 如果没有食用鱼, 食肉鱼无法生存. 设食肉鱼的自然死亡率为 d (d＞0), 则 y / y= ??d. 食用鱼为食肉鱼提供了食物, 致使食肉鱼的死亡率降低, 或者说为食肉鱼提供了增长的条件. ? 设增长率与食用鱼的数量 x (t )成正比, 于是, y / y = ?d + ? x, 常数 ??反映食用鱼对食肉鱼 的供养能力. 根据上述分析可得该模型为? dx ? rx ? ? xy , ? ? dt ? ? d y ? ? dy ? ? xy . ? dt ?4-4这就是著名的捕食--被捕食模型. 容易看出它的两个平衡点是 A (0, 0) 和 B (x0, y0), 其中 x0 = d /??, y0 = r/ ? . 对于平衡点 A(0, 0), p = r ??d, q = ??rd＜0,所以 A (0, 0)不稳定；对于平衡点 B (x0, y0), p = 0, q = rd＞0, 处 于临界状态, 不能用平衡点稳定性判别准则判别. 下面通过求出该模型的解来进行分析. 从微分方程组(4-1)中消去 d t 后, 用分离变量法得到其通解为 xde???x yre???y = c 4-5 式中的 c＞0, 为任意常数, 可由初始条件确定. 下面再来研究方程(4-5)的图形(通常称之为 微分方程组(4-4)的轨线或相轨线). 记 f (x) = xde???x, g (y) = yre?? y, 利用微积分中函数作图的方法作出它们的图形 (图 4-3, 图 4-4), 函数的极值点分别为 x = x0, y = y0, 函数的极大值分别记为 f m 和 gm.图 4 -3图 4 -4当 c＞f m gm 时, 方程(4-5)无解； 当 c = f m gm 时, 方程(4-5)有唯一解：x = x0, y = y0,轨线退化为一点 B (x0, y0) (平衡点)； 当 c＜f m gm 时, 设 c = ??gm (0＜? ＜ f m ＝, 令 y = y0, 注意到 g ( y0) = gm, 由方程(4-5) 可得代数方程 f (x) = ??, 从图 4-3 知道这个方程有两个根 x1 和 x2, 且 x1＜x0＜x2, 于是在图 4-5 上得到这条轨线上的两个点 S1 (x1, y0), S2 (x1, y0). 对于区间 [x1, x2 ] 之外的任意一点 x, 从图 4-3 知 f (x)＜?, 再由 f (x) g ( y) = ? gm 得 g ( y) ＞gm, 与 gm 的定义矛盾. 对于区间(x1, x2)内的任意一点 x, 因为 f (x)＞?,图 4 -5 f (x) g ( y) = ? gm, 所以 g ( y)＜gm . 从图 4-4 知道方程 g ( y) = ???( ? ＜gm )有两个根 y1 和 y2, 且 y1＜ y0＜y2, 于是在图 4-5 上得到同一条轨线上的另外两个点 S3 和 S4. 因为 x 是任意的, 所 以可知这条轨线是封闭曲线. 通过上面的分析可知, 当 c 从 f m gm 变小而趋于零时, 轨线是一族以平衡点 B (x0, y0)为中 ? ? 心, 越来越扩展的封闭曲线(图 4-6). 轨线的方向可由图上四个区域中 x 和 y 的符号决定.图 4 -6 方程(1)的轨线族图 4 -7封闭轨线对应着方程(4-4)的周期解(图 4-7) x = x (t), y = y (t), 设周期为 T, 利用方程(4-4) 可以求出 x (t ) 和 y (t )在一周期内的平均值 x , y .x ? 1 T? x (t )d t ? T ?0T1T(0d??1 dy y dt) d t = d / ??+ [ ln y (T ) ??ln y (0)]/T = d / ??,类似地有 y = r /??. 即 x (t ) 和 y (t )在一周期内的平均值 x , y 为它们在平衡点 B (x0, y0)的值. 对于周期性变化的 x (t ) 和 y (t ), 我们用它们的平均值 x , y 表示其大小.综上所述, 食肉 鱼的数量取决于两个参数 r 和??. 当食用鱼的自然增长率 r 下降时, 食肉鱼的数量将减少, 这 就是说, 在弱肉强食情况下, 如果人们想使强者减少, 应当去降低弱者的繁殖率, 而当食肉 鱼掠取鱼的能力? 提高时, 对食用鱼没有影响, 只会使食肉鱼减少.另一方面, 食肉鱼死亡率 d 上升将导致食用鱼的增多, 而食用鱼对食肉鱼供养能力 ??的提高会使食用鱼减少. 下面用这个模型解释开始提出的为什么战争时捕捞量下降食用鱼反而减少的问题. 上面我们讨论的是在自然环境中食用鱼和鲨鱼数量的变化规律, 在人工捕鱼的情况下, 设表示捕捞能力的系数为 e, 这相当于食用鱼的自然增长率由 r 降为 r ??e, 鲨鱼的死亡率由 d 升为 d + e.在这种改动下上述模型仍然成立, 于是食用鱼和鲨鱼的(平均)数量 x , y 变为： x = (d + e)/?, y = (r ??e)/??. 因此捕捞能力 e 下降, 鲨鱼的数量( y )增加, 而食用鱼的数量( x )减少. 上述分析用于解释使用杀虫剂的影响时也是出乎人们意料的.害虫靠吃农作物生存, 为 被捕食者；益虫是害虫的天敌, 为捕食者.如果某种杀虫剂不仅杀死害虫, 也能杀死益虫, 那 么增强杀虫能力相当于提高上述分析中的捕捞系数 e, (从长期效果上看)会使害虫( x )增多, 益虫( y )减少, 事与原违！这个事实说明, 捕食--被捕食模型需要改进, 改进后的模型如下：? d x1 ? x1 x2 ? ?, ? r1 x 1 ? 1 ? ??1 ? ? ? N1 N2 ? ? dt ? ? ? x1 x2 ? dx 2 ? r2 x 2 ? ? 1 ? ? 2 ? ? ? dt N1 N2 ? ?? ?. ? ?各符号的含义及其稳定性分析留作习题 8.第二节 差分方程模型 在上一节中,我们利用微分方程方法研究了一些连续变化的变量,但现实世 界中常常会遇到一些问题,其变量只在一定的离散集合中取值,例如,在对经济进 行动态分析时,一般总根据一计划周期期末的指标值来判断某经济计划执行得如 何.有些变量虽然可以认为是连续变化或近似连续变化得,但在研究时,并不时时 刻刻去统计它,而只在某些特定时刻才去统计它.这实际上就是差分的基本思想, 利用这种方法,我们可以将微分方程的变量离散化,得到相应的差分方程模型.§4.2.1 差分方程模型及求解这一节中,我们先介绍差分及差分方程的一些基本概念,然后就几类特殊的 差分方程的解法进行讨论.§4.2.1.1 基本概念以 n 表示时间,规定 n 只取非负整数,记 x n 为变量 x 在 n 时刻的取值,则称? x n ? x n? 1 ? x n 为 x n的一阶差分,称 ? 2 x n ? ? ( ? x n ) ? ? x n ?1 ? ? x n ? x n ? 2 ? 2 x n ?1 ? x n 为 阶差分 ? k x n .x n 的二阶差分.类似地,可以定义 x n 的 k由 n , x n 以及 x n 的差分给出的方程称为 x n 的差分方程,其中含有 x n 的最高阶 差分的阶数称为差分方程的阶.差分方程也可以写成不显含差分的形式.例如,二 阶差分方程 ? 2 x n ? ? x n ? x n ? 0 ,可以改写成 x n ? 2 ? x n ? 1 ? x n ? 0 .因此,k 阶差分方程 也可以写成F ( n , x n , x n ?1 , ? , x n ? k ) ? 0(4.2.1)满足差分方程(4.2.1)的序列 x n ? x ( n ) ,称为差分方程(4.2.1)的解.包含 k 个任 意常数的解称为(4.2.1)的通解 , x 0 , x1 , ? , x k ?1 为已知时,称为(4.2.1)的初始条件 , 通 解中的任意常数都由初始条件确定后的解称为(4.2.1)的特解. 形如x n ? k ? a1 ( n ) x n ? k ? 1 ? ? ? a k ( n ) x n ? f ( n )(4.2.2)的差分方程,称为 k 阶线性差分方程.其中 a i ( n ) 为已知系数,且 a k ( n ) ? 0 . 若差分方程(4.2.2)中 f ( n ) ? 0 ,则称差分方程(4.2.2)为 k 阶齐次线性差分方程, 否则称为 k 阶非齐次线性差分方程.§4.2.1.2 差分方程的求解 一、常系数线性齐次差分方程的解析解法 形如x n ? k ? a1 x n ? k ?1 ? ? ? a k x n ? 0(4.2.3)的差分方程,称为 k 阶常系数线性齐次差分方程,其中 a i 为常数,且 a k ? 0 .代数方程? ? a1 ?k k ?1? ? ? ak ? 0(4.2.4)称为差分方程(4.2.3)的特征方程,其根称为特征根. 根据特征根的情况,我们可以得到方程(4.2.3)的通解. 1) 特征值均为单根. 若差分方程(4.2.3)的特征方程(4.2.4),有 k 个相异的特征根 ?1 , ? 2 , ? , ? k ,则x n ? C 1 ?1 ? C 2 ? 2 ? ? ? C k ? kn n n是方程(4.2.3)的一个通解,其中 C i 为任意常数. 如果有一组初始条件x 0 ? u 0 , x1 ? u 1 , ? , x k ? 1 ? u k ? 1则可以确定出一个满足初始条件的特解. 2) 特征值中有重根 若差分方程(4.2.3)的特征方程(4.2.4),有 l 个相异特征根 ?1 , ? 2 , ? , ? l ,重数依次 为 m 1 , m 2 , ? m l .显然 ? m i ? k .则差分方程(4.2.3)的通解为i ?1 lm1 m2 mixn ??C1 j nj ?1?1 ?nj ?1?C2 jnj ?1?2 ? ? ?nj ?1?Cj ?1ijnj ?1?in3) 特征值中有复根 若 差 分 方 程 (4.2.3) 的 特 征 方 程 (4.2.4) 的 特 征 根 中 出 现 一 对 共 轭 复 根?1, 2 ? ? (co s ? ? i sin ? ) ,和相异的 k ? 2 个根 ? 3 , ? 4 , ? , ? k ,则差分方程(4.2.3)的通解为x n ? C 1 ? co s n? ? C 2 ? sin n? ? C 3 ? 3 ? ? ? C k ? kn n n n.注 1:实际上,如果我们把复根的形式用指数表示,即 ?1 ? e ? ? i? , ? 2 ? e ? ? i? .则可以直 接运用情况 1)全为单根的相应结论得到同样的通解,即x n ? C 1 ?1 ? C 2 ? 2 ? ? ? C k ? k ? C 1 en n n n n n ? ? in?? C 2enn ? ? in?? C 3?3 ? ? ? C k ? kn nn? C '1 ? sin n? ? C ' 2 ? co s n? ? C ' 3 ? 3 ? ? ? C ' k ? k其中 C 'i 在复数集中取值. 同理可得,如果出现特征值为复根且为重根的情形,我们可以用同样的方法,将 复根化成指数表达的形式,然后运用情况 2)的结论得到通解. 例 4-11 求 Fibonacci 数,通项为 Fn ? 2 ? Fn ? 1 ? Fn 的通解,并确定满足初始条件F1 ? F 2 ? 1 的特解.1? 2 5 1? 2 5解:该差分方程的特征方程为 ? 2 ? ? ? 1 ? 0 .特征根为 ?1 ?, ?2 ?,是互异的,根据情形 1)的结论,我们可以得到,该差分方程的通解为Fn ? C 1 ( 1? 2 5 ) ? C2(n1? 25)n将初始条件 F1 ? F2 ? 1 代入通解,得C1 ( 1? 2 5 ) ? C2( 1? 2 5 ) ? 1, C 1 ( 1? 2 5 ) ? C2(21? 25) ?12所以 C 1 ?1 5,C2 ? ?1 5,则满足初始条件得特解为1 ? 1? 5 n 1? 5 n? ) ?( ) ? ?( 2 2 5 ? ?Fn ?例 4-12 求差分方程 x n ? 4 ? x n ? 3 ? 3 x n ? 2 ? 5 x n ? 1 ? 2 x n ? 0 的通解,并求满足初始条件x 0 ? 1, x1 ? 0, x 2 ? 1, x 3 ? 2的特解.解:该方程对应的特征方程为? ? ? ? 3? ? 5 ? ? 2 ? 04 3 2特征根为-1,-1,-1,2,出现重根的情形,利用 2)的结论,得方程的通解为x n ? ( C 1 ? C 2 n ? C 3 n )( ? 1) ? C 4 22 n n将初始条件代入通解,得到 C 1 ? 的特解为xn ? 1 5242 52,C2 ??29 52, C3 ?7 52,C4 ?10 52,所以满足初始条件( 4 2 ? 2 9 n ? 7 n )( ? 1) ?2 n10 522n例 4-13 求差分方程 x n ? 2 ? x n ? 0 的通解,并求满足初始条件 x 0 ? 1, x1 ? 2 时的特解. 解:该方程是齐次线性方程,其对应的特征方程为 ? ?1 ? 02对应的特征根为一对共轭复根 ?1 ? i ? co sn? 2?2? i sin?2, ? 2 ? ? i ? co s?2? i sin?2.所以根据情形 1)和 3)的结论,我们都可以得到,该差分方程的通解为x n ? C 1 co s ? C 2 sin n? 2当 x 0 ? 1, x1 ? 2 时,代入 x n 的表达式,求得 C 1 ? 1, C 2 ? 2 ,所以满足初始条件的特解为x n ? co s n? 2 ? 2 sin n? 2二、系数线性非齐次差分方程的解析解法 形如x n ? k ? a1 x n ? k ? 1 ? ? ? a k x n ? f ( n )(4.2.5)的差分方程,称为 k 阶常系数线性齐次差分方程,其中 a i 为常数,且 a k ? 0 , f ( n ) ? 0 . 差分方程(4.2.5)对应的齐次差分方程为(4.2.3). 与微分方程类似,非齐次差分方程的解也满足叠加原理.即非齐次差分方程 (4.2.5)的通解等于对应的齐次差分方程(4.2.3)的通解加上非齐次方程的特解. 非齐次差分方程(4.2.5)特解的求解方法,与非齐次方程微分方程类似,可以通 过 常 数 变 异 法 , 对 一 些 特 殊 的 f (n) 函 数 , 还 可 以 用 比 较 系 数 法 . 例 如 , 当f (n ) ? b pl (n )n时,其中 p l ( n ) 为 n 的 l 次多项式.1)若 b 不是特征根,则非齐次方程(4.2.5)有形如 b n q l ( n ) 的特解, 其中 q l ( n ) 为 n 的 l 次多项式; 2)若 b 是 r 重特征根,则非齐次方程(4.2.5)有形如 b n n r g l ( n ) 的特解, 其中 g l ( n ) 为 n 的 l 次多项式; 进而,通过待定系数法,将带有待定未知常数的特解代回原方程,可以求出q l ( n ) , g l ( n ) ,从而求得非齐次方程的特解.另外,由于特解的随意性,只要能找到一个满足方程的特殊解就可以解决问题, 所以有时对于一些比较简单的差分方程,我们可以直接通过观察得到特解. 例 4-14 求解差分方程 x n ? 1 ? 2 x n ? 1 ,满足初始条件 x 0 ? 1 的解. 解:该方程是一个简单的非齐次问题,其对应的齐次方程为x n ?1 ? 2 x n ? 0对应的特征方程为 ? ? 2 ? 0 ,特征根为 ? ? 2 ,所以齐次方程的通解为xn ? C1 2n通过观察我们发现,如果我们取 x n ? a 始终为常数,代入原差分方程可以找到 a ? 2 a ? 1 ,即 x n ? a ? ? 1 ,为满足非齐次方程的一个特解.根据叠加原理,我们得到非齐次方程的特解为xn ? C1 2 ? 1n代入初始条件 x 0 ? 1 ,得到 C 1 ? 2 .所以满足初始条件的原方程的解为xn ? 2n ?1?1例 4-15 求差分方程 x n ? 2 ? 4 x n ? 1 ? 4 x n ? 2 n 的通解. 解: 该方程为非齐次差分方程,其对应的齐次方程的特征方程为? ? 4? ? 4 ? 02求得特征根为二重根 ?1 ? ? 2 ? 2 ,所以齐次方程的通解为xn ? C1 2 ? C 2 n 2n n根据题述非齐次差分方程的右端 f ( n ) ? 2 n ,又因为 2 是特征方程的二重根,可得非 齐次方程有形如xn ? C 3n 22 n的特解.将 x n 代入原方程,解得 C 3 ?n1 8n.所以非齐次方程的通解为n n 2 n?3xn ? C1 2 ? C 2 n 2 ? xn ? C1 2 ? C 2 n 2 ? n 2三、计算机求解 若 x0, x1, ? , xk?1 已知,则形如 xn+k = g (n; xn, xn+1, ? , xn+k?1 )的差分方程的解可 以在计算机上直接用循环语句实现.这种方法虽然不能求得差分方程的解的具体 表达式,但是可以知道每个时刻的数值解,而且适用于一切有显示表达的差分方 程. 例 4-16 计算差分方程xn ? 4 ? x2 n?3? 3 x n ? 2 ? 5 x n ?1 ? 2 x n ? 2n满足初值条件 x 0 ? 1, x1 ? 0, x 2 ? 1, x 3 ? 2 的解的前 10 项. 解:该差分方程含有非线性项,根据我们之前讨论的求解方法,不能处理这类问题, 我们只能通过通项求解差分方程的有限项的数值解. 在 matlab 中输入命令 &&x(1)=1;x(2)=0;x(3)=1;x(4)=2; %赋初值,注意,matlab 中的数组中元素的下标从 1 开始 k=4; n=10; %差分方程阶数 %求解项数 for i=k+1:n x(i)=-(x(i-1))^2+3*x(i-2)+5*x(i-3)+2*x(i-4)+2^(i-4); %通项公式 end x 运行结果: x= 1 0 1 2 3 6 -7 4 31 -908 即为所求差分方程解的前 10 项. 对于我们之前求解的常系数差分方程,也可以用此方法求解数值解.例如,例 4-14 中的差分方程,我们将上述命令稍做改动,便可以得到它的解的前 20 项. 编写命令 &&x(1)=1; k=1; n=20; for i=k+1:n x(i)=2*x(i-1)+1; end x 运行结果: x = 1 3 7 15 31 63 127 255 511 95 8191 %输出结果 %通项公式 %差分方程阶数 %求解项数 %输出结果
287 1048575 运用我们在例 4-14 中得到的解的通项公式 x n ? 2 n ? 1 ? 1 ,可以计算出第 20 项x1 9 ? 220? 1 ? 1 0 4 8 5 7 5 ,与我们的数值解是一致的.§4.2.2 差分方程稳定性理论简介与微分方程模型类似,差分方程平衡点的稳定性研究也是有重大意义的.在 这一节中,我们先给出差分方程平衡点及其稳定性的概念,再对几类特殊的差分 方程的平衡点的稳定性进行具体研究. 若有常数 a 是差分方程(4.2.1)的解,即 F ( n , a , a , ? , a ) ? 0 ,则称 a 是差分方程 (4.2.1)的平衡点.若对差分方程(4.2.1)由任意初始条件确定的解 xn 都有xn ? a (n ? ? )则称这个平衡点 a 是稳定的. 一、一阶线性差分方程的平衡点和稳定性 一阶线性常系数差分方程x n ?1 ? a x n ? b(4.2.6) .的平衡点可以由 x ? ax ? b 解得,为 x * ?b 1? a注意到,通过变量代换可以把方程(4.2.6)的平衡点的稳定性问题转换为x n ?1 ? a x n ? 0(4.2.7)的平衡点 x * ? 0 的稳定性问题.根据上一节介绍的齐次线性差分方程的解法,可知, 方程(4.2.7)的解可以表示为xn ? ( ? a ) x0n根据稳定性的定义,可知,当且仅当 a ? 1 时,方程(4.2.7)的平衡点才是稳定的, 从而方程(4.2.6)的平衡点也是稳定的. 对于 m 维向量 x ( n ) 和 m ? m 常数矩阵 A 构成的方程组x ( n ? 1) ? A x ( n ) ? 0(4.2.8)其平衡点稳定的条件是 A 的特征值 ? i , i ? 1, 2, ? ,均有?i ? 1即均在复平面上的单位圆内. 二、二阶线性差分方程的平衡点和稳定性 二阶齐次线性差分方程x n ? 2 ? a 1x n ? 1 a x2n ? 0 ?(4.2.9)的平衡点为 x * ? 0 .差分方程(4.2.9)的特征方程为 ? 2 ? a1 ? ? a 2 ? 0 ,记特征根为? 1 , ? 2 ,则方程(4.2.9)的通解可以表示为x n ? C 1? 1 ? C ? 2n n 2(4.2.10)其中常数 C 1 , C 2 可以由初始条件 x 0 , x1 确定.由(4.2.10)可得,当且仅当?i ? 1( i ? 1, 2 )时方程(4.2.9)的平衡点才是稳定的. 与一阶线性方程类似,通过变量代换可得,非齐次线性方程x n ? 2 ? a 1x n ? 1 a x2n ? b ?(4.2.11) 的平衡点的稳定性和方程(4.2.9)相同. 二阶方程的上述结果可以推广到 k 阶线性方程,即 k 阶方程x n ? k ? a1 x n ? k ?1 ? ? ? a k x n ? b的平衡点稳定的充要条件是其对应的特征方程的特征根 ? i ,均有?i ? 1i ?( 1, 2 , ? n. , )(4.2.12)注 2：高阶方程与一阶方程组通过一定的变量代换是可以相互转化的,我们注意 到,(4.2.12)这一条件与一阶方程组(4.2.8)的平衡点稳定性条件是一致的. 三、一阶非线性差分方程的平衡点和稳定性 一阶非线性差分方程x n ?1 ? f ( x n )(4.2.13)的平衡点 x * 可以由代数方程 x ? f ( x ) 解出.为分析 x * 的稳定性,我们采用与微分 方程类似的线性近似的方法,将方程(4.2.13)的右端在 x * 处作 Taylor 展开,只取到 一次项,得到方程(4.2.13)的线性近似方程x n ? 1 ? f ( x ) ? f '( x )( x n ? x )* * *(4.2.14)显然, x * 也是(4.2.14)的平衡点.线性方程(4.2.14)的平衡点的讨论与(4.2.6)相 同,而当 f '( x * ) ? 1 时,方程(4.2.14)与(4.2.13)的平衡点的稳定性相同.于是我们得到 ①当 f '( x * ) ? 1 时,方程(4.2.13)的平衡点 x * 是稳定的; ②当 f '( x * ) ? 1 时,方程(4.2.13)的平衡点 x * 是不稳定的.§4.2.3 市场经济的蛛网模型在自由贸易的集市上你注意过这样的现象吗？一个时期由于猪肉的上市量远大于需求, 销售不畅导致价格下降, 农民觉得养猪赔钱, 于是转而经营其它农副业.过一段时间后猪肉 上市量大减, 供不应求导致价格上涨, 原来的饲养户看到有利可图, 又重操旧业.这样下一个 时期会重现供大于求、价格下降的局面, 在没有外界干扰的情况下, 这种现象将如此循环下 去. 在完全自由竞争市场经济中上述现象通常是不可避免的, 因为商品的价格是由消费者 的需求关系决定的, 商品数量越多价格越低.而下一时段商品的数量由生产者的供应关系决 定, 商品价格越低生产的数量就越小.这样的需求和供应关系决定了市场经济中商品的价格 和数量必然是振荡的.这种振荡越小越好, 如果振荡太大就会影响人民群众的正常生活. 现在我们用差分方程建模, 讨论市场经济趋于稳定的条件, 再用图形方法建立“蛛网模 型”对上述现象进行分析, 对结果进行解释, 然后作适当推广. 记第 n 时段商品数量为 xn, 价格为 yn, n = 1, 2, ?. 这里我们把时间离散化为时段, 1 个 时段相当于商品的 1 个生产周期, 如蔬菜、水果可以是一年, 肉类则是 1 个饲养周期. 在 n 时段商品的价格 yn 取决于数量 xn, 设 yn = f (xn ). 它反映消费者对这种商品的需求关 系, 称需求函数.因为商品的数量越多价格越低, 所以在图 4-8 中用一条下降曲线 f 表示它, f 为需求曲线. 在 n + 1 时段商品的数量 xn+1 由上一时段价格 yn 决定, 用 xn+1 = g ( yn )表示, 它反映生产者的供应关系, 称供应函 数.因为价格越高生产量才越大, 所以在图 4-8 中供应曲线 g 是一条上升曲线. 设图 4-8 中两条曲线 f 和 g 相交于 P0(x0, y0)点, 在 P0 点附近取函数 f 和 g 的线性近似, 即 需求函数( f )：yn ? y0 = ????(xn ??x0), ??＞0 4-13 图 4 -8 供应函数(g)：xn+1 ??x0 = ???( yn ?? y0 ), ??＞0. 4-14 由式(4-13), (4-14)消去 yn 得一阶线性差分方程 xn+1 = ??????xn + (1 + ????) x0, n = 1,2,?. 4-5 因此 x0 是其平衡点, 即 P0 是平衡点, 对 n 递推不难得到 xn+1 = (???????) n x1 +[1? (??????) n ]x0, n = 1,2,?. 由此可得, 平衡点 P0 稳定的条件是：??? ＜1; 不稳定的条件是：??? ＞1. 下面用图形解释此模型. 若对某一个 k 有 xk = x0, 则由式(4-5)递推可得, 当 n≥k 时 xn = x0, 从而 yn = y0, 即商品的 数量和价格将永远保持在 P0(x0, y0)点.但是实际生活中的种种干扰使得 xn, yn 不可能停止在 P0 点上.不妨设 x1 偏离 x0 (如图 4-9, 图 4-10), 我们分析随着 n 的增加 xn, yn 的变化.图 4 -9 P0 点是稳定的图 4 -10 P0 点是不稳定的数量 x1 给定后, 价格 y1 由曲线 f 上的 P1 点决定, 下一时段的数量 x2 由曲线 g 上的 P2 点 决定, 这样得到一系列的点 P1(x1, y1), P2(x2, y1), P3(x2, y2), P4(x3, y2),?,在图 4-9 上这些点将接 箭头所示方向趋向 P0(x0, y0), 表明 P0 是稳定平衡点, 意味着市场经济（商品的数量和价格） 将趋向稳定.但是如果需求函数和供应函数由图 4-10 的曲线所示, 则类似的分析发现, 市场 经济将按照 P1, P2, P3, P4,?, 的规律变化而远离 P0, 即 P0 是不稳定平衡点, 市场经济趋向不 稳定. 由此可见, 需求曲线越平, 供应曲线越陡, 越有利于经济稳定. 图 4-9 和图 4-10 中折线 P1P2P3P4?形似蛛网, 于是这种用需求曲线和供应曲线分析市场 经济稳定性的图示法在经济学中称蛛网模型.实际上, 需求曲线 f 和供应曲线 g 的具体形式通 常是根据各个时段商品的数量和价格的一系列统计资料得到的.一般地说, f 取决于消费者对 这种商品的需要程度和他们的消费水平, g 则与生产者的生产能力, 经营水平等因素有关. 下面再来解释此模型的实际意义. 首先考察参数?, ??的含义, 需求函数 f 的斜率??(取绝对值)表示商品供应量减少 1 个单位 时价格的上涨幅度；供应函数 g 的斜率?表示价格上涨 1 个单位时(下一时期)商品供应增加 量.所以? 的数值反映消费者对商品需求的敏感程度. 如果这种商品是生活必需品, 消费者 处于持币待购状态, 商品数量稍缺, 人们立即蜂拥抢购, 那么??会比较大； 反之, 若这种商品 非必需品, 消费者购物心理稳定, 或者消费水平低下, 则??较小. ? 的数值反映生产经营者 的对商品价格的敏感程度, 如果他们目光短浅, 热衷于追逐一时的高利润, 价格稍有上涨立 即大量增加生产, 那么??会比较大；反之, 若他们素质较高, 有长远的计划, 则???较小. 根据?, ??的意义很容易对市场经济稳定与否的条件作出解释.当供应函数 g, 即???固定 时, ? 越小, 需求曲线越平, 表明消费者对商品需求的敏感程度越小, 越利于经济稳定.当需 求函数 f, 即? 固定时, ???越小, 供应曲线越陡, 表明生产者对价格的敏感程度越小, 越利于 经济稳定.反之, 当??? 较大, 表明消费者对商品的需求和生产者对商品的价格都很敏感, 则 会导致经济不稳定. 当市场经济趋向不稳定时, 政府有两种干预办法：一种办法是控制价格, 无论商品数量 多少, 命令价格不得改变, 于是??= 0, 不管曲线 g 如何, 总是稳定的；另一种办法是控制市 场上的商品数量, 当上市量小于需求时, 政府从外地收购或调拔, 投入市场, 当上市量多于 需求时, 政府收购过剩部分, 于是??= 0, 不管曲线 f 如何, 也总是稳定的. 最后我们将模型稍加推广. 如果生产者的管理水平更高一些, 他们在决定商品生产数量时, 不是仅根据前一时期的 价格, 而是根据前两个时期的价格, 为简单起见不妨设根据二者的平均值(yn + yn?1)/2, 于是 供应函数为 xn+1 = g[(yn + yn?1)/2].在 P0 点附近取线性近似时 供应函数(g)：xn+1 ??x0 = ? [(yn + yn?1)/2? y0], ??＞0. 4-16 又设需求函数仍由(1)式表示, 则由(4-13),(4-16)式得到 2xn+2 + ?? xn+1 + ?? xn = (1 + ?? )x0, n =1,2,?. 4-17 (4-17)式是二阶线性差分方程. P0 点稳定的条件可由其特征方程 2??2 +???? + ??? = 0 的根??1, 2 =? ?? ? (?? ) 2 ? 8?? 4确定.当?? ＞8 时, 显然有??2 =? ?? ?(?? ) 2 ? 8?? 4? ???4＜?2,从而|???2 |＞2, P0 是不稳定的. 当???＜8 时, 可以算出|???1, 2 | =??2, 由|???1, 2 |＜1 得到 P0 点稳定的条件为??? ＜2.这与原有模型中 P0 点稳定的条件??? ＜1 相比, 保持经济稳定的参数??,????的范围放大了(?, ???的含义未变). 可以想到, 这是生产经营者的生产管理水平提高, 对市场经济稳定起着有利 影响的必然结果.参考文献[1]胡良剑,孙晓军.Matlab 数学实验.北京:高等教育出版社,2006. [2]姜启源,谢金星,叶俊.数学模型(第三版).北京:高等教育出版社,2003. [3]刘承平.数学建模方法.北京:高等教育出版社,2002. [4]薛定宇,陈阳泉.高等应用数学问题的 Matlab 求解.北京:清华大学出版社,2006 [5]王高雄,周之铭,朱思铭,王寿松.常微分方程(第二版).北京:高等教育出版社,1983.习题 41. 一角度为 60 度圆锥形漏斗里装着 10cm 高的水, 其下端小孔的面积为 0.6cm2,求这此水流完需用多少 时间？ 2. 某容器温度为 60℃, 将其内的温度计移入另一容器中, 十分钟后读数为 70℃, 又过十分钟后读数为 76℃, 问另一容器的温度为多少度？（提示：温度变化率与温度差成正比） 3. 一大学教师小李从 31 岁开始建立自己的养老基金, 他把已有的积蓄 1 万元也一次性地存入, 已知月 利率为 0.01（以复利计）, 每月存入 300 元, 试问当小李 60 岁退休时, 他的退休基金有多少？又若, 他退休 后每月要从银行提取 1000 元, 试问多少年后他的退休基金将用完？你能否根据你了解的实际情况建立一个 较好的养老基金的数学模型及相应的算法和程序(软件). 4. 某人每天由饭食获取 2500 卡热量, 其中 1200 卡用于新陈代谢, 此外每公斤体重需支付 16 卡热量作 为运动消耗, 其余热量则转化为脂肪, 已知以脂肪形成贮存的热量利用率为 100%, 每公斤脂肪含热量 10000 卡, 问此人的体重如何随时间而变化？ 5. 1650 年世界人口为 5 亿, 当时的年增长率为 3%.用指数增长模型计算什么时候世界人口达到 10 亿. (实际上 1850 年前已超过 10 亿), 1970 年世界人口为 36 亿, 年增长率为 2.1%, 用指数增长模型预测什么时 候世界人口会翻一番(这个结果可信吗？)你对同样的模型得到的两个结果有什么看法. 6. 设某城市共有 n + 1 人, 其中一人出于某种目的编造了一个谣言.该城市具有初中以上文化程度的人 占总人数一半, 这些人只有 1/4 相信这一谣言, 而其他人约有 1/3 会相信.又设凡相信此谣言的人在每单位时 间内传播的平均人数正比于时当未听说此谣言的人数, 而不相信此谣言的人不传播谣言.试建立一个反映谣 言传播情况的微分方程模型. 7. 如果两个种群都能独立生存, 共处时又能相互提供食物, 试建立种群相互依存模型并讨论平衡点的 稳定性, 解释稳定的意义.又如果两个种群都不能独立生存, 试建模以讨论共处的可能性. 8. 请你捕食--被捕食模型加以改进, 并讨论平衡点的稳定性, 解释稳定的意义. 9. 如果在捕食--被捕食模型中, 捕食者掠食的对象只是成年的被捕食者, 而未成年的因体积太小免遭 捕获, 在适当的假设下建立这三者之间关系的模型, 求平衡点. 10. 大陆上物种数目是常数, 各物种独立地从大陆向附近一岛迁移, 岛上物种数量的增加与尚未迁移 的物种数目有关, 而随着迁移物种数的增加又导致岛上物种的减少, 在适当假设下建立岛上物种数的模型, 并讨论稳定状况. 11. 对于技术革新的推广, 在下列几种情况下分别建立模型. ① 推广工作通过已经采用新技术的人进行, 推广速度与已采用新技术的人数成正比, 推广是无限的. ② 总人数有限, 因而推广速度随着尚未采用新技术人数的减少而降底. ③ 在②的前提下还要考虑广告等媒介的传播作用. 12. 在市场经济的蛛网模型中, 继续讨论下列问题： ① 因为一个时段上市的商品不能立即售完, 其数量也会影响到下一段的价格, 所以 yn+1 由 xn+1 和 xn 决 定, 再设 xn+1 仍只取决于 yn, 给出稳定平衡的条件. ② 若除了 yn+1 由 xn+1 和 xn 决定之外, xn+1 也由 yn 和 yn?1 确定, 试分析稳定平衡的条件是否还会放宽. 13. 在按年龄分组的的种群增长模型中, 设一群的动物最高年龄为 15 岁, 每 5 岁一组分成 3 个年龄组, 各组的繁殖率为 b1 = 0, b2 = 4, b3 = 3, 存活率为 s1 = 1/2, s2 = 1/4 , 开始时 3 组各有 1000 只.求 15 年后各组分别 有多少只？以及时间充分长以后种群的增长率和按年龄组的分布. 14. 最优捕鱼策略(1996 年全国大学生数学建模竞赛题 A 题)为了保护人类赖以生存的自然环境, 可再 生资源(如渔业, 林业资源)的开发必须适度, 一种合理, 简化的策略是, 在实际可持续收获的前提下, 追求 最大产量或最佳效益. 考虑对某种鱼(鱼)的最优捕捞策略：假设这种鱼分 4 个年龄组, 称 1 龄鱼, …, 4 龄鱼, 各年龄组每条 鱼的平均重量分别为 5.07, 11.55, 17.86, 22.99(g), 各年龄组自然死亡率均为 0.8(1/年), 这种鱼为季节性集中 产卵繁殖, 平均每条 4 龄鱼的产卵量为 1.109×105(个), 3 龄鱼的卵量为这个数的一半, 2 龄鱼和 1 龄鱼不产 卵, 产卵和孵卵期为每年的最后 4 个月, 卵孵化并成活为 1 龄鱼, 成活率(1 龄鱼条数与产卵总量 n 之比)为 1.22××1011+n). 渔业管理部门, 每年只允许在产卵孵化期前的 8 个月内进行捕捞作业, 如果每年投入的捕捞能力(如渔 船数, 下网次数等)固定不变, 这时单位时间捕捞量将与各年龄组鱼数成正比, 比例系数不妨称为捕捞强度 系数, 通常使用 13mm 网眼的拉网, 这种网只能捕捞 3 龄鱼和 4 龄鱼, 其中两个捕捞强度系数之比为 0.42:1. 渔业上称为这种方式为固定努力量捕捞. ① 建立数学模型分析如何实现可持续捕捞(即每年开始捕捞渔场中各年龄组鱼群条不变)并且在此前 提下得到最高的年收获量(捕捞总重量). ② 某渔业公司承包这种鱼的捕捞业务 5 年, 合同要求 5 年后鱼群的生产能力不受到太大破坏, 已知承 包时各年龄组鱼群的数量分别为：122, 29.7, 10.1, 3.29(×109 条), 如果要用固定努力量的捕捞方式, 该公司 应采取怎样的策略才能使总收获量最高.
更多相关文档}
我要回帖
更多关于 一阶齐次线性微分方程 的文章
更多推荐
- ·河南省人力资源和社会保障厅河南社保卡是不是医保卡卡吗
- ·微纳金属3d打印金属打印技术应用：AFM探针？
- ·三门峡高考状元前100，全省能排多少？
- ·有没有大神能帮我翻译图片文字出这张图片的内容？
- ·小太阳鹦鹉最新法规可以合法开放了吗？
- ·珠海地坪平整度测试仪漆不平整应该怎么处理？
- ·有关复合函数的微分法微分的问题
- ·美国有一半总统有爱尔兰裔美国总统血统吗？
- ·计量泵如何计选择合适的脉冲阻尼器 选型容积？
- ·囊括四海彻仙路云霄 全文 手打中打一动物
- ·聚氨脂环氧砂浆地坪坪材料买哪家的呢
- ·哪位大哥帮忙用名字，夏彩蝶舞夏 何真真为名写首诗，以表达我对她的爱慕之情。
- ·No I cannot be overloadedI want to be onWith me on the
- ·您好，想moto请教您一个问题题
- ·想学铁路造价, 看不懂铁路施工图预算与工程造价,应该怎么入手,先看什么再看什么
- ·有6张信用卡逾期每月还一点都逾期了2个月，能在平安银行工作吗？
- ·SCALE，ASCALE，这两个步指令是什么意思思
- ·二阶一阶非线性微分方程程可不可以没有一阶导数?！！！！很迷！
- ·海苔丝可以高温黄铜烘烤后会变色吗吗?
- ·阜新哪家驻车燃油加热器那个好好？
- ·招商蛇口研究生待遇是国企吗
- ·要是请个长期的物业钟点工申请表谁家比较好啊？在网上看见睿智宏达家政的口碑还可以吗？他家行吗
- ·今天是第二次考科目二考试紧张怎么办三了，好紧张，这几天模拟都没有一次过的，肿么办，真不想考的了😭&
- ·我体内湿气重怎么调理喝*这个可以明显改善吗
- ·河源市看抑郁症去广州协佳耳鼻喉医院医院能看吗
- ·天津二复属医院能滴尿检测数据复核核基制蛋白吗？
- ·前列腺炎怎么引起的有什么症状
- ·岳阳市世纪妇科世纪医院做人流确找老牌专业妇科医院要做人流的话应该怎么办？
- ·肾虚耳鸣可以吃枸杞会引起耳鸣吗？
- ·这种心绞痛症状是什么怎么回事啊？
- ·请问扬州哪家助听器好？老人助听器需要买一个
- ·请问中药怎么治疗多囊卵巢巢用中药调理还是用西药哪个更好
- ·复清发展通丁是什么药
