微分方程模型与一阶常微分方程初值问题数值解

1、 数 学 建 模 电 子 教 案第 3 教学周/星期 五 /第 7、8 节/第 5 次课课 题第三章 微分方程模型与一阶常微分方程初值问题数值解§3.1 一阶微分方程初值问题数值解 §3.2猪的最佳销售时机教学内容1.常微分方程的两个模型2.一阶常微分方程初值问题数值解法3.猪的最佳销售时机问题的模型及实验教学目标1. 了解一阶微分方程的初值问题的两个数值解法：欧拉方法、Runge-kutta(龙格库塔)方法。2. 会利用变化率分析并建立微分方程模型。3. 会用软件Mathematica和MATLAB求解微分方程模型。教学重点1.掌握微分方程数值解法得基本思想.2.了解欧拉

2、方法、利用改进的欧拉公式解一阶微分方程的初值问题的数值解教学难点Runge-kutta (龙格库塔)方法双语教学内容、安排 Differential equation; 微分方程 numerical value solution; 数值解 教学手段、措施板书、 结合多媒体教学作业、后记P69,2教学过程及教学设计备注 §3.1 一阶微分方程初值问题数值解一、两个模型1、饿狼追兔问题现有一只兔子，一只狼，兔子位于狼的正西100米处。假设兔子与狼同时发现对方并一起起跑，兔子往正北60米处的巢穴跑，而狼在追兔子，已知兔子、狼是匀速跑且狼的速度是兔子的两倍。问题是兔子能否安全回到巢穴？解 首

3、先建立坐标系，兔子在O处，狼在A处。由于狼要盯着兔子追，所以狼行走的是一条曲线，且在同一时刻，曲线上狼的位置与兔子的位置的连线为曲线上该点处的切线。设狼的行走轨迹是y=f(x)，则有 ，又因狼的速度是兔子的两倍，所以在相同时间内狼走的距离为兔子走的距离的两倍。假设在某一时刻，兔子跑到(0,h)处，而狼在(x,y)处，则有整理得到下述模型这属于可降阶的二阶微分方程，解得狼的行走轨迹因，所以狼追不上兔子。2、尸体冷却模型受害者的尸体于晚上7:30被发现，法医于晚上8:20赶到凶案现场，测得尸体温度为32.6；一小时后，当尸体即将被抬走时，测得尸体温度为31.4，室温在几个小时内始终保持21.1。此

4、案最大的嫌疑犯张某声称自己是无罪的，并有证人说：“下午张某一直在办公室上班，5:00时打完电话后就离开了办公室”。从张某到受害者家（凶案现场）步行需5分钟，现在的问题是，张某不在凶案现场的证言能否被采信，使他排除在嫌疑犯之外。解：首先应确定凶案的发生时间，若死亡时间在下午5点5分之前，则张某就不是嫌疑犯，否则不能将张某排除。设T(t)表示t时刻尸体的温度，并记晚上8:20为t=0，则T(0)=32.6，T(1)=31.4。假设受害者死亡时体温是正常的，即T=37是要确定受害者死亡的时间，也就是求T(t)=37的时刻，进而确定张某是否是嫌疑犯。人体体温受大脑神经中枢调节。人死亡后体温调节的功能消

5、失，尸体的温度受外界环境温度的影响。假设尸体温度的变化率服从牛顿冷却定律，即尸体温度的变化律与他同周围的温度差成正比。即 k是常数，分离变量积分得： 由T(0)=21.1+a=32.6 得a=11.5；由T(1)=21.1+ae-k=31.4 得e-k115/103，即k=0.11，所以 T(t)=21.1+11.5e-0.11t .当T=37时，有t=-2.95 小时-2小时57分，8小时20分2小时57分5小时23分。即死亡时间大约在下午5:23，因此张某不能被排除在嫌疑犯之外。二、一阶微分方程初值问题数值解1、导入课程：微分方程的定解问题中着重考虑的一阶方程的初值问题函数满足利普希茨条件

6、：（3-1）的解存在并且唯一。常微分方程的解析方法只能用来求解一些特殊类型的方程，实际问题中主要靠数值解法。所谓数值解法，就是寻求解y(x)在一系列离散节点上的近似值相邻两个节点的距离称为步长，节点为初值问题（3-1）的数值解法都采用进步式，即只要给出用已知信息就能给出计算的递推公式。2、欧拉方法的递推公式：它的基本思想是在小区间上用差商代替导数，而方程右端函数中的在小区间的端点上取值，得到方程的近似表达式，称为欧拉公式。（1） 向前欧拉公式：n=0,1,2. （1）被称作向前欧拉公式或显式欧拉公式。（2）向后欧拉公式： （2）被称作向后欧拉公式或隐式欧拉公式。（3）梯形公式：n=0,1,2,

