数学建模第二章课件

1、第二章 微分方程本章学习目的：本章的主要目的在于：学习微分方程模型的建立、求解方法、分析结果及解决实际问题的全过程。1知道求解微分方程的解析法、数值解法以及图形表示解的方法；2理解求微分方程数值解的欧拉方法，了解龙格库塔方法的思想；3熟练掌握使用MATLAB软件的函数求微分方程的解析解、数值解和图形解；4通过范例学习怎样建立微分方程模型和分析问题的思想。2.1 引例在大学物理中，我们曾学习过简谐振动（如：弹簧振子、单摆），那就是一个典型的二阶常微分方程的模型。这里我们讨论“倒葫芦形状容器壁上的刻度问题”。x对于圆柱形状容器壁上的容积刻度，可以利用圆柱体体积公式：，其中容器的直径D为常数，体积V

2、与相对于容器底部的任意高度H成正比，因此在容器壁上可以方便地标出容积刻度。而对于几何形状不规则的容器，比如“倒葫芦形状”的容器壁上如何标出容积刻度呢？如图所示，建立坐标系，由微元法分析可知：，其中x表示高度，直径是高度的函数，记为D（x）。可得微分方程：如果该方程中的函数D(x)无解析表达式，只给出D(x)的部分测试数据，如何求解此微分方程呢？h=0.2;d=0.04,0.11,0.26,0.56,1.04,1.17;x(1)=0;v(1)=0;for k=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)

3、,v=v(1:6),plot(x,v) x = 0 0.2000 0.4000 0.6000 0.8000 1.0000v = 0 0.0011 0.0073 0.0373 0.1469 0.3393x = Columns 1 through 5 0 0.2000 0.4000 0.6000 0.8000 Column 6 1.0000v = Columns 1 through 5 0 0.0011 0.0073 0.0373 0.1469 Column 6 0.33932.2 微分方程模型的建立在工程实际问题中，“改变”、“变化”、“增加”、“减少”等关键词提示我们注意什么量在变化，关键词“

4、速率”、“增长”、“衰变”、“边际的”等常涉及到导数。我们熟悉的速度公式：就是一个简单的一阶微分方程。微分方程是指含有导数或微分的等式。一般形式：常用的建立微分方程的方法有：运用已知物理定律；利用平衡与增长式；运用微元法；应用分析法。2.2.1 运用已知物理定律建立微分方程模型时应用已知物理定律，可事半功倍。例2.1 一个较热的物体置于室温为180C的房间内，该物体最初的温度是600C，3分钟以后降到500C。想知道它的温度降到300C 需要多少时间？10分钟以后它的温度是多少？牛顿冷却（加热）定律：将温度为T的物体放入处于常温 m 的介质中时，T的变化速率正比于T与周围介质的温度差。分析：假

5、设房间足够大，放入温度较低或较高的物体时，室内温度基本不受影响，即室温分布均衡，保持为m，采用牛顿冷却定律是一个相当好的近似。建立模型：设物体在冷却过程中的温度为T(t)，t0， 根据牛顿加热（冷却）定律：，建立微分方程（2.1）其中参数k 0，m=18。2.2.2 利用平衡与增长式 许多研究对象在数量上常常表现出某种不变的特性，如封闭区域内的能量、货币量等。利用变量间的平衡与增长特性，可分析和建立有关变量间的相互关系。此类建模方法的关键是分析并正确描述基本模型的右端，使平衡式成立。例2.2 战斗模型：两方军队交战，希望为这场战斗建立一个数学模型，应用这个模型达到如下目的： 1. 预测哪一方将

6、获胜？ 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 =ayt同理 y =bxt 令，得到微分方程组： （2.2）2.2.3 微元法 基本思想：通过分析研究对象的有关变

7、量在一个很短时间内的变化情况建立微分方程。例2.3 一个高为2m的球体容器里盛了一半的水，水从它的底部小孔流出，小孔的横截面积为1cm2。试求放空容器所需要的时间。解：对孔口的流速做两条假设：（1）t 时刻的流速依赖于此刻容器内水的高度h(t). （2）整个放水过程无能量损失。 分析：放空容器意味着 模型建立：由流体力学知：水从孔口流出的流速Q为通过“孔口横截面的水的体积V对时间t 的变化率”，即（孔口流速公式） （2.3）S孔口横截面积（单位：cm2） h(t) 水面高度（单位：cm） t时间（单位：s）当S=1cm2，有。2.2.4 分析法基本思想：根据对现实对象特性的认识，分析其因果关系

