matlab教学第二章微分方程

1、第二章第二章 微分方程模型微分方程模型2021-12-131本章学习目的：本章学习目的：l学习微分方程模型的建立、求解方法、分析结学习微分方程模型的建立、求解方法、分析结果及解决实际问题的全过程。果及解决实际问题的全过程。l熟练掌握使用熟练掌握使用MATLAB软件的函数求微分方软件的函数求微分方程的解析解、数值解和图形解程的解析解、数值解和图形解 。2021-12-1322.1 引例引例 l对于圆柱形状容器壁上的容积刻度，对于圆柱形状容器壁上的容积刻度，可以利用圆柱体体积公可以利用圆柱体体积公式：式： ，其中容器的直径，其中容器的直径D为常数，体积为常数，体积V与相对于容器底部的与相对于容器底

2、部的任意高度任意高度H成正比，因此在容器壁上成正比，因此在容器壁上可以方便地标出容积刻度。可以方便地标出容积刻度。l而对于几何形状不规则的容器，比如而对于几何形状不规则的容器，比如“倒葫芦形状倒葫芦形状”的容器壁上如何标出的容器壁上如何标出容积刻度呢？容积刻度呢？x4/2HDV2021-12-133l建立坐标系，由微元法分析可知：建立坐标系，由微元法分析可知：l l其中其中x表示高度，直径是高度的函数，记为表示高度，直径是高度的函数，记为D(x)l可得微分方程：可得微分方程：l如果该方程中的函数如果该方程中的函数D(x)无解析表达式，只给出无解析表达式，只给出D(x)的部分测试数据，如何求解此