7、 (3)被称作梯形公式。（4）改进的欧拉公式: (4)被称作改进的欧拉公式。例1、求解初值问题 （3.3）解（1）向前欧拉公式的方法具体形式为：取步长h=0.1，计算结果见表31。0.11.10001.09540.61.50901.48320.21.19181.18320.71.58031.54920.31.27741.26490.81.64981.61250.41.35821.34160.91.71781.67330.51.43511.41421.01.78481.7321（2）用改进的欧拉公式为仍取h=0.1，计算结果比较见表32 0.11.09591.09540.61.48601.483

8、20.21.18411.18320.71.55251.54920.31.26621.26490.81.61531.61250.41.34341.34160.91.67821.67330.51.41641.41421.01.73791.73213、Runge-kutta(龙格库塔)方法它的基本思想是将在x处进行泰勒展开，并取其前面几项来近似而得到公式：由上式产生的迭代公式：若，则称以上述公式为P阶公式，P的大小反映出精度的高低。二阶Rungekutta公式(5)其中为待定系数，它满足：上式被称作是二阶RK公式。四阶Rungekutta公式 (6)其中待定系数共13个，由于计算复杂，直接给出一组的

9、值，得 (7)称作四阶R-K公式。例3.2 设取步长h=0.2，从x=0直到x=1用四阶龙格库塔方法求解初值问题(3-3)。解 四阶龙格-库塔公式表3-3列出计算结果，表中仍表示准确解。 表3-30.21.18321.18320.41.61251.34160.61.48331.48320.81.61251.61251.01.73211.7321比较例3.3和例3.2的计算结果，显然以龙格库塔方法的精度为高。值得指出：龙格库塔方法的推导基于泰勒展开方法，因而它要求所求的解具有较好的光滑性质。如果解的光滑性差，使用四阶龙格库塔方法求得的数值解，其精度可能不如改进的欧拉法，实际计算时，我们应当针对问

10、题的具体特点选择合适的算法。 §3.2猪的最佳销售与天然气储量问题一、猪的最佳销售时机问题：1、问题的提出：对于猪的商业性饲养和销售，人们总是希望获得最大的利润，在市场需求不变的情况下，如果我们不考虑猪的饲养技术、水平，猪的类型等因素的影响，那么影响销售利润的主要因素，就是销售时机问题，由于随着猪的生长，单位时间消耗的饲养费用逐渐增多，而猪的体重增长却逐渐变慢，因此对猪的饲养时间过长是不合算的。假定一头猪在开始饲养时的重量为x0，在饲养后任意时刻t的重量为x(t)，对于某一品种的猪，它的最大重量假定为X0，猪的最小出售体重为xs，相应的饲养时间为ts。一头猪从开始饲养到时刻t所需的费

11、用为y(t)，同时我们假定反映猪体重变化速度的参数为，猪在达到最大体重后，单位时间的饲养费为y，反映饲养费用变化大小的参数为，请根据上面的假设，建立起猪的最佳销售时机的数学模型，并用下面所给的数据验证你的模型。假设X0=200（kg），xs=75（kg），=0.5（kg/天），猪的市场销售价设为c=6元/kg，=1.5（元/天），=1（元/天），x0=5kg。2、问题分析：由于猪在进行饲养时已具有一定的体重，而其体重的增加随饲养时间的延长逐渐减慢，因此由Logistic模型可得；又由于猪的体重增加，单位时间消耗的饲养费用就越多，达到最大体重后，饲养费用为常数，所以有，因此，得到微分方程：求解可

12、得： （1）养猪能否获利，主要看猪从出生到ts时，如果出售是否可以获利，因此，获利的充要条件为： （2）其中c0为仔猪的价格。由（1）式可得：，解之可得：，将（1），（2）式代入可得： （3）所以只要（3）式成立，饲养就会获利。设猪的最佳出售时机为t*，由（1）式求导可得： （4）由盈亏平衡原理：即单位时间内由猪增加体重所获得的利润与消耗的饲养费用相等，可得由（4）式可得：，解之可得。若时，故猪应在时出售。若时，故猪应在时出售（因为猪必须长到xs）。猪的最佳销售时机问题的计算。下面给出Mathematica计算程序：%最佳销售时机ch411%文件名：ch411.mX=200.0;xs=75.0