8、，找出反映内部机理的规律。例2.4 独家广告模型 广告是调整商品销售的强有力的手段，广告与销售量之间有什么内在联系？如何评价不同时期的广告效果？解：1分析广告的效果，可做如下的条件假设：商品的销售速度会因广告而增大，当商品在市场上趋于饱和时，销售速度将趋于一个极限值；商品销售率（销售加速度）随商品销售速度的增高而降低。 2符号说明A(t) t 时刻的广告费用S(t) t 时刻商品的销售速度； M 销售饱和水平，即销售速度的上限； 衰减因子，广告作用随时间的推移而自然衰减的速度，0；p 响应系数，表征A(t) 对 S(t) 的影响力。3模型建立选择如下广告策略，t时刻的广告费用为： 建立微分方程

9、： （2.4）模型分析：是否与假设相符？2.3 微分方程求解方法2.3.1 解析解法解析解法只能解决一些特殊微分方程，这些方法主要针对：一阶特殊的微分方程：如使用分离变量法、方程变换法、线性方程的常数变易法或公式法求解。二阶或高阶常系数线性微分方程的特征根法。在高等数学的教程中有专门介绍。下面着重介绍微分方程的数值解法。232 数值解法微分方程的数值解法是解决某些实际问题中经常使用的方法。设待求解的定解问题为求该问题数值解法的基本过程如下：引入自变量取值点序列，定义 为步长，常用定步长(与n无关，为常数)，其精确解记为，一般难以得到。为了寻求的近似值，设想根据一定的原理，结合当前得到近似解，近

10、似地表示该点或前一点的导数值，由此推出计算的迭代公式。因此数值解法一般只能得到微分方程的近似解。下面介绍两个微分方程中最常用的数值解法。1.欧拉方法这是一种最简单的解微分方程的数值方法：就是在小区间xn, xn+1上用差商代替微商，可以得到近似的表达式若f(x,y)中的x取左端点，结合已经得到的y(xn)的近似值（数值解）yn，即，有y(xn+1)的近似值为 , n = 0,1,这就是求解微分方程的显式欧拉公式。也称向前欧拉公式。向前欧拉法计算简单，易于计算，但精度不高，收敛速度慢若f(x,y)中的x取右端点，可得向后欧拉公式如下：yn+1 = yn + h f (xn +1, yn +1)

11、，n = 0,1,称为隐式公式，因为要得出数值解yn+1，就必须求解这个非线性方程，计算比较困难。如果用将向前和向后欧拉公式加以平均，可得到梯形公式：该法的计算精度比向前和向后欧拉法都高，但计算和向后欧拉法一样困难。改进的欧拉算法：（1）先用向前欧拉法算出yn+1的预测值，（2）将预测值代入梯形公式的右端作为校正，得到yn+1，n=1,2,该式称为改进欧拉公式。例2.5 求解微分方程y = -y +x +1, y(0) = 1, 取步长h = 0.1和0.001。分别用三种数值解法求解，并结合其精确解，对求解误差进行分析比较。解 这是一个一阶线性微分方程，可用解析解法得到其精确解y = x +

12、 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)。x1(1)=0;y1(1)=1;y2(1)=1;y3(1)=1;h=0.1;for k=1:10x1(k+1)=

13、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

14、 = 0 0.1000 0.2000 0.3000 0.4000 0.5000 0.6000 0.7000 0.8000 0.9000 1.0000y = Columns 1 through 5 1.0000 1.0048 1.0187 1.0408 1.0703 Columns 6 through 10 1.1065 1.1488 1.1966 1.2493 1.3066 Column 11 1.3679y1 = 1.0000 1.0000 1.0100 1.0290 1.0561 1.0905 1.1314 1.1783 1.2305 1.2874 1.3487y2 = 1.0000 1.

15、0091 1.0264 1.0513 1.0830 1.1209 1.1645 1.2132 1.2665 1.3241 1.3855y3 = 1.0000 1.0048 1.0186 1.0406 1.0701 1.1063 1.1485 1.1963 1.2490 1.3063 1.3676图中，蓝色曲线是精确解，黑色曲线是向前欧拉法曲线，红色曲线是向后欧拉法曲线，绿色“* ”号为梯形法曲线。计算结果如下： 表2-1当h = 0.1时精确解向前欧拉法后退欧拉法梯形法011110.11.004811.00911.00480.21.01871.01001.02641.01860.31.0408