3、微分方程呢？的部分测试数据，如何求解此微分方程呢？dxxDdV2)(410)0()(412VxDdxdV2021-12-134lh=0.2;ld=0.04,0.11,0.26,0.56,1.04,1.17;lx(1)=0;v(1)=0;lfor k=1:5lx(k+1)=x(k)+h;lv(k+1)=v(k)+(h/2)*(pi/4)*(d(k)2+d(k+1)2);lendlx=x(1:6),v=v(1:6),lplot(x,v) lx =l Columns 1 through 5 l 0 0.2000 0.4000 0.6000 0.8000l Column 6 l 1.0000lv =l

4、 Columns 1 through 5 l 0 0.0011 0.0073 0.0373 0.1469l Column 6 l 0.339300.10.20.30.40.50.60.70.80.9100.050.10.150.20.250.30.352021-12-1352.2 微分方程模型的建立微分方程模型的建立 l在工程实际问题中，在工程实际问题中，“改变改变”、“变化变化”、“增加增加”、“减少减少”等关键词提示我们注意什么量在变化，关键词等关键词提示我们注意什么量在变化，关键词“速率速率”、“增长增长”、“衰变衰变”、“边际的边际的”等常涉及到等常涉及到导数。导数。l我们熟悉的速度公

5、式：我们熟悉的速度公式： 就是一个简单的一阶微分方就是一个简单的一阶微分方程。程。l微分方程是指含有导数或微分的等式。微分方程是指含有导数或微分的等式。l一般形式：一般形式：l常用的建立微分方程的方法有：运用已知物理定律；利用常用的建立微分方程的方法有：运用已知物理定律；利用平衡与增长式；运用微元法；应用分析法。平衡与增长式；运用微元法；应用分析法。vdtdy).,(0),()1()()(nnnyyyxfyyyyxF或：2021-12-1362.2.1 运用已知物理定律运用已知物理定律l例例2.1 一个较热的物体置于室温为一个较热的物体置于室温为180C的房间内，的房间内，该物体最初的温度是该

6、物体最初的温度是600C，3分钟以后降到分钟以后降到500C。想。想知道它的温度降到知道它的温度降到300C 需要多少时间？需要多少时间？10分钟以后分钟以后它的温度是多少？它的温度是多少？l牛顿冷却（加热）定律牛顿冷却（加热）定律：将温度为：将温度为T的物体放入处于的物体放入处于常温常温 m 的介质中时，的介质中时，T的变化速率正比于的变化速率正比于T与周围介与周围介质的温度差。质的温度差。l分析：假设房间足够大，放入温度较低或较高的物体分析：假设房间足够大，放入温度较低或较高的物体时，室内温度基本不受影响，即室温分布均衡，保持时，室内温度基本不受影响，即室温分布均衡，保持为为m，采用牛顿冷

7、却定律是一个相当好的近似。，采用牛顿冷却定律是一个相当好的近似。2021-12-137l建立模型：设物体在冷却过程中的温度为建立模型：设物体在冷却过程中的温度为T(t)，t0， l根据牛顿加热（冷却）定律：根据牛顿加热（冷却）定律： ，建立微分方程建立微分方程 l其中参数其中参数k 0，m=18。）成正比与（mTdtdT60)0()(TmTkdtdT2021-12-1382.2.2 利用平衡与增长式利用平衡与增长式 l许多研究对象在数量上常常表现出某种不变的特性，许多研究对象在数量上常常表现出某种不变的特性，如封闭区域内的能量、货币量等。利用变量间的平衡如封闭区域内的能量、货币量等。利用变量间

8、的平衡与增长特性，可分析和建立有关变量间的相互关系。与增长特性，可分析和建立有关变量间的相互关系。l此类建模方法的关键是分析并正确描述基本模型的右此类建模方法的关键是分析并正确描述基本模型的右端，使平衡式成立。端，使平衡式成立。l例例2.2 战斗模型：两方军队交战，希望为这场战斗建战斗模型：两方军队交战，希望为这场战斗建立一个数学模型，应用这个模型达到如下目的：立一个数学模型，应用这个模型达到如下目的： l1. 预测哪一方将获胜？预测哪一方将获胜？ l2. 估计获胜的一方最后剩下多少士兵？估计获胜的一方最后剩下多少士兵？ l3. 计算失败的一方开始时必须投入多少士兵才能赢得计算失败的一方开始时

9、必须投入多少士兵才能赢得这场战斗？这场战斗？ 2021-12-139l解：模型建立：解：模型建立：l设设 x(t)： t 时刻时刻X方存活的士兵数方存活的士兵数l y(t)： t 时刻时刻Y方存活的士兵数方存活的士兵数l假设：假设：l1）双方所有士兵不是战死就是活着参加战斗，）双方所有士兵不是战死就是活着参加战斗，x(t)与与y(t)都是连续变量；都是连续变量； l2）Y方军队的一个士兵在单位时间内杀死方军队的一个士兵在单位时间内杀死X 方军队方军队 a 名名士兵；士兵； l3）X 方军队的一个士兵在单位时间内杀死方军队的一个士兵在单位时间内杀死Y方军队方军队 b 名名士兵；士兵；lt 时间内

10、时间内X军队减少的士兵数军队减少的士兵数 = t 时间内时间内Y军队消灭军队消灭对方的士兵数对方的士兵数 l即有即有 x =aytl同理同理 y =bxt 令0t， 得得 )0()0(bbxdtdyaaydtdx2021-12-13102.2.3 微元法微元法 l基本思想：通过分析研究对象的有关变量在一个基本思想：通过分析研究对象的有关变量在一个很短时间内的变化情况建立微分方程。很短时间内的变化情况建立微分方程。l本章引例本章引例2021-12-13112.2.4 分析法分析法l基本思想：根据对现实对象特性的认识，分析其因果基本思想：根据对现实对象特性的认识，分析其因果关系，找出反映内部机理的

11、规律。关系，找出反映内部机理的规律。l例例2.4 独家广告模型独家广告模型 广告是调整商品销售的强有力广告是调整商品销售的强有力的手段，广告与销售量之间有什么内在联系？如何评的手段，广告与销售量之间有什么内在联系？如何评价不同时期的广告效果？价不同时期的广告效果？l解：解：1分析广告的效果，可做如下的条件假设：分析广告的效果，可做如下的条件假设：l商品的销售速度会因广告而增大，当商品在市场上商品的销售速度会因广告而增大，当商品在市场上趋于饱和时，销售速度将趋于一个极限值；趋于饱和时，销售速度将趋于一个极限值；l商品销售率（销售加速度）随商品销售速度的增高商品销售率（销售加速度）随商品销售速度的

12、增高而降低。而降低。 2021-12-1312l2符号说明符号说明lA(t) t 时刻的广告费用时刻的广告费用lS(t) t 时刻商品的销售速度；时刻商品的销售速度； lM 销售饱和水平，即销售速度的上限；销售饱和水平，即销售速度的上限； l 衰减因子，广告作用随时间的推移而自然衰减因子，广告作用随时间的推移而自然衰减的速度，衰减的速度，0；lp 响应系数，表征响应系数，表征A(t) 对对 S(t) 的影响力。的影响力。2021-12-1313l3模型建立模型建立l选择如下广告策略，选择如下广告策略，t时刻的广告费用为：时刻的广告费用为：l建立微分方程：建立微分方程：l模型分析：是否与假设相符

13、？模型分析：是否与假设相符？ttAtA, 00,)()()(1)(tSMtStpAdtdS2021-12-13142.3 微分方程求解方法微分方程求解方法 l2.3.1 解析解法解析解法l解析解法只能解决一些特殊微分方程，这些方解析解法只能解决一些特殊微分方程，这些方法主要针对：法主要针对：l一阶特殊的微分方程：如使用分离变量法、方一阶特殊的微分方程：如使用分离变量法、方程变换法、线性方程的常数变易法或公式法求程变换法、线性方程的常数变易法或公式法求解。解。l二阶或高阶常系数线性微分方程的特征根法。二阶或高阶常系数线性微分方程的特征根法。l在高等数学的教程中有专门介绍。在高等数学的教程中有专门

14、介绍。l下面着重介绍微分方程的数值解法。下面着重介绍微分方程的数值解法。2021-12-13152.3.2 数值解法数值解法l微分方程的数值解法是解决某些实际问题中经常使用的方微分方程的数值解法是解决某些实际问题中经常使用的方法。设待求解的定解问题为法。设待求解的定解问题为l求该问题数值解法的基本过程如下：求该问题数值解法的基本过程如下：l引入自变量取值点序列引入自变量取值点序列 ，定义，定义 为步长，常为步长，常用定步长用定步长( 与与n无关，为常数无关，为常数)，其精确解记为，其精确解记为 ，一般难以得到。为了寻求一般难以得到。为了寻求 的近似值的近似值 ，设想根据一，设想根据一定的原理，

15、结合当前得到近似解，近似地表示该点或前一定的原理，结合当前得到近似解，近似地表示该点或前一点的导数值，由此推出计算点的导数值，由此推出计算 的迭代公式。因此数值解法的迭代公式。因此数值解法一般只能得到微分方程的近似解一般只能得到微分方程的近似解 。下面介绍两个微分。下面介绍两个微分方程中最常用的数值解法。方程中最常用的数值解法。00)(),(yxyyxfdxdy nx1nnnxxhnh)(nxy)(nxynyny ny2021-12-13161.欧拉方法欧拉方法l这是一种最简单的解微分方程的数值方法：就是在小区间这是一种最简单的解微分方程的数值方法：就是在小区间xn, xn+1上用差商代替微商

16、，可以得到近似的表达式上用差商代替微商，可以得到近似的表达式l若若f(x,y)中的中的x取左端点取左端点 ，结合已经得到的，结合已经得到的y(xn)的近似的近似值（数值解）值（数值解）yn，即，即 ，有，有y(xn+1)的近似值为的近似值为l n = 0,1, l这就是求解微分方程的这就是求解微分方程的显式欧拉公式显式欧拉公式。也称。也称向前欧拉公式向前欧拉公式。l向前欧拉法计算简单，易于计算，但精度不高，收敛速度向前欧拉法计算简单，易于计算，但精度不高，收敛速度慢慢),()()(1yxfhxyxynnnnyxy)(),(1nnnnyxfhyynx2021-12-1317l若若f(x,y)中的

17、中的x取右端点，可得取右端点，可得向后欧拉公式向后欧拉公式如下：如下：lyn+1 = yn + h f (xn +1, yn +1) ，n = 0,1,l称为隐式公式，因为要得出数值解称为隐式公式，因为要得出数值解yn+1，就必须求解这个，就必须求解这个非线性方程，计算比较困难。非线性方程，计算比较困难。l如果用将向前和向后欧拉公式加以平均，可得到如果用将向前和向后欧拉公式加以平均，可得到梯形公式梯形公式：l该法的计算精度比向前和向后欧拉法都高，但计算和向后该法的计算精度比向前和向后欧拉法都高，但计算和向后欧拉法一样困难。欧拉法一样困难。),(),(21111nnnnnnyxfyxfhyy20

18、21-12-1318l改进的欧拉算法：改进的欧拉算法：l（1）先用向前欧拉法算出）先用向前欧拉法算出yn+1的预测值的预测值 ，ll（2）将预测值代入梯形公式的右端作为校正，得到）将预测值代入梯形公式的右端作为校正，得到yn+1l ，n=1,2,l该式称为该式称为改进欧拉公式改进欧拉公式。1ny),(1nnnnyxfhyy),(),(21111nnnnnnyxfyxfhyy2021-12-1319l例例2.5 求解微分方程求解微分方程ly = -y +x +1, y(0) = 1, 取步长取步长h = 0.1和和0.001。l分别用三种数值解法求解，并结合其精确解，对求解误差进分别用三种数值解

19、法求解，并结合其精确解，对求解误差进行分析比较。行分析比较。l解解 这是一个一阶线性微分方程，可用解析解法得到其精确解这是一个一阶线性微分方程，可用解析解法得到其精确解y = x + e-x。2021-12-1320l三种数值解如下：三种数值解如下：l 1) 向前欧拉法：迭代公式为向前欧拉法：迭代公式为 yn+1 = （1-h）yn + hxn + h，n=0,1,.。其中。其中y0= y(0) = 1。l2) 后退欧拉法：由后退欧拉法隐式公式得后退欧拉法：由后退欧拉法隐式公式得lyn+1 = yn + h(-yn+1+xn+1+1)，变形为，变形为lyn+1 = (yn + hxn+1 +

20、h)/(1+h)。l3) 梯形法：将隐式梯形公式转化为显示迭代公式梯形法：将隐式梯形公式转化为显示迭代公式如下：如下：lyn+1 = (yn + (h/2)*(- yn + xn + xn+1 +2)/(1+h/2)。2021-12-1321lx1(1)=0;y1(1)=1;y2(1)=1;y3(1)=1;h=0.1;lfor k=1:10lx1(k+1)=x1(k)+h;ly1(k+1)=(1-h)*y1(k)+h*x1(k)+h;ly2(k+1)=(y2(k)+h*x1(k+1)+h)/(1+h);ly3(k+1)=(y3(k)+(h/2)*(-y3(k)+x1(k)+x1(k+1)+2)

21、/(1+h/2);lendlx=0:0.1:1;ly=x+exp(-x);lx1=x1(1:11),y=y(1:11),y1=y1(1:11),y2=y2(1:11),y3=y3(1:11),lplot(x,y,x1,y1,k:,x1,y2,r-,x1,y3,g*) l程序中，程序中，x1为自变量，为自变量，y为精确解，为精确解，y1、y2、y3分别为向分别为向前欧拉法、后退欧拉法和梯形法的解。结果如下前欧拉法、后退欧拉法和梯形法的解。结果如下：2021-12-1322表表2-1当当h = 0.1时时nx精确解精确解向前欧拉法向前欧拉法后退欧拉法后退欧拉法梯形法梯形法011110.11.004

22、811.00911.00480.21.01871.01001.02641.01860.31.04081.02901.05131.04060.41.07031.05611.08301.07010.51.10651.09051.12091.10630.61.14881.13141.16451.14850.71.19661.17831.21321.19630.81.24931.23051.26651.24900.91.30661.28741.32411.306311.36791.34871.38551.36762021-12-132300.10.20.30.40.50.60.70.80.9111.0

23、51.11.151.21.251.31.351.4图中，蓝色曲线是精确解，黑色曲线是向前欧拉法曲线，图中，蓝色曲线是精确解，黑色曲线是向前欧拉法曲线，红色曲线是向后欧拉法曲线，绿色红色曲线是向后欧拉法曲线，绿色“* ”号为梯形法曲线。号为梯形法曲线。24表表2-2 当当h = 0.001时时 nx精确解精确解向前欧拉法向前欧拉法后退欧拉法后退欧拉法梯形法梯形法011110.11.00481.00481.00491.00480.21.01871.01861.01881.01870.31.04081.04071.04091.04080.41.07031.07021.07051.07030.51.1

24、0651.10641.10671.10650.61.14881.14861.14901.14880.71.19661.19641.19681.19660.81.24931.24911.24951.24930.91.30661.30641.30681.306611.36791.36771.36811.36792021-12-132500.20.40.60.811.21.411.051.11.151.21.251.31.351.4 计算结果表明：当步长计算结果表明：当步长h =0.1时，它们的前两位有效数字时，它们的前两位有效数字是精确的；当步长是精确的；当步长h =0.001时，它们的前四位有效

25、数字是精确时，它们的前四位有效数字是精确的。说明在迭代中，步长的。说明在迭代中，步长h越小，计算结果越精确。越小，计算结果越精确。 2021-12-1326l衡量求解公式好坏的一个主要标准是求解公式的精度，因衡量求解公式好坏的一个主要标准是求解公式的精度，因此引入局部此引入局部截断误差和阶数截断误差和阶数的概念。的概念。l定义定义2.1 假定假定 没有误差，即没有误差，即 ，用数值方法计，用数值方法计算算 的误差的误差 ，称为该数值方法计算时的，称为该数值方法计算时的局部截断误差。局部截断误差。l为估计欧拉公式的局部截断误差，先将精确解为估计欧拉公式的局部截断误差，先将精确解 在在 处作泰勒级

26、数展开处作泰勒级数展开l对于向前欧拉公式，在对于向前欧拉公式，在 的假定下可记作的假定下可记作l式中式中h的最低阶项称为局部截断误差主项，它对于的最低阶项称为局部截断误差主项，它对于h是是2阶阶的。的。ny)(nnxyy 1ny11)(nnyxy)(1nxynx)()(! 2)()()(321hOxyhxyhxyxynnnn )()()(,()(1nnnnnnxyhxyxyxhfxyy)(nnxyy )()()(2)(232111hOhOxyhyxynnnn 2021-12-1327l同理，对于向后欧拉公式的局部截断误差同理，对于向后欧拉公式的局部截断误差是是 ，梯形公式和改进欧拉公式的局部，

27、梯形公式和改进欧拉公式的局部截断误差是截断误差是 。l求解微分方程的数值计算方法的精度是由局部截求解微分方程的数值计算方法的精度是由局部截断误差中断误差中h的阶定义的：如果一个方法的局部截的阶定义的：如果一个方法的局部截断误差为断误差为 ，则称该方法具有，则称该方法具有p阶精度。阶精度。l因此向前和向后欧拉方法的精度为因此向前和向后欧拉方法的精度为1阶，梯形公阶，梯形公式和改进欧拉公式的精度为式和改进欧拉公式的精度为2阶。阶。)(21hOn)(31hOn)(1phO2021-12-1328l数值方法具有越高阶的精度，得到的解就越精确。数值方法具有越高阶的精度，得到的解就越精确。我们发现，向前和

28、向后欧拉公式各用了区间一个我们发现，向前和向后欧拉公式各用了区间一个端点的导数，具有端点的导数，具有1阶精度，而梯形和改进欧拉阶精度，而梯形和改进欧拉公式用了两个端点的导数取平均，得到了公式用了两个端点的导数取平均，得到了2阶精阶精度。因此就引导人们利用内的若干点的导数，对度。因此就引导人们利用内的若干点的导数，对它们作线性组合得到平均斜率，由此得到更高阶它们作线性组合得到平均斜率，由此得到更高阶的精度，这就是的精度，这就是龙格龙格-库塔方法库塔方法的基本思路。的基本思路。l在在MATLAB软件中含有数值求解的系统函数，其软件中含有数值求解的系统函数，其实现原理就是龙格库塔方法。实现原理就是龙

29、格库塔方法。 2021-12-13292.4 Matlab软件求解软件求解 l2.4.1 解析解解析解l用MATLAB命令ldsolve(eqn1,eqn2, .) 求常微分方程（组）的解析解。l其中eqni表示第i个微分方程，Dny表示y的n阶导数,默认的自变量为t。2021-12-1330l1. 微分方程微分方程l例例2.6 求解一阶微分方程 l(1) 求通解l输入：ldsolve(Dy=1+y2) l输出：lans =ltan(t+C1) 21ydxdy（2）求特解输入：dsolve(Dy=1+y2,y(0)=1,x) 指定初值为1，自变量为x输出：ans =tan(x+1/4*pi)

30、2021-12-1331l例例2.7 求解二阶微分方程求解二阶微分方程 l输入输入:ldsolve(D2y+(1/x)*Dy+(1-(1/2)2/x2)*y=0,y(pi/2)=2,Dy(pi/2)=-2/pi,x) lans =2(1/2)*pi(1/2)/x(1/2)*sin(x) l化简输出结果，输入：化简输出结果，输入：lpretty(ans) 1/2 1/2l 2 pi sin(x)l -l 1/2l x2/1/2)2/(2)2/(0)(222 nyyynxyxyxxxysin2即：即： 2021-12-13322微分方程组微分方程组l例2.8 求解 df/dx=3f+4g; dg/

31、dx=-4f+3g。l（1）通解: lf,g=dsolve(Df=3*f+4*g,Dg=-4*f+3*g) lf =lexp(3*t)*(C1*sin(4*t)+C2*cos(4*t)lg =lexp(3*t)*(C1*cos(4*t)-C2*sin(4*t) l特解：lf,g=dsolve(Df=3*f+4*g,Dg=-4*f+3*g,f(0)=0,g(0)=1) lf =lexp(3*t)*sin(4*t)lg =lexp(3*t)*cos(4*t) 2021-12-13332.4.2 数值解数值解l在微分方程在微分方程(组组)难以获得解析解的情况下，可以用难以获得解析解的情况下，可以用M

32、atlab方便地求出数值解。格式为方便地求出数值解。格式为:lt,y = ode23(F,ts,y0)， t,y = ode45(F,ts,y0)l注意：注意：l微分方程的形式：微分方程的形式：y = F(t, y)，t为自变量，为自变量，y为因变为因变量（可以是多个，如微分方程组）；量（可以是多个，如微分方程组）；lt, y为输出矩阵，分别表示自变量和因变量的取值；为输出矩阵，分别表示自变量和因变量的取值；lF代表微分方程组的函数名（代表微分方程组的函数名（m文件，必须返回一个文件，必须返回一个列向量）；列向量）；lts的取法有几种，（的取法有几种，（1）ts=t0, tf 表示自变量的取值

33、表示自变量的取值范围，（范围，（2）ts=t0,t1,t2,tf,则输出在指定时刻则输出在指定时刻t0,t1,t2,tf处给出，（处给出，（3）ts=t0:k:tf,则输出在区间则输出在区间t0,tf的等分点给出；的等分点给出；ly0为初值条件。为初值条件。2021-12-1334l比如例例2.7的数值解：l解：令y1=y,y2=y1,将二阶微分方程转化为一阶微分方程组ly1=y2ly2=-y2/x+(n/x)2-1)y1l首先建立M-文件函数:l function f=jie3(x,y)l f=y(2);-y(2)/x+(1/2)2/x2-1)*y(1);l计算:lx,y=ode23(jie

34、3,pi/2,pi,2,-2/pi) 2/1/2)2/(2)2/(0)(222 nyyynxyxyx2021-12-1335lx =l 1.5708l 1.6074l 1.7645l 1.9215l 2.0786l 2.2357l 2.3928l 2.5499l 2.7069l 2.8640l 3.0211l 3.1416ly =l 2.0000 -0.6366l 1.9758 -0.6869l 1.8518 -0.8879l 1.6982 -1.0631l 1.5192 -1.2108l 1.3193 -1.3293l 1.1032 -1.4174l 0.8756 -1.4744l 0.64

35、16 -1.5002l 0.4060 -1.4951l 0.1735 -1.4602l 0.0002 -1.4140 2021-12-1336y1=y(:,1);y2=y(:,2);plot(x,y1,x,y2,r)2021-12-13371.41.61.822.22.42.62.833.2-2-1.5-1-0.500.511.522.5 范例：地中海鲨鱼问题范例：地中海鲨鱼问题 l意大利生物学家意大利生物学家Ancona曾致力于鱼类种群相互制约关曾致力于鱼类种群相互制约关系的研究，他从第一次世界大战期间（系的研究，他从第一次世界大战期间（1914年年7月月28日日1918年年11月月11日）

36、日）,地中海各港口捕获的几种鱼类地中海各港口捕获的几种鱼类捕获量百分比的资料中，发现鲨鱼的比例有明显增加捕获量百分比的资料中，发现鲨鱼的比例有明显增加（见下表）。（见下表）。年代年代1914191419151915191619161917191719181918百分比百分比11.911.921.421.422.122.121.221.236.436.4年代年代1919191919201920192119211922192219231923百分比百分比27.327.316.016.015.915.914.814.819.719.72021-12-1338l战争为什么使鲨鱼数量增加？是什么原因？战

37、争为什么使鲨鱼数量增加？是什么原因？l因为战争使捕鱼量下降，食用鱼增加，显然鲨鱼也随之增因为战争使捕鱼量下降，食用鱼增加，显然鲨鱼也随之增加。加。 l但为何鲨鱼的比例大幅增加呢？生物学家但为何鲨鱼的比例大幅增加呢？生物学家Ancona无法解无法解释这个现象，于是求助于著名的意大利数学家释这个现象，于是求助于著名的意大利数学家V.Volterra，希望建立一个食饵希望建立一个食饵捕食者系统的数学模型，定量地回答捕食者系统的数学模型，定量地回答这个问题。这个问题。l 1、符号说明：、符号说明：lx1(t), x2(t)分别是食饵、捕食者（鲨鱼）在分别是食饵、捕食者（鲨鱼）在t时刻的数时刻的数量；量

38、； lr1, r2是食饵、捕食者的固有增长率；是食饵、捕食者的固有增长率；l1是捕食者掠取食饵的能力，是捕食者掠取食饵的能力， l 2是食饵对捕食者的供养能力；是食饵对捕食者的供养能力；2021-12-1339l2、基本假设：、基本假设：l 捕食者的存在使食饵的增长率降低，假设捕食者的存在使食饵的增长率降低，假设x1降低的降低的程度与捕食者数量程度与捕食者数量x2成正比，即成正比，即l食饵对捕食者的数量食饵对捕食者的数量x2起到增长的作用，起到增长的作用， 其程度其程度与食饵数量与食饵数量x1成正比，即成正比，即l综合以上综合以上和和，得到如下模型：，得到如下模型：l模型一：不考虑人工捕获的情

39、况模型一：不考虑人工捕获的情况)(21111xrxdtdx)(12222xrxdtdx)()(1222221111xrxdtdxxrxdtdx2021-12-1340l该模型反映了在没有人工捕获的自然环境中食饵与捕食者该模型反映了在没有人工捕获的自然环境中食饵与捕食者之间的制约关系，没有考虑食饵和捕食者自身的阻滞作用，之间的制约关系，没有考虑食饵和捕食者自身的阻滞作用，是是Volterra提出的最简单的模型。提出的最简单的模型。l给定一组具体数据，用给定一组具体数据，用matlab软件求解。软件求解。 l食饵：食饵： r1= 1, 1= 0.1, x10= 25; l捕食者捕食者(鲨鱼鲨鱼)：r2=0.5, 2=0.02, x20= 2;l编制程序如下编制程序如下l1、建立、建立m-文件文件shier.m如下：如下：l function dx=shier(t,x) l dx=zeros(2,1); l dx(1)=x(1)*(1-0.1*x(2); l dx(2)=x(2)*(-0.5+0.02*x(1);2021-12-1341l2、在命令窗口执行如下程序：、在命令窗口执行如下程序： l t,x=ode4