13、;x0=5.0;c=6.0;alpha=0.5;bet=1.0;gama=1.5;temp=gama*X/(X-xs)-(c*alpha+bet);t=iftemp<0,X/alpha*Log(c*alpha+bet)*(X-x0)/(gama*X),X/alpha*Log(X-x0)/(X-xs)执行后输出Out1=True,382.205,177.874即猪的最佳出售时间为饲养到382天左右。用多媒体某些类型的导弹对目标追击的数学模型于此模型类似。牛顿冷却定律：即尸体温度的变化律与他同周围的温度差成正比。记将向前和向后欧拉公式加以平均得到。迭代过程分为两步：由向前欧拉公式算出的预测值

14、再把它带入梯形公式的右端，作为较正。此方程的解析解为，按这个解析式子算出的准确值同近似值相比较欧拉方法的精度很差。改进欧拉法明显改善了精度在区间内仍取2个点，仿照改进的欧拉公式得注意：虽然四阶龙格库塔方法的计算量)比改进的欧拉方法大一倍，但由于这里放大了步长(h=0.2)，表33和表32所耗费的计算量几乎相同。利用多媒体此处介绍Logistic模型此处用Mathematica计算程序演示课 题第三章 微分方程模型与一阶常微分方程初值问题数值解§3.3天然气储量问题 §3.4 最优捕鱼策略教学内容微分方程的模型与实验： 1.天然气储量问题的模型与实验 2.最优捕鱼策略教学目标

15、1、会利用变化率分析并建立微分方程模型。2、会用软件Mathematica和MATLAB求解微分方程模型。教学重点会利用变化率分析并建立微分方程模型。教学难点引入变量，并得出变量之间的关系双语教学内容、安排 教学手段、措施 多媒体教学，结合板书作业、后记P69 3、4教学过程及教学设计备注§3.3天然气储量问题一、 天然气储量问题1 、 问题的提出:天然气资源是现代社会重要的基础能源之一，应合理的开发和利用。对开发公司而言，准确地预测天然气的产量和可采储量，始终是一项重要而又艰难的工作，下面是某天然气公司在19571976年20年间对某气田产量的统计资料。 某气田1957年至1976

16、年产量表年 度1957 1958 1959 1960 1961 1962 1963 1964 1965 1966产量（108m3） 19 43 59 82 92 113 138 148 151 157年 度1967 1968 1969 1970 1971 1972 1973 1974 1975 1976产量（108m3） 158 155 137 109 89 79 70 60 53 45试根据所给的数据资料，建立起该气田产量的预测模型，并验证所建立模型的合理性。2、模型假设及符号说明（1）假设该气田的产量是连续的，没有阶段性停产现象。（2）假设所提供的数据是正常生产情况下，气田的产量。（3）假

17、设没有因意外事故或自然灾害造成停产或减产的情况。（4）假设r（t）为天然气产量的增长率。（5）假设N（t）为天然气田的累积产量。（6）假设Q为气田的年产量。（7）假设气田的可采储量为NR，相对应的开发时间tR。3、问题分析及数学模型根据所给的实际问题，预测气田的产量和可采储量，在这方面，目前国内外的方法很多，但各种预测方法中有一种简单而实用的指数增长模型，它是借鉴英国人口学家Malthus（马尔萨斯）于1798年提出的人口增长模型，而得到的。若假设天然气产量的增长率为r（t），它是时间t的连续函数。气田的累积产量设为N（t），则它们满足如下的关系：而气田的年产量，于是上述方程变为：有统计资料显

18、示，气田的每年产量与累积产量之比与气田的开发时间t存在如下关系，即 （5）或写成：，其中：假设气田的可采储量为NR，相对应的开发时间为tR，由上面的分析可得到方程：解之可得 ： （6）对上式求导预测模型为 （7）4 、 模型的分析与计算（1）模型（3）中a , b的计算由于，所以关键问题在于求出A，B的值。由（1）式，设，其中为第i年的产量，为第i年之前的累积产量，时间t以年为单位，则由所给数据可得根据线性回归的最小二乘估计，令为使L（A，B）最小，取L分别关于A，B的偏导数，并令它们为0。解此方程可得 其中： 从而。（6）模型（7）中NR的计算对（2）式两边取常用对数可得 ： （8）其中：。

