unit14_4.doc

14.2 微分方程模型实验学习目标1. 会求放射性废物处理问题中的微分方程；2. 能求解阻滞增长人口模型中的微分方程；3. 会进行油气产量和可采储量预测问题中的计算；4. 能求解价格竞争问题和除雪机除雪问题中的微分方程。一、 放射性废物处理问题1. 问题 美国原子能委员会以往处理浓缩的放射性废料的方法，一直是把它们装入密封的圆桶里，然后扔到水深为90多米的海底。生态学家和科学家们表示担心，怕圆桶下沉到海底时与海底碰撞而发生破裂，从而造成核污染。原子能委员会分辨说这是不可能的。为此工程师们进行了碰撞实验，发现当圆桶下沉速度超过12.2m/s与海底相撞时，圆桶就可能发生碰裂。这样为避免圆桶碰裂，需要计算一下圆桶沉到海底时速度是多少？这时已知圆桶重量为239.46kg，体积为0.2058m3， 海水密度为1035.71 kg /m3。如果圆桶速度小于12.2m/s，就说明这种方法是安全可靠的，否则就要禁止用这种方法来处理放射性废料。假设水的阻力与速度大小成正比例，其正比例常数k=0.6。2. 问题建立的模型 圆桶的位移和速度分别满足下面的微分方程： （1） （2） 考虑受到的阻力，这时圆桶的速度应满足如下的微分方程： （3）3. 模型的解答 根据方程（1）和（2），加上初始条件v|t=0=0，s|t=0=0，以及题设的初始数据。通过Mathematica就可以求出圆桶的位移和速度的方程，具体步骤如下： m=239.46；v=0.2058；g=9.8；p=1035.71；k=0.6；ChopDSolvem st = = m g p g w k st， s0 = = 0，s0 = = 0， st，t DSolvem* vt = = m g p g w k vt， v0 = = 0，vt，t 得出位移的方程为： s（t）= - 171511+429.744t+171511e-0.00250564t （4） 速度的方程为： v（t）= 429.744 - 429.744e-0.00250564t （5） 通过方程（4）及s（t）= 90m，利用下面Mathematica程序： FindRoot90 = = - 171511+429.744t + 171511/Exp0.00250564t ， t，13求出时间t = 13.002s，再把它代入方程（5），求出圆桶的速度为v=13.7729m/s。 显然，此时圆桶的速度已超过12.2m/s，可以得出这种处理废料的方法不合理，美国原子能委员会已经禁止用这种方法来处理放射性废料。 根据方程（3），加上初始条件v| t=0 = 0，可利用下面的Mathematica程序：DSolvem vt= = m g-p g v k （vt）2， v0 = = 0，vt，t 求出圆桶的速度：v（t）=。 若把题设中F，k（仍设为0.6的话）和m的值代入上式，可得圆桶的速度为： v（t）=20.7303Tanh（0.0519426t） 这时若该速度要小于12.2m/s，那么利用Mathematica可得圆桶的运动时间就不能超过13s，位移不能超过84.8m。对应的Mathematica源程序如下： FindRoot20.7303*Tanh0.0519426t= =12.2， t，12 （*求时间*）Integrate20.7303*Tanh0.0519426t， t，0，13 （*求位移*）从这个模型，我们也可以得出原来处理核废料的方法是不合理的。二、 阻滞增长人口模型 假定人口数量p（t）满足微分方程p=ap-bp2（这是最为著名的Logistic人口模型），以下使用Mathematica解这个方程，源程序如下： 首先求出满足初始条件 p（t0）= p0 的解的一般表达式。 In1：=DSolvept = = a pt - b pt2，pt0 = = p0，pt，t Out1=p（t）= 再求人口增长速度最快的时期，即求使p(t)=0的时刻t，由于解出的表达式比较复杂，下面没有列出，但在Out4中给出该时刻的人口值。 In2：= p = pt / .% 1； In3：= SolveDp，t，2 = = 0，t； In4：= Simplifyp / . % 1； Out4= 以下代入具体数据，其中p0是1961年的世界人口数量，求出人口极限值大约是98.6亿，In7是绘制p（t）曲线的语句，得到如图14-5所示的图形。 In5：= a = 0.029；b = 2.941*10( - 4)；p0 = 30.6；t0 = 0； Limitp，t Out6= 98.6059。In7：= Plotp，t，-100，500，PlotRangeAll，AxesOrigin0，0图14-5 人口增长曲线三、 油气产量和可采储量的预测问题1 问题 根据某气田19571976年共20个年度的产气量数据（见下表），建立该气田的产量预测模型，并将预测值与实际值进行比较。年份1957195819591960196119621963产量108m31943598292113138年份1964196519661967196819691970产量108m3148151157158155137109年份197119721973197419751976产量108m38979706053452 问题分析与建立模型 我们将指数增长模型用于油气产量预测，并试着假设增长率r随时间t变化，即r是t的函数，从而得到油气田的累积产量Np与开发时间t的关系： 如果开发时间t以年为单位，则油气田的年产量Q=，方程可改写成 现在问题的关键是寻找油气产量的增长率r(t)了。1995年有人通过对国内外一些油气田开发资料的统计研究，得出结论：油气田的产量与累积产量之比（Q/Np），与其开发时间t存在着较好的半对数关系，即 Log=A-Bt 或写成 其中a = 10A，b = ln10B = 2.303B。 设油气田的可采储量为NR，相对应的开发时间为tR，由此，便得到预测油气产量的微分方程 这是一阶线性齐次微分方程，其解为 由于tR很大，。所以得到预测油气田累积产量的模型为 对上式求导，即得油气田年产量的预测模型为 为了确定油气田的可采储量N，我们对前一式两边取常用对数： 其中 3 计算过程 第一步：根据油气田实际生产数据，利用线性回归求得截距A和斜率B，进而计算出a、b之值。 第二步：计算出不同时间的x（=e-bt）和logNp，并进行Np与x的线性回归，求得截距a和斜率。 第三步：计算出油气田的可采储量NR =10。 第四步：将a、b和NR的值代入，即得预测油气田的累积产量和年产量的计算公式。 第五步：利用所得计算公式，计算相应年份累积产量Np和年产量Q的预测值。 Mathematica源程序如下： In1：=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； In2：=data2=Table0，n，20；i=2； Fordata21=data11，i=20，i+， data2i=data2i-1= data1iIn5：=data3=Table0，m，20，n，2；Fori=1，i=20，i+， data3i=i，Log10，data1i/data2i In6：=Fitdata3，1，t，t Out6= - 0.0215995-0.0809426t In7：= aa = - 0.0215995；bb=0.0809426； a=10aa；b=Log10.0*bb；c=a/b； In12：=data4=Table0，m，20，n，2； Fori=1，i=20，i+，data4i=Exp-b*i， Log10，data2i In13：=Fitdata4，1，x，x Out13= 3.36832 - 2.35678x In14：=alpha=3.36832；beta=2.35678；Nr=10alpha； In17：=Qpt_ ：=a*Nr*Exp-c*Exp-b*t-b*t； In18：=Npt_ ：=Nr*Exp-c*Exp-b*t； In19：=data5=TableQpt_ ，t，1，20； In20：=data6=TableNpt_ ，t，1，20；In21：=compdata1=Table0，m，20，n，2； Fori=1，i=20，i+，compdata1i= data1i，data5i；In22：=compdata2=Table0，m，20，n，2； Fori=1，i=20，i+，compdata2i= data2i，data6i； In23：=ListPlotdata1； In24：=PlotQpt，t，0，20； In25：=ListPlotdata2； In26：=PlotNpt，t，0，20； In27：=Show%23，%24 Out27：= Graphics（略） In28：=Show%25，%26 Out28：= Graphics（略） In29：=MatrixFormcompdata1 Out29：= （略） In30：=MatrixFormcompdata2 Out30：= （略）实际值与预测值对照表年份T(a)Q(108m3/a）Np(108m3)实际值预测值实际值预测值19571190266471901958243045457620693561959359068603121012611719604820935272030207160196159201171862950312745196261130136899408044020719637138015089754605846261964814801524926940739855196591510159935845089955219661015701561161002010579701967111580148242116001210430196812155013758013150135352019691313701252821452014850501970141090112298156101603860197115890993511650017096601972167908694717290180275019731770075409179901883840197418600649141859019539101975195305553419120201404019762045047265195702065350从上面图表中可以看出，预测结果是令人满意的。四、 价格竞争问题 位于同一条公路的甲、乙两个加油站彼此竞争激烈，同一城市中还有其它加油站。当乙站突然宣布降价后，甲站需要根据乙站的售价来调整自己的售价，以求获得尽可能高的利润。 下面引入一些指标： x 甲站的销售价格（便士/升）； y 乙站的销售价格（便士/升）； w 成本价（便士/升）； L 价格战前甲站的日销售量（升/日）； P 标准销售价格（便士/升）； e 甲站的利润（便士/日）。 问题化为已知y，求x使e最大。 设甲站的日销售量线性依赖于双方价格的变动，于是得到公式： In1：=L1=L - a（x-y）-b（p - y）+c（p - x）； e =（x - w）L1； 其中L1是甲站的日销售量，而e是利润。已知y时，e是x的函数，求e的最大值点就是求驻点： In3：=e1=（e）x Out3= L+c（p - x）+（-a - c）（- w + x）- b（p - y）- a（x - y） In4：= x = x/.Solvee1= =0，x 1 Out4= In5：=SimplifyeOut5= 其中Out4是驻点，Out5是利润的最大值。接下来代入具体数值进行计算：In6：= L= 20000；p = 40；w = 30； y = 37.0，38.0，39.0； a = 4000；b = c = 1000； x Out9= 35.5，36.，36.5 In10：= e Out10= 151250.，180000.，211250. In11：= y = 39.0；x = 35，37； e Out12= 200000.，210000. 这里一次给出3个y值，计算结果表明，当y = 37时，求出x = 35.5 ,总之必须低于对方的价格，结论合理。Out12表明，当y = 39时，求出的x = 36.5确实是最大值。五、 除雪机除雪问题 有条10km长的公路，由一台除雪机负责除雪。每当路面雪的平均厚度达到0.5m时，除雪机开始工作。但是雪仍在下，路面雪的厚度在不断增加，除雪机的前进速度会不断降低。当雪的厚度达到1.5m时，除雪机将无法工作。问除雪机能否将整条公路的积雪清除？当然这与降雪的速度有关，以下在一些合理的假设下进行讨论和计算。 已知：在无雪的路面上除雪机的行驶速度为10m/s；雪下了1时，雪最大时路面积雪的厚度以0.1cm/s的速度增加，前半小时雪越下越大，后半小时越下越小。 假定除雪机的速度v随雪的厚度h线性变化，利用已知条件可得v=10（1-2h/3）。而h是时间t的函数，假定前半小时h(t)匀速增加而后半小时匀速减少，可得 （单位：m/s）， 再积分得到： ，注意h（t）也是分段函数。 下面可以利用Mathematica进行计算了。In1：= h1 = 0.001（s /1800 UnitStep1800 - s + （2 s /1800）UnitSteps - 1800）； h = Integrateh1，s，0，t； v = 10（1 - 2h /3）； h / . t1800 Out4= 0.9 In5：= h / . t3600 Out5= 1.8 在In1中的h1是h(t)，对它求定积分得到雪的厚度h(t)。Out4给出t=1800 s时雪的厚度是0.9m，Out5表明t=3600 s时雪的厚度是1.8m。 除雪机从雪的厚度是0.5m时开始工作，直到雪的厚度是1.5m时停止，以下求出它开始和停止工作的时间，再积分得到它前进的距离。 In6：=FullSimplifyh，t1800 Out9= - 2.7777810-7 (-6145.58+t)(-1054.42+t) In10：=Solve%= =1.5，t Out10= t2560.77，t4639.23 In11：= t1 = t / . %1 Out11=2560.77 In12：= Integratev，t，t0，t1 Out12=3859.94 由于函数h(t)是分段的（表达式中含有函数UnitStep），解方程的函数Solve对它无能为力，只好分段求解。In6通过有条件的化简，成功地得到当t 1800时h(t)的表达式。接下来，使用Solve求除雪机开始工作的时刻t0，并从解集中提取出合理的答案。同样再得到除雪机停止工作的时刻t1，最后In12求除雪机前进的距离，Out12表明除雪机只能清除大约4km的积雪。 使用Solve求解不方便，如改用求数值解的函数FindRoot，则容易的多： In13：=FindRooth= =0.5，t，900 Out13=t13