机理分析建模概要

机理分析建模法———成都大学机理分析是根据对现实对象特性的认识，分析其因果关系，找出反映内部机理的规律。机理分析方法立足于揭示事物内在规律与问题相关的物理、化学、经济等方面的知识；通过对数据和现象的分析对事物内在规律做出的猜想(模型假设)。模型特点：有明确的物理或现实意义对现实对象的认识来源：机理分析建模常用方法：常微分方程偏微分方程逻辑方法比例方法代数方法目录常微分方程建模微分方程的建立微分方程的求解逻辑方法建模一微分方程建模当实际问题需寻求某个变量y随另一变量t的变化规律y=y(t)，且直接求很困难时，可以建立关于未知变量、未知变量的导数以及自变量的方程(即变量满足的微分方程)。在实际问题中，“改变”、“变化”、“增加”、“减少”等关键词提示我们注意什么量在变化；关键词“速率”、“增长”“衰变”，“边际的”，常涉及到导数。这些都是建立微分方程模型的关键。建立常微分方程模型的常用方法：运用已知物理定律利用平衡与增长式运用微元法运用分析法(一)微分方程的建立建立微分方程模型时应用已知物理定律，可事半功倍。例1.1一个较热的物体置于室温为180C的房间内，该物体最初的温度是600C，3分钟以后降到500C。想知道它的温度降到300C需要多少时间？10分钟以后它的温度是多少？牛顿冷却(加热)定律：将温度为T的物体放入处于常温m的介质中时，T的变化速率正比于T与周围介质的温度差。1、运用已知物理定律分析：假设房间足够大，放入温度较低或较高的物体时，室内温度基本不受影响，即室温分布均衡，保持为m，采用牛顿冷却定律是一个相当好的近似。建立模型：设物体在冷却过程中的温度为T(t)(t≥0)，“T的变化速率正比于T与周围介质的温度差”翻译成数学语言也就是：。成正比与mTdtdT60)0()(TmTkdtdT建立微分方程其中参数k0，m=18，求得一般解为ln(T－m)=－kt+c)0(tcemTkt或代入条件，求得c=42，,最后得2116ln31k)0(4218)(2116ln31tetTt该物体温度降至300C需要8.17分钟。结果：)(3.394218)10(0102116ln31CeT2、利用平衡与增长式许多研究对象在数量上常常表现出某种不变的特性，如封闭区域内的能量、货币量等。利用变量间的平衡与增长特性，可分析和建立有关变量间的相互关系.续人口增长模型对某地区时刻t的人口总数P(t)，除考虑个体的出生、死亡，再进一步考虑迁入与迁出的影响。在很短的时间段Δt内，关于P(t)变化的一个最简单的模型是：{Δt时间内的人口增长量}={Δt内出生人口数}－{Δt内死亡人口数}+{Δt内迁入人口数}－{Δt内迁出人口数}{Δt时间内的净改变量}={Δt时间内输入量}－{Δt时间内输出量}不同的输入、输出情况对应不同的差分或微分方程。更一般地输入量：含系统外部输入及系统内部产生的量；输出量：含流出系统及在系统内部消亡的量。此类建模方法的关键是分析并正确描述基本模型的右端，使平衡式成立。例1.2(战斗模型)两方军队交战，希望为这场战斗建立一个数学模型，应用这个模型达到如下目的：1.预测哪一方将获胜？2.估计获胜的一方最后剩下多少士兵？3.计算失败的一方开始时必须投入多少士兵才能赢得这场战斗？模型建立设:x(t)—t时刻X方存活的士兵数；y(t)—t时刻Y方存活的士兵数；假设：1)双方所有士兵不是战死就是活着参加战斗，x(t)与y(t)都是连续变量。2)Y方军队的一个士兵在单位时间内杀死X方军队a名士兵；3)X方军队的一个士兵在单位时间内杀死Y方军队b名士兵；{Δt时间内X军队减少的士兵数}={Δt时间内Y军队消灭对方的士兵数}平衡式：即有：Δx=－ayΔt，同理：Δy=－bxΔt令Δt0，得到微分方程组：)0()0(bbxdtdyaaydtdx基本思想：通过分析研究对象的有关变量在一个很短时间内的变化情况。3、微元法例1.3一个高为2米的球体容器里盛了一半的水，水从它的底部小孔流出，小孔的横截面积为1平方厘米。试求放空容器所需要的时间。2米对孔口的流速做两条假设：1．t时刻的流速v依赖于此刻容器内水的高度h(t)。2．整个放水过程无能量损失。分析：放空容器容器内水的体积为零容器内水的高度为零模型建立：由水力学知：水从孔口流出的流量Q为通过“孔口横截面的水的体积V对时间t的变化率”，即ghSdtdVQ262.0S—孔口横截面积(单位：平方厘米)h(t)—水面高度(单位：厘米)t—时间(单位：秒)当S=1平方厘米，有)1(262.0dtghdVh(t)h+Δh在[t，t+Δt]内，水面高度h(t)降至h+Δh(Δh0),容器中水的体积的改变量为ΔV=V(h)－V(h+Δh)=－πΔh[3(r12+r22)＋o(Δh)]≈－πr2Δh＋o(Δh)r1r2222200)100(100hhhr记令Δt0,得dV=－πr2dh，(2)比较(1)、(2)两式得微分方程如下：100)200(262.002thdhhhdtgh积分后整理得)31000700000(265.42523hhgt(0≤h≤100)令h=0，求得完全排空需要约2小时58分。基本思想：根据对现实对象特性的认识，分析其因果关系,找出反映内部机理的规律。例1.4(独家广告模型)广告是调整商品销售的强有力的手段,广告与销售量之间有什么内在联系？如何评价不同时期的广告效果？分析：广告的效果，可做如下的条件假设：1.商品的销售速度会因广告而增大，当商品在市场上趋于饱和时，销售速度将趋于一个极限值；4、分析法2.商品销售率(销售加速度)随商品销售速度的增高而降低；3.选择如下广告策略，t时刻的广告费用为：.,0;0,)(ttAtA建模：M—销售饱和水平，即销售速度的上限；λ(＞0)—衰减因子，广告作用随时间的推移而自然衰减的速度。直接建立微分方程记S(t)—t时刻商品的销售速度；)())(1)((tSMtStpAdtdS称p为响应系数，表征A(t)对S(t)的影响力。模型分析：是否与前三条假设相符？)())(()(tStSMMtApdtdS假设1市场余额假设2销售速度因广告作用增大,同时又受市场余额的限制。改写模型求解常微分方程模型的常用方法：微分方程的数值解微分方程的定性分析(二)微分方程的求解1.1常微分方程数值解的定义：在生产和科研中所处理的微分方程往往很复杂,且大多得不出一般解。而实际问题中对初值问题的求解，一般是要求得到在若干个点上满足规定精确度的近似值，或者得到一个满足精确度要求的便于计算的表达式。1、常微分方程的数值解0000121212'(,)()(),(),,(),,,.nnnyfxyxyxyxxxxxyxyxyxyyy即：对常微分方程，其数值解是指由初始点开始的若干离散的处的值，即对，求出准确值的相应近似值1.2建立数值解法的一些途径：用差商代替导数使用数值积分使用泰勒公式数值公式的精度100,0,1,2,,1,'(,)()iixxhinyfxyyxy设则可用以下离散化方法求解微分方程用差商代替导数（欧拉法）若步长h较小，则有hxyhxyxy)()()('故有公式：100(,)0,1,2,,-1()iiiiyyhfxyinyyx对方程y’=f(x,y),两边由xi到xi+1积分，并利用梯形公式，有：11111()()(,())[(,())(,())]2iixiixiiiiiiyxyxftytdtxxfxyxfxyx实际应用时，与欧拉公式结合使用：,2,1,0)],(),([2),()(11)1(1)0(1kyxfyxfhyyyxhfyykiiiiikiiiii(1)()1111120kkkiiiiiyyyyy（）对于已给的精确度，当满足时，取，然后继续下一步的计算．故有公式：)()],(),([200111xyyyxfyxfhyyiiiiii使用数值积分（改进的欧拉法）以此方法为基础，有龙格-库塔法、线性多步法等方法。使用泰勒公式数值公式的精度当一个数值公式的截断误差可表示为O（hk+1）（其中k为正整数，h为步长）时，称它是一个k阶公式．k越大，则数值公式的精度越高．•欧拉法是一阶公式，改进的欧拉法是二阶公式．•龙格-库塔法有二阶公式和四阶公式．•线性多步法有四阶亚当斯外插公式和内插公式．1.3用MATLAB软件求常微分方程的数值解[t，x]=solver(’f’,ts,x0,options)ode45ode23ode113ode15sode23s由待解方程写成的M文件名ts=[t0,tf]，t0、tf为自变量的初值和终值函数的初值ode23：组合的2/3阶龙格–库塔–费尔贝格算法ode45：运用组合的4/5阶龙格–库塔–费尔贝格算法自变量值函数值用于设定误差限(缺省时设定相对误差10-3,绝对误差10-6),命令为：options=odeset（’reltol’,rt,’abstol’,at）,rt，at：分别为设定的相对误差和绝对误差．1．在解含n个未知数的方程组时，x0和x均为n维向量，M文件中的待解方程组应以x的分量形式写出．2．使用MATLAB软件求数值解时，高阶微分方程必须等价地变换成一阶微分方程组．注意:0)0(,2)0(0)1(1000'222xxxdtdxxdtxd例：'121,yyxy解：令则微分方程变成一阶微分方程组：0)0(,2)0()1(1000211221'22'1yyyyyyyy1．建立M文件vdp1000．m如下：functiondy=vdp1000(t,y)dy=zeros(2,1);dy(1)=y(2);dy(2)=1000*(1-y(1)^2)*y(2)-y(1);2．取t0=0，tf=3000，输入命令：[T,Y]=ode15s('vdp1000',[03000],[20]);plot(T,Y(:,1),'-')3．结果如图：050010001500200025003000-2.5-2-1.5-1-0.500.511.52例1.5导弹追踪问题设位于坐标原点的甲舰向位于x轴上点A(1,0)处的乙舰发射导弹，导弹头始终对准乙舰．如果乙舰以最大的速度v0(常数)沿平行于y轴的直线行驶，导弹的速度是5v0，求导弹运行的曲线方程．乙舰行驶多远时，导弹将它击中？解法一（解析法）假设t时刻导弹的位置为P(x(t),y(t))，乙舰位于),1(0tvQ.由于导弹头始终对准乙舰，故此时直线PQ就是导弹的轨迹曲线弧OP在点P处的切线，即有xytvy1'0即yyxtv')1(0(1)又根据题意，弧OP的长度为AQ的5倍，即2001'd5xyxvt(2)由(1),(2)消去t,整理得模型:(3)'151)1(2yyx初值条件为:0)0(y0)0('y其解即为导弹的运行轨迹:245)1(125)1(855654xxy当1x时245y，即当乙舰航行到点)245,1(处时被导弹击中.被击中时间为:00245vvyt.若v0=1,则在t=0.21处被击中.1．建立M文件eq1．mfunctiondy=eq1(x,y)dy=zeros(2,1);dy(1)=y(2);dy(2)=1/5*sqrt(1+y(1)^2)/(1-x);2．取x0=0，xf=0．9999，建立主程序chase1．m如下：x0=0,xf=0.9999[x,y]=ode15s('eq1',[x0xf],[00]);plot(x,y(:,1),'b．')holdont=0:0.01:2;plot(1,t,'b*')2151''