19、由（8）式和所给的数据，建立回归方程（同上），可求得、，从而计算出油田的可采储量（略）。（3）模型的求解将根据上述求得的a，b，NR的值代入模型（6）（7）式便可计算出相应年份累积产量N（t）和年产量Q的预测值。气田储量的mathematica计算程序：%气田储量ch42%文件名：ch42.mdata1=19.0,43.0,59.0,82.0,92.0,113.0,138.0,148.0,151.0,157.0,158.0,155.0,137.0,109.0,89.0,79.0,70.0,60.0,53.0,45.0;data2=Table0,n,20;i=2;Fordata21=data11

20、,i<=20,i+,data2i=data2i-1+data1idata3=Table0,m,20,n,2;Fori=1,i<=20,i+,data3i=i,Log10,data1i/data2iFitdata3,1,t,taa=-0.0215995;bb=0.0809426;a=10aa;b=Log10.0*bb;c=a/b;data4=Table0,m,20,n,2;Fori=1,i<=20,i+,data4i=Exp-b*i,Log10,data2iFitdata4,1,x,xalpha=3.36832;bet=2.35678;Nr=10alpha;Qpt_:=a*Nr

21、*Exp-c*Exp-b*t-b*t;Npt_:=Nr*Exp-c*Exp-b*t;data5=TableQpt,t,1,20;data6=TableNpt,t,1,20;compdata1=Table0,m,20,n,2;Fori=1,i<=20,i+,compdata1i=data1i,data5i;compdata2=Table0,m,20,n,2;Fori=1,i<=20,i+,compdata2i=data2i,data6i;f1=ListPlotdata1;f2=PlotQpt,t,0,20;f3=ListPlotdata2;f4=PlotNpt,t,0,20;Show

22、f1,f2Showf3,f4MatrixFormcompdata1MatrixFormcompdata2执行后输出A=-0.0215995 B=0.0809426NR=10alpha alpha=3.36832实际值与预测值对照表年份T（a）Q（108m3/a）NP（108m3）实际值预测值实际值预测值19571958195919601961196219631964196519661967196819691970197119721973197419751976123456789101112131415161718192019.043.059.082.092.0113.0138.0148.015

23、1.0157.0158.0155.0137.0109.089.079.070.060.053.045.026.67445.45768.60393.527117.186136.899150.897152.492159.935156.116148.242137.580125.282112.29899.35186.94775.40964.91455.53447.26519.062.0121.0203.0295.0408.0546.0694.0845.01002.01160.01315.01452.01561.01650.01729.01799.01859.01912.01957.069.365126

24、.117207.160312.745440.207584.626739.855899.5521057.9701210.4301353.5201485.0501603.8601709.6601802.7501883.8401953.9102014.0402065.350从结果看计算比较精确。§3.3 最优捕鱼策略1 问题的提出：这是1996年全国大学生数学建模竞赛的A题，问题如下：为保护人类赖以生存的自然环境，可再生资源（如渔业、林业等资源）的开发必须适度。一种合理、简化的策略是，在实现可持续收获的前提下，追求最大产量或最佳效益。考虑对某种鱼的最优捕捞策略：假设这种鱼分4个年龄组：称一

25、龄鱼、二龄鱼、三龄鱼、四龄鱼。各年龄组每条鱼的平均重量分别为5.07，11.55，17.86，22.99（克）；各年龄组鱼的自然死亡率均为0.8（1/年）；这种鱼为季节性集中产卵繁殖，平均每条4龄鱼的产卵量为1.109×105（个），3龄鱼的产卵量为这个数的一半，2龄鱼和1龄鱼不产卵，产卵和孵化期为每年的最后4个月；卵孵化并成活为1龄鱼，成活率为（1龄鱼条数与产卵总量n之比）1.22×1011/(1.22×1011+n)。渔业管理部门规定，每年只允许在产卵孵化期的前8个月内进行捕捞作业。如果每年投入的捕捞能力（如鱼船数、下网次数等）固定不变，这时单位时间捕捞量将与

26、各年龄组鱼群条数成正比，比例系数不防称为捕捞强度系数。通常使用13mm网眼的拉网，这种网只能捕捞3龄鱼和4龄鱼，其两个捕捞强度系数之比为0.42:1。渔业上称这种方式为固定努力量捕捞。（1）建立数学模型分析如何实现可持续捕捞（即每年开始捕捞时渔场中各年龄组鱼群条数不变），并且在此前提下得到最高的年收获量（捕捞总重量）。（2）某渔业公司承包这种鱼的捕捞业务5年，合同要求5年后鱼群的生产能力不能受到太大破坏。已知承包时各年龄组鱼群数量分别为：122，29.7，10.1，3.29（×109条）。如果仍用固定努力量的捕捞方式，该公司采用怎样的策略才能使总收获量最高。2、 模型假设及符号说明（