16、1.02901.05131.04060.41.07031.05611.08301.07010.51.10651.09051.12091.10630.61.14881.13141.16451.14850.71.19661.17831.21321.19630.81.24931.23051.26651.24900.91.30661.28741.32411.306311.36791.34871.38551.3676当h0.001时x1(1)=0;y1(1)=1;y2(1)=1;y3(1)=1;h=0.001;for k=1:1000x1(k+1)=x1(k)+h;y1(k+1)=(1-h)*y1(k)

17、+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);n=1;for k=1:11x1(k)=x1(n);y1(k)=y1(n);y2(k)=y2(n);y3(k)=y3(n);n=n+100;endx1=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 = Col

18、umns 1 through 5 0 0.1000 0.2000 0.3000 0.4000 Columns 6 through 10 0.5000 0.6000 0.7000 0.8000 0.9000 Column 11 1.0000y = Columns 1 through 5 1.0000 1.0048 1.0187 1.0408 1.0703 Columns 6 through 10 1.1065 1.1488 1.1966 1.2493 1.3066 Column 11 1.3679y1 = Columns 1 through 5 1.0000 1.0048 1.0186 1.04

19、07 1.0702 Columns 6 through 10 1.1064 1.1486 1.1964 1.2491 1.3064 Column 11 1.3677y2 = Columns 1 through 5 1.0000 1.0049 1.0188 1.0409 1.0705 Columns 6 through 10 1.1067 1.1490 1.1968 1.2495 1.3068 Column 11 1.3681y3 = Columns 1 through 5 1.0000 1.0048 1.0187 1.0408 1.0703 Columns 6 through 10 1.106

20、5 1.1488 1.1966 1.2493 1.3066 Column 11 1.3679表2-2 当h = 0.001时 精确解向前欧拉法后退欧拉法梯形法011110.11.00481.00481.00491.00480.21.01871.01861.01881.01870.31.04081.04071.04091.04080.41.07031.07021.07051.07030.51.10651.10641.10671.10650.61.14881.14861.14901.14880.71.19661.19641.19681.19660.81.24931.24911.24951.2493

21、0.91.30661.30641.30681.306611.36791.36771.36811.3679计算结果表明：当步长h =0.1时，它们的前两位有效数字是精确的；当步长h =0.001时，它们的前四位有效数字是精确的。说明在迭代中，步长h越小，计算结果越精确。通过进一步的计算，还可以发现：迭代离开初始点越远，误差越大。如上图所示，中间的曲线表示精确解曲线。2 误差和阶衡量求解公式好坏的一个主要标准是求解公式的精度，因此引入局部截断误差和阶数的概念。定义2.1 假定没有误差，即，用数值方法计算的误差，称为该数值方法计算时的局部截断误差。为估计欧拉公式的局部截断误差，先将精确解在处作泰勒级

22、数展开（2.12）对于向前欧拉公式，在的假定下可记作（2.13）由式（2.12）和（2.13）之差可得到式中h的最低阶项称为局部截断误差主项，它对于h是2阶的。同理，对于向后欧拉公式的局部截断误差是，梯形公式和改进欧拉公式的局部截断误差是。求解微分方程的数值计算方法的精度是由局部截断误差中h的阶定义的：如果一个方法的局部截断误差为，则称该方法具有p阶精度。因此向前和向后欧拉方法的精度为1阶，梯形公式和改进欧拉公式的精度为2阶。3. 龙格-库塔方法数值方法具有越高阶的精度，得到的解就越精确。我们发现，向前和向后欧拉公式各用了区间一个端点的导数，具有1阶精度，而梯形和改进欧拉公式用了两个端点的导数

23、取平均，得到了2阶精度。因此就引导人们利用内的若干点的导数，对它们作线性组合得到平均斜率，由此得到更高阶的精度，这就是龙格-库塔方法的基本思路。Runge-Kutta方法的导出对于常微分方程的初值问题的解，在区间上使用微分中值定理，有，其中。即设则，即（2.5）k可以认为是在区间上的平均斜率，只要使用不同方法给出在区间内平均斜率的近似值k，就可得到不同的数值计算方法。如：为向前欧拉方法；为向后欧拉方法。1）. 二阶R-K方法一般形式为（2.15）若满足，式（2.15）具有2阶的精度。如果时，就得到改进的欧拉方法。2）. 三阶R-K方法（2.16）式（2.16）具有3阶的精度。3）. 四阶经典R

