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

微分方程模型与一阶常微分方程初值问题数值解 数 学 建 模 电 子 教 案 第 3 教学周/星期 五 /第 7、8 节/第 5 次课 第三章 微分方程模型与一阶常微分方程初值问题数值解 课 题?3.1 一阶微分方程初值问题数值解 ?3.2猪的最佳销售时机 1.常微分方程的两个模型 教学内容 2.一阶常微分方程初值问题数值解法 3.猪的最佳销售时机问题的模型及实验 1. 了解一阶微分方程的初值问题的两个数值解法:欧拉方法、Runge-kutta(龙 格,库塔)方法。 教学目标 2. 会利用变化率分析并建立微分方程模型。 3. 会用软件Mathematica和MATLAB求解微分方程模型。 1.掌握微分方程数值解法得基本思想. 教学重点 2.了解欧拉方法、利用改进的欧拉公式解一阶微分方程的初值问题的数值解 教学难点 Runge-kutta (龙格,库塔)方法 双语教学 Differential equation; 微分方程 内容、安排 numerical value solution; 数值解 教学手段、板书、 结合多媒体教学 P69,2 作业、后记 教学过程及教学设计 备注 ?3.1 一阶微分方程初值问题数值解 一、两个模型 1、饿狼追兔问题 现有一只兔子，一只狼，兔子位于狼的正西100米处。假设兔子与狼同时发现对方 并一起起跑，兔子往正北60米处的巢穴跑，而狼在追兔子，已知兔子、狼是匀速跑且 狼的速度是兔子的两倍。问题是兔子能否安全回到巢穴, y 解 首先建立坐标系，兔子在O处， B -60 狼在A处。由于狼要盯着兔子追，所以 狼行走的是一条曲线，且在同一时刻， y=f(x) 曲线上狼的位置与兔子的位置的连线为 曲线上该点处的切线。设狼的行走轨迹 h C(x,y) ,y,0y,0是y=f(x)，则有 ， x x,100x,100 O A(100,0) 又因狼的速度是兔子的两倍，所以 在相同时间内狼走的距离为兔子走的距离的两倍。假设在某一时刻，兔子跑到 用多媒体 (0,h)处，而狼在(x,y)处，则有 1 数 学 建 模 电 子 教 案 第 3 教学周/星期 五 /第 7、8 节/第 5 次课 hy,,,()fx,,某些类型的导0,x ,1002弹对目标追击,,21()hfxdx,，,x,的数学模型于 整理得到下述模型 此模型类似。 2,,,, 2()1()xfxfx,，, , ,ff(100)0,(100)0,,,, 这属于可降阶的二阶微分方程，解得狼的行走轨迹 311200 22 fxxx()10,,，303 200因，所以狼追不上兔子。 f(0)60,, 3 2、尸体冷却模型 受害者的尸体于晚上7:30被发现，法医于晚上8:20赶到凶案现场，测得尸体温度 为32.6?;一小时后，当尸体即将被抬走时，测得尸体温度为31.4?，室温在几个小 时内始终保持21.1?。此案最大的嫌疑犯张某声称自己是无罪的，并有证人说:“下牛顿冷却定律: 午张某一直在办公室上班，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?的时刻，进而确定张某是否是嫌疑犯。 人体体温受大脑神经中枢调节。人死亡后体温调节的功能消失，尸体的温度受外界 环境温度的影响。假设尸体温度的变化率服从牛顿冷却定律，即尸体温度的变化律与他 同周围的温度差成正比。即 dT k是常数， ,,,kT(21.1)dt ,kt分离变量积分得: Ttae()21.1,， 由T(0)=21.1+a=32.6 得a=11.5;由T(1)=21.1+ae-k=31.4 得e-k,115/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、导入课程:微分方程的定解问题中着重考虑的一阶方程的初值问题 ,y,f(x,y) , (3,1), y(x),y00, 2 数 学 建 模 电 子 教 案 第 3 教学周/星期 五 /第 7、8 节/第 5 次课 函数满足利普希茨条件: f(x,y)f(x,y),f(x,y),Ly,y(3,2)(3-1)的解存在并且唯一。 y,y(x) 常微分方程的解析方法只能用来求解一些特殊类型的方程，实际问题中主要靠数值 解法。所谓数值解法，就是寻求解y(x)在一系列离散节点 x,x,？,x,x,？12nn，1上的近似值相邻两个节点的距离称为步长，节点为 y,y,？,y,y,？.hxx,,12nn，1nn，1 初值问题(3-1)的数值解法都采用进步式，即只要给出用 x,x，nh,n,0,1,2,？.n0 已知信息就能给出计算的递推公式。 yy,y,y,？nn,1n,2n，1 2、欧拉方法的递推公式: yxyx(),()/n，1n它的基本思想是在小区间上用差商代替导数，而方,,x,xxy nn，1h 程右端函数中的在小区间的端点上取值，得到方程的近似表达式，称为欧拉,,x,xnn，1 公式。 (1) 向前欧拉公式: n=0,1,2„„. (1) y,y，hf(x,y)n，1nnn 记y(x),ynn被称作向前欧拉公式或显式欧拉公式。 (2)向后欧拉公式: (2) yyhfxyn,，,(,)0,1,2nnnn，，，111 被称作向后欧拉公式或隐式欧拉公式。 (3)梯形公式:: h n=0,1,2,„ „(3) ,,y,y，f(x,y)，f(x,y)将向前和向后n，1nnnn，1n，12 欧拉公式加以被称作梯形公式。 平均得到。 (4)改进的欧拉公式: ,,y,y，hf(x,y)n，1nnn迭代过程分为, (4) n,0,1,2......,,h，，两步:由向前欧y,y，f(x,y)，f(x,y),n，1nnnn，!n，1,,2，，拉公式算出, 的预测值y被称作改进的欧拉公式。 n，1, 再把它带y n，1 入梯形公式的x2,/yyx,,0,,1，，,右端，作为较例1、求解初值问题 (3.3) y, 正。 ,y(0),1, x2 n解(1)向前欧拉公式的方法具体形式为: y,y，hy,()n，1nny n 取步长h=0.1，计算结果见表3,1。 此方程的解析 解为 ，y,1，2x 3 数 学 建 模 电 子 教 案 第 3 教学周/星期 五 /第 7、8 节/第 5 次课 按这个解析式 子算出的准确 yyxy(x)xy(x)nnnnnn 同近值y(x)0.1 1.1000 1.0954 0.6 1.5090 1.4832 n 似值相比较y0.2 1.1918 1.1832 0.7 1.5803 1.5492 n 欧拉方法的精0.3 1.2774 1.2649 0.8 1.6498 1.6125 度很差。 0.4 1.3582 1.3416 0.9 1.7178 1.6733 0.5 1.4351 1.4142 1.0 1.7848 1.7321 (2)用改进的欧拉公式为 h,()yykk,，， nn，112,2, ,22xxnn，1,,()kykyhk,,,，, 121nnyyhk，,nn1,改进欧拉法明 仍取h=0.1，计算结果比较见表3,2 显改善了精度 yyxy(x)xy(x) nnnnnn 0.1 1.0959 1.0954 0.6 1.4860 1.4832 0.2 1.1841 1.1832 0.7 1.5525 1.5492 0.3 1.2662 1.2649 0.8 1.6153 1.6125 0.4 1.3434 1.3416 0.9 1.6782 1.6733 0.5 1.4164 1.4142 1.0 1.7379 1.7321 3、Runge-kutta(龙格,库塔)方法 它的基本思想是将在x处进行泰勒展开，并取其前面几项来近似 y(x，h)y(x，h) 而得到公式: y(x，h),y(x)，h,,(x,y(x),h) 由上式产生的迭代公式: y,y，h,(x,y,t)n,0,1,2,......n，1nnnp，1若 ，则称以上述公式为P阶公式，P,,y(x，h),y(x)，h,(x,y(x),h),o(h) 的大小反映出精度的高低。 二阶Runge,kutta公式 yy,h(k,k),，，,n1n1122， ,kf(x,y),,1nn (5) ,kf(x,h,y,hk),0,,,1,，，,, 2nn1, ,1,,,,其中为待定系数，它满足: ,,,,,,,，,1,,,,1 12122,2 上式被称作是二阶R,K公式。 四阶Runge,kutta公式 4 数 学 建 模 电 子 教 案 第 3 教学周/星期 五 /第 7、8 节/第 5 次课 yyhkkkk,，，，，(),,,,,nn，111223344,在区间kfxykfxhyhk,,，，(,),(,),,,12111nnnn (6) ,内仍,,x,xnn，1kfxhyhkhk,，，，(,),,,322132nn,取2个点，仿照,kfxhyhkhkhk,，，，，(,),,,,43415263nn,改进的欧拉公其中待定系数共13个，由于计算复杂，直接给出一组的值，得 ,,,,,,,,,,iiiiii式得 h, yykkkk,，，，，(22)nn，11234,6 ,hh, kfxykfxyk,,，，(,),(,) (7) ,121nnnn22, hh, kxykkxhyhk,，，,，，(,),(,)3243nnnn,22, 称作四阶R-K公式。 例 3.2 设取步长h=0.2，从x=0直到x=1用四阶龙格,库塔方法求解初值问题(3-3)。 解 四阶龙格-库塔公式 , , ,h ,yyKKKK,，，，，(22)nn，112346 , ,22xxh，h, nn ,,,，,KyKyK,,121nnhy2 n,，yKn12, ,22()xhxh，，hnn KyKKyhK,，,,，,,,3243nnh2yhK，,n3 yK，n2,,2 表3-3列出计算结果，表中仍表示准确解。 yy(x)nn 注意:虽然四阶 表3-3 龙格,库塔方 yxy(x)nnn 法的计算量)比0.2 1.1832 1.1832 改进的欧拉方0.4 1.6125 1.3416 法大一倍，但由0.6 1.4833 1.4832 于这里放大了0.8 1.6125 1.6125 步长(h=0.2)，1.0 1.7321 1.7321 表3,3和表3比较例3.3和例3.2的计算结果，显然以龙格,库塔方法的精度为高。 ,2所耗费的计值得指出:龙格,库塔方法的推导基于泰勒展开方法，因而它所求的解具有较好的 算量几乎相同。 光滑性质。如果解的光滑性差，使用四阶龙格,库塔方法求得的数值解，其精度可能不 如改进的欧拉法，实际计算时，我们应当针对问题的具体特点选择合适的算法。 5 数 学 建 模 电 子 教 案 第 3 教学周/星期 五 /第 7、8 节/第 5 次课 ?3.2猪的最佳销售与天然气储量问题 一、猪的最佳销售时机问题: 1、问题的提出:对于猪的商业性饲养和销售，人们总是希望获得最大的利润，在 市场需求不变的情况下，如果我们不考虑猪的饲养技术、水平，猪的类型等因素的影响， 那么影响销售利润的主要因素，就是销售时机问题，由于随着猪的生长，单位时间消耗 的饲养费用逐渐增多，而猪的体重增长却逐渐变慢，因此对猪的饲养时间过长是不合算 的。 利用多媒体 假定一头猪在开始饲养时的重量为x，在饲养后任意时刻t的重量为x(t)，对于0 某一品种的猪，它的最大重量假定为X，猪的最小出售体重为x，相应的饲养时间为0s t。一头猪从开始饲养到时刻t所需的费用为y(t)，同时我们假定反映猪体重变化速度s 的参数为α，猪在达到最大体重后，单位时间的饲养费为y，反映饲养费用变化大小的 参数为λ，请根据上面的假设，建立起猪的最佳销售时机的数学模型，并用下面所给的 数据验证你的模型。 此处介绍 假设X=200(kg)，x=75(kg)，α=0.5(kg/天)，猪的市场销售价设为c=6元/kg，0sLogistic模型 γ=1.5(元/天)，λ=1(元/天)，x=5kg。 0 2、问题分析:由于猪在进行饲养时已具有一定的体重，而其体重的增加随饲养时 dxx,,间的延长逐渐减慢，因此由Logistic模型可得;又由于猪的体重增加，,(1) dtX0 单位时间消耗的饲养费用就越多，达到最大体重后，饲养费用为常数γ，所以有 dyx,,,，因此，得到微分方程: ,,(1) dtX0 dxx,,,(1,), dtX0, ,dyx,,,,(1,) , dtX0, ,x(0),x0, (0),0y, , ,t,X0,,,,x(t)X(Xx)e 000,求解可得: (1) ,, ,t,,X0,y(t),t,(X,x)(1,e) 00,,, 养猪能否获利，主要看猪从出生到t时，如果出售是否可以获利，因此，获利的s 充要条件为: xc,xc，y(t) (2) s00s 其中c为仔猪的价格。 0 , ,tX0由(1)式可得:， x,X,e(X,x) 000s 6 数 学 建 模 电 子 教 案 第 3 教学周/星期 五 /第 7、8 节/第 5 次课 XX,x000解之可得:t,，将(1)，(2)式代入可得: ln s,X,x0s X,x0xc,xc，x,x,X (3) ,(),(),ln ss0000X,xs 所以只要(3)式成立，饲养就会获利。 设猪的最佳出售时机为t*，由(1)式求导可得: , ,,txdxX00,,,,(1)e dtX,0 (4) ,,,t,xdy X00,,,,(1,)e,dtX 0, 由盈亏平衡原理:即单位时间内由猪增加体重所获得的利润与消耗的饲养费用相 等，可得 dxdy c, dtdt ,,,t,t00xx00XX 00cee由(4)式可得:，解之可得,,(1,),,,,(1,)XX 00 ,,X(c，)(X,x)000 。 t,ln0,,X0 ,,,XX(c，)(X,x)0000 若时，故猪应在时出售。 ,c,，,t,tt,,ln0s,,X,xX0s0 XX,x,X0000 t,,t,若时，故猪应在时出售(因为猪必,c,，,lnt,ts0s,X,xX,x0s0s 须长到x)。 s 此处用 猪的最佳销售时机问题的计算。 Mathematica计下面给出Mathematica计算程序: 算程序演示 %最佳销售时机ch411 %文件名:ch411.m X=200.0;xs=75.0;x0=5.0;c=6.0; alpha=0.5;bet=1.0;gama=1.5; temp=gama*X/(X-xs)-(c*alpha+bet); t=if[temp<0,X/alpha*Log[(c*alpha+bet)*(X-x0)/(gama*X)], X/alpha*Log[(X-x0)/(X-xs)]] 执行后输出 Out[1]=[True,382.205,177.874] 即猪的最佳出售时间为饲养到382天左右。 7 数 学 建 模 电 子 教 案 第 3 教学周/星期 五 /第 7、8 节/第 5 次课 第三章 微分方程模型与一阶常微分方程初值问题数值解 课 题 ?3.3天然气储量问题 ?3.4 最优捕鱼策略 微分方程的模型与实验: 1.天然气储量问题的模型与实验 教学内容 2.最优捕鱼策略 1、会利用变化率分析并建立微分方程模型。 教学目标 2、会用软件Mathematica和MATLAB求解微分方程模型。 教学重点 会利用变化率分析并建立微分方程模型。 教学难点 引入变量，并得出变量之间的关系 双语教学内 容、安排 教学手段、措 多媒体教学，结合板书 施 P69 3、4 作业、后记 教学过程及教学设计 备注 ?3.3天然气储量问题 一、 天然气储量问题 1 、 问题的提出:天然气资源是现代社会重要的基础能源之一，应合理的开发和利用。对开发公 司而言，准确地预测天然气的产量和可采储量，始终是一项重要而又艰难的工作，下面是某天然气公 司在1957~1976年20年间对某气田产量的统计资料。 某气田1957年至1976年产量表 年 度 1957 1958 1959 1960 1961 1962 1963 1964 1965 1966 83产量(10m) 19 43 59 82 92 113 138 148 151 157 年 度 1967 1968 1969 1970 1971 1972 1973 1974 1975 1976 83 产量(10m) 158 155 137 109 89 79 70 60 53 45 试根据所给的数据资料，建立起该气田产量的预测模型，并验证所建立模型的合理性。 介绍2、模型假设及符号说明 Malthus(马(1)假设该气田的产量是连续的，没有阶段性停产现象。 尔萨斯)人口(2)假设所提供的数据是正常生产情况下，气田的产量。 增长模型 8 数 学 建 模 电 子 教 案 第 3 教学周/星期 五 /第 7、8 节/第 5 次课 (3)假设没有因意外事故或自然灾害造成停产或减产的情况。 (4)假设r(t)为天然气产量的增长率。 (5)假设N(t)为天然气田的累积产量。 (6)假设为气田的年产量。 Q (7)假设气田的可采储量为N，相对应的开发时间t。 RR 3、问题分析及数学模型 根据所给的实际问题，预测气田的产量和可采储量，在这方面，目前国内外的方法很多，但各种 预测方法中有一种简单而实用的指数增长模型，它是借鉴英国人口学家Malthus(马尔萨斯)于1798 年提出的人口增长模型，而得到的。 若假设天然气产量的增长率为r(t)，它是时间t的连续函数。气田的累积产量设为N(t)，则 它们满足如下的关系: dN(t) ,r(t)N(t) dt dN(t)Q而气田的年产量，于是上述方程变为: Q,,r(t) dtN(t) Q有统计资料显示，气田的每年产量与累积产量之比与气田的开发时间t存在如下关系，即 N(t) Q (5) lg,A,Bt N(t) Q,btA或写成:，其中: ,a,ea,10,b,ln10,B,2.303B N(t) 假设气田的可采储量为N，相对应的开发时间为t，由上面的分析可得到方程: RR dN(t),,bt aeN(t),,, dt, ,N(t)N,,RR a,bt(,e)b解之可得 : (6) N(t),Ne R a,bt,e(),btb对上式求导预测模型为 (7) Q,a,N,e,e R 4 、 模型的分析与计算 (1)模型(3)中a , b的计算 QiA y,lg由于，所以关键问题在于求出A，B的值。由(1)式，设，a,10,b,ln10,B,2.303BiNi 其中为第i年的产量，N为第i年之前的累积产量，时间t以年为单位，则由所给数据可得 Qii 9 数 学 建 模 电 子 教 案 第 3 教学周/星期 五 /第 7、8 节/第 5 次课 y,A,Bti,1,2,？,20ii 根据线性回归的最小二乘估计，令 202 L(A,B),(y,A，Bt),iii,1 为使(，)最小，取分别关于，的偏导数，并令它们为0。 LABLAB 20,,L ,,,，,2(yABt)0,,ii,A,,i1 ,20,L ,,,，,2(yABt)t0,iii,,B ,,i1 s,ty ˆ,,B,s,tt 解此方程可得 ,202011 ,ˆˆA,y，tB(),,ii,2020 ii,,11, 20202020201122 其中: s,t,(y)s,ty,(t),(y),,,,,tyiiiittii2020ii,,ii,,i, 11111 ˆAˆ从而。 a,b,B10,2.303 (6)模型(7)中N的计算 R 对(2)式两边取常用对数可得 : (8) logN(t),,,,x a,bt其中:。由(8)式和所给的数据，建立回归方程(同上)，可求得α、,,logN,,,,x,e Rb , β，从而计算出油田的可采储量(略)。 N,10R (3)模型的求解 将根据上述求得的，，的值代入模型(6)(7)式便可计算出相应年份累积产量()和年abNNtR 产量Q的预测值。 用 mathematica气田储量的mathematica计算程序: 软件演示计%气田储量ch42 算程序 %文件名:ch42.m data1={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=Table[0,{n,20}]; i=2; For[data2[[1]]=data1[[1]],i<=20,i++,data2[[i]]=data2[[i-1]]+data1[[i]]] 10 数 学 建 模 电 子 教 案 第 3 教学周/星期 五 /第 7、8 节/第 5 次课 data3=Table[0,{m,20},{n,2}]; For[i=1,i<=20,i++,data3[[i]]={i,Log[10,data1[[i]]/data2[[i]]]}] Fit[data3,{1,t},t] aa=-0.0215995; bb=0.0809426; a=10^aa; b=Log[10.0]*bb; c=a/b; data4=Table[0,{m,20},{n,2}]; For[i=1,i<=20,i++,data4[[i]]={Exp[-b*i],Log[10,data2[[i]]]}] Fit[data4,{1,x},x] alpha=3.36832; bet=2.35678; Nr=10^alpha; Qp[t_]:=a*Nr*Exp[-c*Exp[-b*t]-b*t]; Np[t_]:=Nr*Exp[-c*Exp[-b*t]]; data5=Table[Qp[t],{t,1,20}]; data6=Table[Np[t],{t,1,20}]; compdata1=Table[0,{m,20},{n,2}]; For[i=1,i<=20,i++, compdata1[[i]]={data1[[i]],data5[[i]]}]; compdata2=Table[0,{m,20},{n,2}]; For[i=1,i<=20,i++, compdata2[[i]]={data2[[i]],data6[[i]]}]; f1=ListPlot[data1]; f2=Plot[Qp[t],{t,0,20}]; f3=ListPlot[data2]; f4=Plot[Np[t],{t,0,20}]; Show[f1,f2] Show[f3,f4] MatrixForm[compdata1] MatrixForm[compdata2] 执行后输出 A=-0.0215995 B=0.0809426 N=10^alpha alpha=3.36832 R 实际值与预测 值的图形曲线 很直观的看出 两者的接近程 度 11 数 学 建 模 电 子 教 案 第 3 教学周/星期 五 /第 7、8 节/第 5 次课 实际值与预测值对照表 8383 Q(10m/a) N(10m) T P年份 (a) 实际值 预测值 实际值 预测值 1957 1 19.0 26.674 19.0 1958 2 43.0 45.457 62.0 69.365 1959 3 59.0 68.603 121.0 126.117 1960 4 82.0 93.527 203.0 207.160 1961 5 92.0 117.186 295.0 312.745 1962 6 113.0 136.899 408.0 440.207 1963 7 138.0 150.897 546.0 584.626 1964 8 148.0 152.492 694.0 739.855 1965 9 151.0 159.935 845.0 899.552 1966 10 157.0 156.116 1002.0 1057.970 1967 11 158.0 148.242 1160.0 1210.430 1968 12 155.0 137.580 1315.0 1353.520 1969 13 137.0 125.282 1452.0 1485.050 1970 14 109.0 112.298 1561.0 1603.860 1971 15 89.0 99.351 1650.0 1709.660 1972 16 79.0 86.947 1729.0 1802.750 1973 17 70.0 75.409 1799.0 1883.840 1974 18 60.0 64.914 1859.0 1953.910 1975 19 53.0 55.534 1912.0 2014.040 1976 20 45.0 47.265 1957.0 2065.350 从结果看计算比较精确。 ?3.3 最优捕鱼策略 1 问题的提出: 这是1996年全国大学生数学建模竞赛的A题，问题如下: 为保护人类赖以生存的自然环境，可再生资源(如渔业、林业等资源)的开发必须适度。一种合 理、简化的策略是，在实现可持续收获的前提下，追求最大产量或最佳效益。 考虑对某种鱼的最优捕捞策略: 用多媒体进假设这种鱼分4个年龄组:称一龄鱼、二龄鱼、三龄鱼、四龄鱼。各年龄组每条鱼的平均重量分 行教学 别为5.07，11.55，17.86，22.99(克);各年龄组鱼的自然死亡率均为0.8(1/年);这种鱼为季节 5 性集中产卵繁殖，平均每条4龄鱼的产卵量为1.109×10(个)，3龄鱼的产卵量为这个数的一半，2 龄鱼和1龄鱼不产卵，产卵和孵化期为每年的最后4个月;卵孵化并成活为1龄鱼，成活率为(1龄 1111鱼条数与产卵总量n之比)1.22×10/(1.22×10+n)。 渔业管理部门，每年只允许在产卵孵化期的前8个月内进行捕捞作业。如果每年投入的捕捞 能力(如鱼船数、下网次数等)固定不变，这时单位时间捕捞量将与各年龄组鱼群条数成正比，比例 12 数 学 建 模 电 子 教 案 第 3 教学周/星期 五 /第 7、8 节/第 5 次课 系数不防称为捕捞强度系数。通常使用13mm网眼的拉网，这种网只能捕捞3龄鱼和4龄鱼，其两个 捕捞强度系数之比为0.42:1。渔业上称这种方式为固定努力量捕捞。 (1)建立数学模型分析如何实现可持续捕捞(即每年开始捕捞时渔场中各年龄组鱼群条数不变)， 并且在此前提下得到最高的年收获量(捕捞总重量)。 (2)某渔业公司承包这种鱼的捕捞业务5年，要求5年后鱼群的生产能力不能受到太大破 9坏。已知承包时各年龄组鱼群数量分别为:122，29.7，10.1，3.29(×10条)。如果仍用固定努力 量的捕捞方式，该公司采用怎样的策略才能使总收获量最高。 2、 模型假设及符号说明 (1)假设只考虑一种鱼的繁殖和捕捞，鱼群增长过程中不考虑鱼的迁入与迁出。 (2)假设各年龄组的鱼在一年内的任何时间都会发生自然死亡，产卵可在后四个月内任何时间发生。 (3)假设3、4龄鱼全部具有生殖能力，或者虽然雄鱼不产卵，但平均产卵量掩盖了这一差异。 (4)假设产卵期鱼的自然死亡率发生于产卵之后。 (5)假设各年龄组的鱼经过一年后，即进入高一级的年龄组，但4龄鱼经过一年后仍视为4龄鱼。 (6)假设对鱼的捕捞用固定努力量捕捞方式，每年的捕捞强度系数保持不变，且捕捞只在前八个月 进行。 (7)假设t时刻i龄鱼的数量为。 N(t)(i,1,2,3,4)i(k)(k)(8)假设第k年初i龄的数量为;第k年底i龄鱼的数量为(i=1，2，3，4)。 NNi0i1 (9)假设鱼的自然死亡率为r;4龄鱼的平均产卵量为c。 (10)假设第i龄鱼的平均重量为M(i=1，2，3，4)。 i k(11)假设第k年度鱼的产卵总量为。 Q (12)假设对第i龄鱼的捕捞强度系数为b;对i龄鱼的年捕捞量为a(i=1，2，3，4)。 ii (13)假设年总收获量为M，即M=Ma+Ma。 33445 (14)假设5年的总收获量为MM，即。 MM,M,i,i1 3 问题分析及数学模型 5由已知条件，可得 r,0.8;c,1.109，10 M,5.07;M,11.55;M,17.86;M,22.991234 ，(E为捕捞努力量)，r为自然死亡率，在[t，t+Δt]内，根据死亡 b,b,0;b,0.42E;b,E1234 率的定义，由于不捕捞1、2龄鱼，所以 介绍MalthusN(t),N(t，,t)dN(t)1iii r,lim,,,i,1,2(马尔萨斯),,0t,tN(t)N(t)dtii 人口增长模变形可得 型 dN(t) ,i,,rN(t)i,dt (1) i,1,2 , ,N(t),Nii0 t,0, ,rt解得。 N(t),Ne,i,1,20ii 对于3、4龄鱼由于捕捞在前8个月进行，因此，前8个月，捕捞与死亡均影响鱼的变化，因而 微分方程变形为 13 数 学 建 模 电 子 教 案 第 3 教学周/星期 五 /第 7、8 节/第 5 次课 dN(t)2,i,,[r，H(t,)b]N(t)ii, dt3 (2) ,i,3,4, , N(t),Nii0t,0, 其中 2, 0,t,,2,3 H(t,),, 23,1,t, ,3, 由(2)式解得 2,[r，H(t,)b]ti3 (3) N(t),Nei,3,4 ii0 (k)(k) 令为i龄鱼在k年底时的数量，为i龄鱼在k年初时的数量，得到1、2龄鱼k年末的NNi1i0 数量为 (k)(k),r,N,Ne,1110 (4) ,()()kk,r,N,Ne2120, 此处要用板对于3，4龄鱼，在第8个月末数量由(3)式可得 2书细讲 ,,(r，b)3(k)(k)3NNe,,,3031 , (4) ,2,(r，b)3 ,(k)(k)3NNe,,4041, 在后4个月，对于3、4龄鱼，只有死亡率起作用，因而微分方程为 dN(t),i ,,rN(t)i,dt i,3,4 ,(k)N(t)N,2,,1ii,t 3, ()k,rt解得: N(t),N,e,ii1 因而，3、4龄在k年末的数 212,(rb)rrb,，,,,33(k)(k)(k)333NNeeNe,,, ,313030,,4 () ,212 (rb)rrb,，,,,44,(k)(k)(k)333NNeeNe,,,414040, ,,将(4)式与(4)合在一起，得到在k年底各龄鱼的数量为 (k)(k),r,N,Ne1110, (k)(k),rN,Ne,2120 ,2 (5) ,r,b,3(k)(k)3 ,NNe,3130,2 ,,rb4(k)(k),3N,Ne4140, 由假设，到年底第i龄鱼全部转化为i+1龄鱼(i=1，2，3)，同时由卵孵化产生1龄鱼，于是 14 数 学 建 模 电 子 教 案 第 3 教学周/星期 五 /第 7、8 节/第 5 次课 得到 (k)(k,1), N,N2011,(k)(k,1) ,NN,3021,(k)(k,1)(k,1) (6) ,N,N，N403141, (k,1)d,Q,(k)N,10 (k,1),，dQ, (k,1)11其中，为年度总产卵量，且 Qd,1.22，10k,1 221r,b,b,c3411 (1)(1)(1)(1)(1)k,k,k,k,k,333 Q,(N，N)c,(e,N，e,N),e34304022 22,b,b34(k,1)(k,1) 33 (7) ,32529.55eN，65059.1e,N3040 此外，我们还求得每年对3、4龄鱼的总捕捞重量为 2 3M,(MbN(t)，MbN(t))dt 333444,0 22 ,(r，b),(r，b)MbMb34(k)(k)334433,[1,e]N，[1,e]N3040 r，br，b34 ,Ma，Ma3344 22(rb)(rb),，,， bb34(k)(k)3433 其中 a,[1,e]N,a,[1,e]N330440 r，br，b34 4 模型的分析与计算 (1)年度产量最优模型及其计算 (k)(k,1)为了实现可持续捕捞，即要求，在此前提下获得最高年收获量。结合基本模型，即N,Ni0i0 可得到年度产量最优模型。 此处用 maxM,max{Ma，Ma}3344MATLAB软件 ,演示计算程 kk,rkk,r()(),N,Ne,N,Ne序 11102120,22 ,,r,b,r,b34kkkk()()33N,Ne,N,Ne,31304140 ,(k)(k,1)(k)(k,1)N,N,N,N,20113021 ,k,(1)d,Qkk,k,k()(1)(1)(),N,N，N,N,40314110(k,1),d，Q ,22,b,b34,k,k,k,(1)(1)(1)33Q,32529.55Ne，65059.1Ne3040, (k)(k,1)(k)(k)(k)利用，也就是，，等与时间k无关，化简得 QN,NNNi0i1i0i0 15 数 学 建 模 电 子 教 案 第 3 教学周/星期 五 /第 7、8 节/第 5 次课 maxM,max{Ma，Ma}3344 d,QN,10d，Q 22,r，b,r，b(3)()3433N,eN，eN401040 22,r,b,b24433Q,32529.55e，65059.1eN40 下面给出MATLAB计算程序: %最优捕鱼策略ch331 %文件名:ch331.m x=sym('x'); b3=0.42*x; d=1.22*10^11; 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)))*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*10^11; 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+2/3*b3))/(1-exp(-(r+2/3*b4))) ; M3=17.86;M4=22.99; b3 b4 M=M3*a3+M4*a4; 16 数 学 建 模 电 子 教 案 第 3 教学周/星期 五 /第 7、8 节/第 5 次课 Max=M N1=N10 %各年龄组鱼的数量 N2=N1*exp(-r) N3=N2*exp(-r) N4=N40 执行后输出 b3 = 7.3359 b4 = 17.4664 Max =3.8886e+011 各年龄组数为 N1 = 1.1958e+011 N2 = 5.3730e+010 N3 = 2.4142e+010 N4 = 8.1544e+007 (2)承包期总产量模型及其计算 由于承包时鱼群的年龄组成尚未处于可持续捕捞的状态，因此不能简单地采用前面所得到的最优可持续捕捞的结果。问题中寻求收获最高的捕捞策略的条件是鱼群的生产力不能受到太大的破坏，我们可以把这个条件理解为捕捞五年后鱼群的年龄组成尽量接近于可持续捕捞鱼群的年龄结构。这时最优收获的问题就化为寻找一个收获的策略使得五年的鱼产量最高，而五年后鱼群的年龄构成尽量接近。 首先考虑无捕捞时鱼群的变化情况，在年度最优模型中，令，可以得到无捕捞时鱼群b,b,034 繁殖达到稳定状况的鱼群分布。此时，各年龄组鱼条数为 *11*10N,1.2199，10,N,5.4815，10,12 *10*10N,2.4630，10,N,2.0097，10,34 **其次，以为初值，用基本模型(5)~(7)式进行递推计算，经过N,N,N,N,N,N,0112234**六年左右的时间，鱼群即达到稳定状态。类似地，以为初值时进行递N,N,0,N,N,N,N123344推计算，可经过二年至三年的时间，鱼群即达到稳定状态。 我们认为要使生产能力不受到太大的破坏，只需保证5年后，即第6年初3、4龄鱼的数量保持在一定范围内即可。另外，由于生物学界认同的观点;只有种群减至低于最大持续产量水平时，资源才算开发过度。为此我们构造目标函数 (6)2(6)2 Y,[M(N,N)]，[M(N,N)]33443040 其中，N，N为可持续捕捞时的鱼群数量。 345 根据以上分析，令为五年总产量之和，我们可得承包期总产量模型为 MM,M,i,i1 17 数 学 建 模 电 子 教 案 第 3 教学周/星期 五 /第 7、8 节/第 5 次课 maxZ,max{MM,Y} , (k)k,r(k)k,r,N,Ne,N,Ne11102120,22,,r,b,r,b34kkkk()()33N,Ne,N,Ne,31304140 ,kk,kk,()(1)()(1)N,N,N,N,20113021 ,(k,1)d,Qkk,k,k()(1)(1)(),N,N，N,N,40314110k,(1),d，Q ,22,b,b34,k,k,k,(1)(1)(1)33,，Q32529.55Ne65059.1Ne3040, 用MATLAB计算结果为: b4 = 17.5772 MM = 1.6055e+012 Y = 5.9828e+009 z =5995e+012 第六年初各年龄组鱼的数量 N61 = 1.1941e+011 N62 = 5.3640e+010 N63 = 2.4127e+010 N64 = 7.7658e+007 每年的鱼捕捞量分别为 Mh1 = 1.0e+011 *2.3440 2.1485 3.9617 3.7783 3.8222 18 数 学 建 模 电 子 教 案 第 3 教学周/星期 五 /第 7、8 节/第 5 次课 19