27、1）假设只考虑一种鱼的繁殖和捕捞，鱼群增长过程中不考虑鱼的迁入与迁出。（2）假设各年龄组的鱼在一年内的任何时间都会发生自然死亡，产卵可在后四个月内任何时间发生。（3）假设3、4龄鱼全部具有生殖能力，或者虽然雄鱼不产卵，但平均产卵量掩盖了这一差异。（4）假设产卵期鱼的自然死亡率发生于产卵之后。（5）假设各年龄组的鱼经过一年后，即进入高一级的年龄组，但4龄鱼经过一年后仍视为4龄鱼。（6）假设对鱼的捕捞用固定努力量捕捞方式，每年的捕捞强度系数保持不变，且捕捞只在前八个月进行。（7）假设t时刻i龄鱼的数量为。（8）假设第k年初i龄的数量为；第k年底i龄鱼的数量为（i=1，2，3，4）。（9）假设鱼的自

28、然死亡率为r；4龄鱼的平均产卵量为c。（10）假设第i龄鱼的平均重量为Mi（i=1，2，3，4）。（11）假设第k年度鱼的产卵总量为。（12）假设对第i龄鱼的捕捞强度系数为bi；对i龄鱼的年捕捞量为ai（i=1，2，3，4）。（13）假设年总收获量为M，即M=M3a3+M4a4。（14）假设5年的总收获量为MM，即。3 问题分析及数学模型由已知条件，可得，（E为捕捞努力量），r为自然死亡率，在t，t+t内，根据死亡率的定义，由于不捕捞1、2龄鱼，所以变形可得 （1）解得。对于3、4龄鱼由于捕捞在前8个月进行，因此，前8个月，捕捞与死亡均影响鱼的变化，因而微分方程变形为 （2）其中由（2）式解得

29、 （3）令为i龄鱼在k年底时的数量，为i龄鱼在k年初时的数量，得到1、2龄鱼k年末的数量为 （4）对于3，4龄鱼，在第8个月末数量由（3）式可得 （）在后4个月，对于3、4龄鱼，只有死亡率起作用，因而微分方程为 解得：因而，3、4龄在k年末的数 （）将（4）式与（）合在一起，得到在k年底各龄鱼的数量为 （5）由假设，到年底第i龄鱼全部转化为i+1龄鱼（i=1，2，3），同时由卵孵化产生1龄鱼，于是得到 （6）其中，为年度总产卵量，且 （7）此外，我们还求得每年对3、4龄鱼的总捕捞重量为 其中 4 模型的分析与计算（1）年度产量最优模型及其计算为了实现可持续捕捞，即要求，在此前提下获得最高年收获

30、量。结合基本模型，即可得到年度产量最优模型。利用，也就是，等与时间k无关，化简得下面给出MATLAB计算程序：%最优捕鱼策略ch331%文件名：ch331.mx=sym('x'); b3=0.42*x;d=1.22*1011;r=0.8;q=d*exp(-(3*r+2/3*b3)*(32529.55*exp(r)+65059.1*exp(-2/3*x)/(1-exp(-(r+2/3*x);N10=d*q/(d+q);N40=d*q/(d+q)*exp(-(3*r+2/3*b3)/(1-exp(-(r+2/3*x);a3=b3/(r+b3)*(1-exp(-2/3*(r+b3)*

31、d*q/(d+q)*exp(-2*r);a4=x/(r+x)*(1-exp(-2/3*(r+x)*d*q/(d+q)*exp(-(3*r+2/3*b3)/(1-exp(-(r+2/3*x);M3=17.86;M4=22.99;M=M3*a3+M4*a4;M1=-M;M10=char(M);M11=char(M1);fplot(M10,0,100)b4=fmin(M11,0,100);b3=0.42*b4;d=1.22*1011;r=0.8;q=d*exp(-(3*r+2/3*b3)*(32529.55*exp(r)+65059.1*exp(-2/3*b4)/(1-exp(-(r+2/3*b4);N10=d*q/(d+q);N40=d*q/(d+q)*exp(-(3*r+2/3*b3)/(1-exp(-(r+2/3*b4);a3=b3/(r+b3)*(1-exp(-2/3*(r+b3)*d*q/(d+q)*exp(-2*r);a4=b4/(r+b4)*(1-exp(-2/3*(r+b4)*d*q/(d+q)*exp(-(3*r