24、-K方法（2.17）式（2.17）具有4阶精度。在MATLAB软件中含有数值求解的系统函数，其实现原理就是龙格库塔方法，MATLAB求解方法将在2.4节中介绍。2.4 Matlab软件求解前面已经介绍了微分方程（组）的数值解法。这些方法计算工作量大，需要借助于计算机实现。下面简单介绍Matlab在此领域的应用。2.4.1 解析解用MATLAB命令dsolve(eqn1,eqn2, .) 求常微分方程（组）的解析解。其中eqni表示第i个微分方程，Dny表示y的n阶导数,默认的自变量为t。1. 微分方程例2.6 求解一阶微分方程 (1) 求通解输入：dsolve(Dy=1+y2) 输出：ans

25、=tan(t+C1) （2）求特解输入：dsolve(Dy=1+y2,y(0)=1,x) 指定初值为1，自变量为x输出：ans =tan(x+1/4*pi) 例2.7 求解二阶微分方程 输入:dsolve(D2y+(1/x)*Dy+(1-(1/2)2/x2)*y=0,y(pi/2)=2,Dy(pi/2)=-2/pi,x) ans =2(1/2)*pi(1/2)/x(1/2)*sin(x) 化简输出结果，输入：pretty(ans) 1/2 1/2 2 pi sin(x) - 1/2 x 即：2微分方程组例2.8 求解 df/dx=3f+4g; dg/dx=-4f+3g。（1）通解:f,g=ds

26、olve(Df=3*f+4*g,Dg=-4*f+3*g) f =exp(3*t)*(C1*sin(4*t)+C2*cos(4*t)g =exp(3*t)*(C1*cos(4*t)-C2*sin(4*t) 特解：f,g=dsolve(Df=3*f+4*g,Dg=-4*f+3*g,f(0)=0,g(0)=1) f =exp(3*t)*sin(4*t)g =exp(3*t)*cos(4*t) 2.4.2 数值解在微分方程(组)难以获得解析解的情况下，可以用Matlab方便地求出数值解。格式为:t,y = ode23(F,ts,y0,options,p1,p2,.)注意： 微分方程的形式：y = F(

27、t, y)，t为自变量，y为因变量（可以是多个，如微分方程组）； t, y为输出矩阵，分别表示自变量和因变量的取值； F代表微分方程组的函数名（m文件，必须返回一个列向量）； ts的取法有几种，（1）ts=t0, tf 表示自变量的取值范围，（2）ts=t0,t1,t2,tf,则输出在指定时刻t0,t1,t2,tf处给出，（3）ts=t0:k:tf,则输出在区间t0,tf的等分点给出； y0为初值条件； options用于设定误差限（缺省是设定相对误差是10(-3)，绝对误差是10(-6)）； p1,p2,.用于传递附加的参数值。ode23是微分方程组数值解的低阶方法，ode45为中阶方法，与

28、ode23类似。比如例2.7的数值解：解：令y1=y,y2=y1,将二阶微分方程转化为一阶微分方程组y1=y2y2=-y2/x+(n/x)2-1)y1首先建立M-文件函数: function f=jie3(x,y) f=y(2);-y(2)/x+(1/2)2/x2-1)*y(1);计算:x,y=ode23(jie3,pi/2,pi,2,-2/pi) x = 1.5708 1.6074 1.7645 1.9215 2.0786 2.2357 2.3928 2.5499 2.7069 2.8640 3.0211 3.1416y = 2.0000 -0.6366 1.9758 -0.6869 1.8

29、518 -0.8879 1.6982 -1.0631 1.5192 -1.2108 1.3193 -1.3293 1.1032 -1.4174 0.8756 -1.4744 0.6416 -1.5002 0.4060 -1.4951 0.1735 -1.4602 0.0002 -1.4140 作图程序：y1=y(:,1);y2=y(:,2);plot(x,y1,x,y2,r),gtext(y1),gtext(y2) ? Error using = gtextInterrupted命令gtext（）在MATLAB下执行，则不会有错误提示，会在图形窗口出现十字线，其交点是括号内字符串的位置，移动鼠

