数学建模第二章课件

1第二章微分方程本章学习目的：本章的主要目的在于：学习微分方程模型的建立、求解方法、分析结果及解决实际问题的全过程。1．知道求解微分方程的解析法、数值解法以及图形表示解的方法；2．理解求微分方程数值解的欧拉方法，了解龙格——库塔方法的思想；3．熟练掌握使用MATLAB软件的函数求微分方程的解析解、数值解和图形解；4．通过范例学习怎样建立微分方程模型和分析问题的思想。§2.1引例在《大学物理》中，我们曾学习过简谐振动（如：弹簧振子、单摆）0222xdtxd，那就是一个典型的二阶常微分方程的模型。这里我们讨论“倒葫芦形状容器壁上的刻度问题”。2对于圆柱形状容器壁上的容积刻度，可以利用圆柱体体积公式：4/2HDV，其中容器的直径D为常数，体积V与相对于容器底部的任意高度H成正比，因此在容器壁上可以方便地标出容积刻度。而对于几何形状不规则的容器，比如“倒葫芦形状”的容器壁上如何标出容积刻度呢？如图所示，建立坐标系，由微元法分析可知：dxxDdV2)(41，其中x表示高度，直径是高度的函数，记为D（x）。可得微分方程：0)0()(412VxDdxdV如果该方程中的函数D(x)无解析表达式，只给出D(x)的部分测试数据，如何求解此微分方程呢？h=0.2;d=[0.04,0.11,0.26,0.56,1.04,1.17];x3x(1)=0;v(1)=0;fork=1:5x(k+1)=x(k)+h;v(k+1)=v(k)+(h/2)*(pi/4)*(d(k)^2+d(k+1)^2);endx=x(1:6),v=v(1:6),plot(x,v)x=00.20000.40000.60000.80001.0000v=00.00110.00730.03730.14690.3393x=Columns1through500.20000.40000.60000.8000Column61.0000v=Columns1through500.00110.00730.03730.1469Column60.3393400.10.20.30.40.50.60.70.80.9100.050.10.150.20.250.30.35§2.2微分方程模型的建立在工程实际问题中，“改变”、“变化”、“增加”、“减少”等关键词提示我们注意什么量在变化，关键词“速率”、“增长”、“衰变”、“边际的”等常涉及到导数。我们熟悉的速度公式：vdtdy就是一个简单的一阶微分方程。微分方程是指含有导数或微分的等式。一般形式：).,,,,(0),,,,()1()()(nnnyyyxfyyyyxF或：5常用的建立微分方程的方法有：运用已知物理定律；利用平衡与增长式；运用微元法；应用分析法。2.2.1运用已知物理定律建立微分方程模型时应用已知物理定律，可事半功倍。例2.1一个较热的物体置于室温为180C的房间内，该物体最初的温度是600C，3分钟以后降到500C。想知道它的温度降到300C需要多少时间？10分钟以后它的温度是多少？牛顿冷却（加热）定律：将温度为T的物体放入处于常温m的介质中时，T的变化速率正比于T与周围介质的温度差。分析：假设房间足够大，放入温度较低或较高的物体时，室内温度基本不受影响，即室温分布均衡，保持为m，采用牛顿冷却定律是一个相当好的近似。建立模型：设物体在冷却过程中的温度为T(t)，t≥0，根据牛顿加热（冷却）定律：）成正比与（mTdtdT，建立微分方程60)0()(TmTkdtdT（2.1）其中参数k0，m=18。2.2.2利用平衡与增长式许多研究对象在数量上常常表现出某种不变的特性，如封闭区域内的能量、货币量等。利用变量间的平衡与增长特性，可分析和建立有关变量间的相互关系。6此类建模方法的关键是分析并正确描述基本模型的右端，使平衡式成立。例2.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令0t，得到微分方程组：)0()0(bbxdtdyaaydtdx（2.2）72.2.3微元法基本思想：通过分析研究对象的有关变量在一个很短时间内的变化情况建立微分方程。例2.3一个高为2m的球体容器里盛了一半的水，水从它的底部小孔流出，小孔的横截面积为1cm2。试求放空容器所需要的时间。解：对孔口的流速做两条假设：（1）t时刻的流速依赖于此刻容器内水的高度h(t).（2）整个放水过程无能量损失。分析：放空容器意味着容器内水的高度为零容器内水的体积为零模型建立：由流体力学知：水从孔口流出的流速Q为通过“孔口横截面的水的体积V对时间t的变化率”，即ghSdtdVQ262.0（孔口流速公式）（2.3）S—孔口横截面积（单位：cm2）h(t)—水面高度（单位：cm）t—时间（单位：s）当S=1cm2，有dtghdV262.0。2.2.4分析法基本思想：根据对现实对象特性的认识，分析其因果关系，找出反映内部机理的规律。8例2.4独家广告模型广告是调整商品销售的强有力的手段，广告与销售量之间有什么内在联系？如何评价不同时期的广告效果？解：1．分析广告的效果，可做如下的条件假设：①商品的销售速度会因广告而增大，当商品在市场上趋于饱和时，销售速度将趋于一个极限值；②商品销售率（销售加速度）随商品销售速度的增高而降低。2．符号说明A(t)—t时刻的广告费用S(t)—t时刻商品的销售速度；M—销售饱和水平，即销售速度的上限；λ—衰减因子，广告作用随时间的推移而自然衰减的速度，λ＞0；p—响应系数，表征A(t)对S(t)的影响力。3．模型建立选择如下广告策略，t时刻的广告费用为：ttAtA,00,)(建立微分方程：)())(1)((tSMtStpAdtdS（2.4）模型分析：是否与假设相符？§2.3微分方程求解方法2.3.1解析解法解析解法只能解决一些特殊微分方程，这些方法主要针对：9一阶特殊的微分方程：如使用分离变量法、方程变换法、线性方程的常数变易法或公式法求解。二阶或高阶常系数线性微分方程的特征根法。在高等数学的教程中有专门介绍。下面着重介绍微分方程的数值解法。2．3．2数值解法微分方程的数值解法是解决某些实际问题中经常使用的方法。设待求解的定解问题为00)(),,(yxyyxfdxdy求该问题数值解法的基本过程如下：引入自变量取值点序列nx，定义1nnnxxh为步长，常用定步长(nh与n无关，为常数)，其精确解记为)(nxy，一般难以得到。为了寻求)(nxy的近似值ny，设想根据一定的原理，结合当前得到近似解，近似地表示该点或前一点的导数值，由此推出计算ny的迭代公式。因此数值解法一般只能得到微分方程的近似解ny。下面介绍两个微分方程中最常用的数值解法。1.欧拉方法10这是一种最简单的解微分方程的数值方法：就是在小区间[xn,xn+1]上用差商代替微商，可以得到近似的表达式),()()(1yxfhxyxynn若f(x,y)中的x取左端点nx，结合已经得到的y(xn)的近似值（数值解）yn，即nnyxy)(，有y(xn+1)的近似值为),(1nnnnyxfhyy,n=0,1,…这就是求解微分方程的显式欧拉公式。也称向前欧拉公式。向前欧拉法计算简单，易于计算，但精度不高，收敛速度慢若f(x,y)中的x取右端点1nx，可得向后欧拉公式如下：yn+1=yn+hf(xn+1,yn+1)，n=0,1,…称为隐式公式，因为要得出数值解yn+1，就必须求解这个非线性方程，计算比较困难。如果用将向前和向后欧拉公式加以平均，可得到梯形公式：)],(),([21111nnnnnnyxfyxfhyy该法的计算精度比向前和向后欧拉法都高，但计算和向后欧拉法一样困难。改进的欧拉算法：（1）先用向前欧拉法算出yn+1的预测值1ny，11),(1nnnnyxfhyy（2）将预测值1ny代入梯形公式的右端作为校正，得到yn+1)],(),([21111nnnnnnyxfyxfhyy，n=1,2,…该式称为改进欧拉公式。例2.5求解微分方程y'=-y+x+1,y(0)=1,取步长h=0.1和0.001。分别用三种数值解法求解，并结合其精确解，对求解误差进行分析比较。解这是一个一阶线性微分方程，可用解析解法得到其精确解y=x+e-x。三种数值解如下：（h=0.1）1)向前欧拉法：迭代公式为yn+1=（1-h）yn+hxn+h，n=0,1,...。其中y0=y(0)=1。2)后退欧拉法：由后退欧拉法隐式公式得yn+1=yn+0.1(-yn+1+xn+1+1)，变形为yn+1=(yn+hxn+1+h)/(1+h)。3)梯形法：将隐式梯形公式转化为显示迭代公式如下：yn+1=(yn+(h/2)*(-yn+xn+xn=1+2))/(1+h/2)。12x1(1)=0;y1(1)=1;y2(1)=1;y3(1)=1;h=0.1;fork=1:10x1(k+1)=x1(k)+h;y1(k+1)=(1-h)*y1(k)+h*x1(k)+h;y2(k+1)=(y2(k)+h*x1(k+1)+h)/(1+h);y3(k+1)=(y3(k)+(h/2)*(-y3(k)+x1(k)+x1(k+1)+2))/(1+h/2);endx=0:0.1:1;y=x+exp(-x);x1=x1(1:11),y=y(1:11),y1=y1(1:11),y2=y2(1:11),y3=y3(1:11),plot(x,y,x1,y1,'k:',x1,y2,'r--',x1,y3,'g*')程序中，x1为自变量，y为精确解，y1、y2、y3分别为向前欧拉法、后退欧拉法和梯形法的解。结果如下：x1=00.10000.20000.30000.40000.50000.60000.70000.80000.90001.0000y=Columns1through51.00001.00481.01871.04081.0703Columns6through101.10651.14881.19661.24931.306613Column111.3679y1=1.00001.00001.01001.02901.05611.09051.13141.17831.23051.28741.3487y2=1.00001.00911.02641.05131.08301.12091.16451.21321.26651.32411.3855y3=1.00001.00481.01861.04061.07011.10631.14851.19631.24901.3063141.367600.10.20.30.40.50.60.70.80.9111.051.11.151.21.251.31.351.4图中，蓝色曲线是精确解，黑色曲线是向前欧拉法曲线，红色曲线是向后欧拉法曲线，绿色“*”号为梯形法曲线。计算结果如下：表2-1当h=0.1时nx精确解向前欧拉法后退欧拉法梯形法011110.11.004811.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.