30、标可移动该交点，鼠标单击一下就可将字符串固定在那里。例2.9 求解一个经典的范得波（Van Der pol）微分方程：解 形式转化：令。则范得波方程转化为一阶微分方程组： 。编写M文件如下（必须是M文件表示微分方程组，一般地，M文件的名字与函数名相同）： function dot1=vdpol(t,y); dot1=y(2); (1-y(1)2)*y(2)-y(1);在命令窗口写如下命令:t,y=ode23(vdpol,0,20,1,0),y1=y(:,1);y2=y(:,2);plot(t,y1,t,y2,-);title(Van Der Pol Solution );xlabel(Time

31、,Second);ylabel(y(1)andy(2) 执行：注：Van der Pol方程描述具有一个非线性振动项的振动子的运动过程。最初,由于它在非线性电路上的应用而引起广泛兴趣。一般形式为。243 图形解无论是解析解还是数值解，都不如图形解直观明了。即使是在得到了解析解或数值解的情况下，作出解的图形，仍然是一件深受欢迎的事。这些都可以用Matlab方便地进行。1.图示解析解如果微分方程（组）的解析解为：y=f (x)，则可以用Matlab函数fplot作出其图形：fplot(fun,lims)其中：fun给出函数表达式；lims=xmin xmax ymin ymax限定坐标轴的大小。例

32、如fplot(sin(1/x), 0.01 0.1 -1 1) 2.图示数值解设想已经得到微分方程（组）的数值解(x,y)。可以用Matlab函数plot(x,y)直接作出图形。其中x和y为向量（或矩阵）。2.5 范例：地中海鲨鱼问题意大利生物学家Ancona曾致力于鱼类种群相互制约关系的研究，他从第一次世界大战期间（1914年7月28日1918年11月11日）,地中海各港口捕获的几种鱼类捕获量百分比的资料中，发现鲨鱼的比例有明显增加（见下表）。年代19141915191619171918百分比11.921.422.121.236.4年代19191920192119221923百分比27.31

33、6.015.914.819.7战争为什么使鲨鱼数量增加？是什么原因？因为战争使捕鱼量下降，食用鱼增加，显然鲨鱼也随之增加。 但为何鲨鱼的比例大幅增加呢？生物学家Ancona无法解释这个现象，于是求助于著名的意大利数学家V.Volterra，希望建立一个食饵捕食者系统的数学模型，定量地回答这个问题。 1、符号说明：x1(t), x2(t)分别是食饵、捕食者（鲨鱼）在t时刻的数量； r1, r2是食饵、捕食者的固有增长率；1是捕食者掠取食饵的能力， 2是食饵对捕食者的供养能力；2、基本假设： 捕食者的存在使食饵的增长率降低，假设降低的程度与捕食者数量成正比，即食饵对捕食者的数量x2起到增长的作用，

34、 其程度与食饵数量x1成正比，即综合以上和，得到如下模型：模型一：不考虑人工捕获的情况 该模型反映了在没有人工捕获的自然环境中食饵与捕食者之间的制约关系，没有考虑食饵和捕食者自身的阻滞作用，是Volterra提出的最简单的模型。给定一组具体数据，用matlab软件求解。 食饵： r1= 1, 1= 0.1, x10= 25; 捕食者(鲨鱼)：r2=0.5, 2=0.02, x20= 2;编制程序如下1、建立m-文件shier.m如下： function dx=shier(t,x) dx=zeros(2,1); dx(1)=x(1)*(1-0.1*x(2); dx(2)=x(2)*(-0.5+0

35、.02*x(1);2、在命令窗口执行如下程序： t,x=ode45(shier,0:0.1:15,25 2); plot(t,x(:,1),-,t,x(:,2),*),grid on 图中，蓝色曲线和绿色曲线分别是食饵和鲨鱼数量随时间的变化情况，从图中可以看出它们的数量都呈现出周期性，而且鲨鱼数量的高峰期稍滞后于食饵数量的高峰期。模型（二） 考虑人工捕获的情况假设人工捕获能力系数为e，相当于食饵的自然增长率由r1 降为r1-e，捕食者的死亡率由r2 增为 r2+e，因此模型(一)修改为：设战前捕获能力系数e=0.3, 战争中降为e=0.1, 其它参数与模型(一)的参数相同。观察结果会如何变化？1）当e=0.3时：2）当e=0.1时：分别求出两种情况下鲨鱼在鱼类中所占的比例。即计算画曲线：plot(t,p1(t),t,p2(t)