数学实验第三章.ppt

1、第三章 微分方程,1. 微分方程概述,2. 用数值法解微分方程-MATLAB解微分方程的原理,3. 图形法解微分方程,4. MATLAB解微分方程,5. 微分方程建模及求解实例,3.1 导 言,在许多实际问题中，量x与量y之间的函数关系，常常以y(x)的导数y(x)形式表现出来.生物种群以及人口数量y(x)中用到是“速率y(x)”、“增长率”，在放射性问题中用到的是“衰变”,在经济学中用到的是“边际”等，因此在求解y(x)的方程中通常包含了y(x)方面的信息. 微分方程是求解未知函数的一类方程.解微分方程是认识函数关系的重要方法.微分方程是研究变化规律的有利工具，在科技、工程、经济管理、生态、

2、环境、人口等各个领域中有着广泛的应用。,微分方程是由未知函数y(x),导数y(n)(x),自变量x构成的方程 例如: y(x)=y(x) y(x)*x=1 y(2)(x)+y(x)+y(x)+x+1=0 y(x)=snnx+x2+1,y(x)=exp(x),y(x)=lnx,要建立一个微分方程就是寻找函数、函数对变量的导数和变量之间的一种平衡关系。而建立这种平衡关系,是需要根据物理的、非物理的原理或规律、作出适当的假设。,建立微分方程只是解决问题的第一步,通常还需要求出微分方程的解,进一步,要把解放在实际现象中进行解释和加以检验. 如果能得到微分方程的解析解固然便于对问题的分析和应用,但得到解

3、析解往往只局限于一些特殊类型的微分方程,如常系数的或特殊函数形式的微分方程等,而绝大多数的微分方程等都是求不出解析解的,如变系数的、非线性函数形式的微分方程等，这就需要借助于数值解法。数值解法是求解微分方程的一个十分重要的手段。,本章主要介绍几个有关微分方程的实例，微分方程（组）的数值解法和 图解法，用MATLAB软件求解微分方程的解析解和数值解以及用图形的形 式表示解曲线，对解的性态进行分析和检验，最后通过范例展现解决微分 方程实际问题的全过程。,认 识 动 态 过 程 的 一 般 规 律,数 值 解 法 的 地 位,NEXT,3.2 引例：单摆运动,这是一个我们熟悉的物理模型，图3.1中一

4、根（无弹性）细线，一端固定，另一端悬挂一个质量为m的小球，在重力作用下小球处于竖直的平衡位置。若使小球偏离平衡位置一个角度 ，它会沿圆弧摆动。在不考虑空气阻力的情况下，小球作一定周期的简谐运动。研究函数(t),m,小球所受的重力为mg 小球在圆周切线方向上的受力为mgsnn t时刻小球的角速度为(t) t时刻小球的角加速度为(t) t时刻小球的线加速度为l*(t) t时刻小球在圆周切线方向上的受力为ml*(t),mg,1）该微分方程是线性的还是非线性的？是否存在解析解？ 2）如果不存在解析解，能否求出其近似解？,想,ml (t) =mgsnn(t) (0)=0 (0)=0,因此,利用牛顿第二定

5、律,可得到如下的微分方程,这是一个通过力学分析,建立微分方程的实例,注意:,我们的方程是 l (t) =gsnn(t) 不是 l (t) =gsnn(t),另一个实例:倒葫芦形状容器壁上的刻度问题 对于圆柱形状容器壁上的容积刻度,可以利用圆柱体体积公式： V=D2h/4其中直径D为常数。由于体积V与(相对容器底部的)任意高度h 的函数关系明确，因此在容器壁 上可以方便地标出容积刻度。而对于几何形状不是规则的容器，比如倒葫芦形状的容器壁上如何标出容积刻度呢？下面是某特殊容器经测量得到的几个部分容器高度h与直径D的关系,表3.1,h 0 0.2 0.4 0.6 0.8 1.0 D 0.04 0.1

6、1 0.26 0.56 1.04 1.17,根据表3.1中的数据，可以拟合出倒葫芦形状容器图3.2,h,用微元思想分析:在高度 x 处,给出高度微量x,使这一段容器近似于圆柱体,因此有 V=1/4D(x)2x , V(0)=0 两边同除以x ，得 V/x =1/4D(x)2 , V(0)=0 两边对x 取极限，得微分方程 dV/dx =1/4D(x)2 , V(0)=0 其中x表示高度，D(x)是高度为x时的容器直径，表格3.1给出了6组测量数据,(3.1),1）此微分方程中含有表格表示的离散函数，是否存在解析解？ 2）如果不存在解析解，如何求解或求近似解？通过数值解可以 得到近似程度较高的V

7、(x).,想,x,o,x,图3.2 倒葫芦形状容器,如图3.2所示，建立坐标系，根据测量数据，得,o,一般形式为 f (x,y(x),y(x),.,y (n) (x)=0,或,微分方程在解决工程技术、社会和经济问题中有着十分广泛的应用。在解决实际问题时，必须经过两个重要的阶段，一是微分方程的建立，它取决于对实际问题的深入认识，适当合理地简化并将之提炼成一个数学问题(上一节的两个例题)。另一个是方程的求解及结果分析。前一个阶段，我们将在“数学建模”课程中学习一些原则、方法和范例；第二个阶段即当前，我们主要简述关于微分方程求解的系列方法。,3.3 微分方程模型,y(n)(x)= f (x,y(x)

8、,y(x),.,y (n-1) (x),例如: y(x)=y(x)+snnx+x2+1,法则 f , 自变量 x 的取值范围是已知 , y(x) 待求,back,3.4 微分方程求解方法,在微积分的教程中，有一章专门介绍微分方程基本知识及求解方法。 但它们主要介绍解析解法，这些方法能够求解一些比较特殊类型 的微分方 程，而对大量较为复杂的微分方程，却难以求出解析解，这就需要使用数值 解法，在实际问题中数值解法是求解微分方程（组）的一种常用方法。本节 我们重点介绍数值解法。,3.4.1 数值解法,设待求解的函数 y(x) , 。满足微分方程,法则 f, 定义域均是已知的, 初始条件是给定的。 y

9、(x) 、 y(x) 均是未知的,例如: y(x)=y(x)+x2+1,y(0)=1.其中x的取值范围为0,5.求函数y(x) f 法则体现为加法,乘方. 例如：二元函数 g(s,t)=(snns)*ln(sqrt(t)+s2)-1/t. g法则是明确的 微分方程为 y(x)=g(x,y(x) ,y(1)=0 这就是本段要讨论的一类微分方程,给出自变量取值点列 xn ，步长为h=xn-xn-1.,x0,x1,x2,xn,xn+1,y0,y(x0)= y0 是微分方程的初始条件 利用y0和x1 我们来求y(x1)的近似值y1 利用y1和x2 我们来求y(x2)的近似值y2 利用y2和x3 我们来

10、求y(x3)的近似值y3 一般的利用yn和xn+1我们来求y(xn+1)的近似值yn+1,x3,设 数值解法的思路：,注意:我们没有得到函数y(x),而是对于取定的自变量序列 x0, x1, x2, xn, xn+1 得到 y(x0), y(x1), y(x2), y(xn) , y(xn+1) 的近似值 y0, y1 ,y2 , yn , yn+1 但是,近似程度是可控的。步长h越小,精度越高,数值解法一般只能得到微分方程的近似解,想,已知的,想求的,得到的,x0,x1,1.欧拉方法：欧拉法是一种最简单的数值解法。,曲,切,它的基本思路是,h =,在x0 , x1 上, y(x1) y(x0

11、)+h* f(x0 , y(x0) ,右边记为y1 于是有 y(x1) y1 在x1,x2 上, y(x2) y(x1)+h* f(x1,y(x1) y1 +h* f(x1, y1) 记y2= y1 +h* f(x1, y1) 于是有 y(x2) y2 在x2 , x3 上, y(x3) y(x2)+h* f(x2,y(x2) y2 +h* f(x1, y2) 一般的,我们有迭代公式 yn+1= yn+h* f(xn,yn) n=0,1, (3.4) y0 , y1 , y2 , yn , yn+1 就叫做微分方程的数值解,实验页,(3.5),若在小区间 上用差商代替导数y(x)时，取 (y(

12、xn+1)-y(xn) /h y(xn+1) 可得另一个欧拉公式,公式(3.4)被称为显式欧拉公式。也称为向前欧拉公式,y(xn+1) yn+1= yn+h* f(xn+1,yn+1),xn,xn+1,这一公式不易直接使用,其意义在于:与向前欧拉法结合产生改进的欧拉公式 yn+1=yn+(h/2)*(k1+k2) k1=f(xn,yn), k2=f( xn+1 , yn+h*k1 ).,(3.7),xn,xn+1,x0,x1,曲,前切,后切,割线斜率代替区间左端点切线斜率：,割线斜率代替区间右端点切线斜率：,一般地有,两式相加除以2,再把 yn+1 的近似值代入,整理得,yn+1 = yn+

13、(h/2)*(k1+k2) k1=f(xn,yn), k2=f( xn+1 , yn+h*k1 ).,(3.7),+,此式体现的近似，在几何上，明显地好于向前欧拉法。,例3.1 求解以下微分方程并带初值条件 y(x)=-y(x)+x+1 ,y(0)=1,取步长h=0.01 分别用向前欧拉法和改进欧拉法求解，并结合其解析解，对求解误差进行绘图比较。 解 首先用解析法 dsolve(Dy=-y+x+1,y(0)=1,x) 得 y(x)=x+exp(-x),绘图 x=0:0.01:1; y=x+exp(-x); plot(x,y) hold on,微分方程右边的二元函数为： f(x,y(x)=-y(

14、x)+x+1 因此 f(u,v)=-v+u+1 f(a,a+b)=-(a+b)+a+1 f(s+t,s-t)=-(s-t)+(s+t)+1,向前欧拉法: 因为 yn+1= yn+h* f(xn,yn) ,所以 yn+1= yn+h*(-yn + xn +1 ),f ( x , y(x)= - y(x) +x+1,h=0.01 ; x=0: h :1; % 共有snze(x,2)个分量 y(1)=1; % y的第一个分量 for n=1:snze(x,2) y(n+1)= y(n)+h*(-y(n) + x(n) +1 ) ; % y的其他分量 end plot(x,y),用改进的欧拉法: 公式

15、为 yn+1=yn+(h/2)*(k1+k2) k1=f(xn,yn), k2=f(xn+1,yn+h*k1),x=0:h:1; y(1)=1; for n=1:100 k1=-y(n)+x(n)+1; k2=-(y(n)+h*k1)+x(n+1)+1; y(n+1)=y(n)+(h/2)*(k1+k2); end plot(x,y),yn+1 =yn +h*f(xn,yn) 补充例题: y(x)=y+x2+1 ; y(0)=1 对比其解析解和数值解 y(x)=1/y ; y(1)=1 对比其解析解和数值解 y(x)=x/y2 ; y(1)=1 对比其解析解和数值解,计算结果表明,当步长h=0

16、. 01时, 三种方法的绘图结果是非常接近的 在迭代中，步长h越小，计算结果越精确。从迭代过程来看，迭代离开初 始点越远，误差越大。如图3.3,图3.3数值解曲线与精确解曲线比较图,因为y(x)=f(x,y(x),所以y(x)=fx+fy*y = fx+fy* f可以表达为x,y(x) 的函数,使用归纳法可证y(n)(x)也可以表达为x,y(x) 的函数.故,取泰勒公式前面若干项可近似得到公式 y(x+h) y(x)+h*(x,y(x),h),2.龙格-库塔方法 关于微分方程 y(x)=f(x,y(x),龙格-库塔方法（简称R-K方法）：由上式产生的迭代公式,则称以上迭代公式为p阶公式，p的大

17、小反映了截断误差的大小。 要得到一个p阶公式，关键在于如何选取 使之满足阶 的要求。常用的R-K公式,若,由泰勒公式:,其中c介于x与x+h之间,设对x值已求得y(x),在x点处,取步长h，现在来求y(x+h),yn+1=yn+ h*(xn,yn,h) n=0,1,.,误 差 估 计,1）2阶公式,中点公式,改进欧拉公式,2）3 阶公式,例题:y=y/x2 用系统命令dsolve和龙格-库塔3阶公式求解,在MATLAB软件中含有数值求解的系统函数，其原理就是龙格库塔方法,返回实验页,dsolve(Dy=y/x2,y(1)=1,x),y=1/exp(-1)*exp(-1./x);,x=0.01:

18、0.01:1;,解为：y(x)=1/exp(-1)*exp(-1/x),plot(x,y,g),hold on,考虑龙格-库塔方法,微分方程y=y/x2, y(1)=1 f(x,y)=y/x2,h=0.01; x=0.01:h:1; y(1)=1; for n=1:99 k1=h*y(n)/x(n)2; k2=h*(y(n)+1/2*k1)/(x(n)+1/2*h)2; k3=h*(y(n)-k1+2*k2)/(x(n)+h)2; y(n+1)=y(n)+1/6*(k1+4*k2+k3); end plot(x,y,r),用数值解法(龙格库塔3阶公式) 解微分方程,问题：,2) 对应问题的迭代

19、公式,解决方法： 1) 对函数 y= y(x) 的自变量的取值，给出点列,步长取为常量h,3) 编程的一般形式 h= 步长值 ; x= 初始值 : h : 终止值 ; y(1)= 初始值; for n= 1: (snze(x)-1) k1=h* f ( x(n) , y(n) ) ; k2=h* f ( x(n)+0.5*h , y(n)+0.5*k1 ) ; k3=h* f ( x(n)+h , y(n)-k1+2*k2 ) ; y(n+1)= y(n)+(1/6)*(k1+4*k2+k3) ; end plot(x,y),back,补充例题: y(x)=y+x2+1 ; y(0)=1 求数

20、值解 y(x)=1/y +x ; y(1)=1 求数值解 y(x)=x/y2 ; y(1)=1 求数值解,3.4.2 图解法,前面介绍的解析解法，尽管能够找到精确解，但是其适应面太窄；数值解法适应面广，却要损失精度，并且只能得到一些离散点处解的近似值；下面的图解法能从新的角度展示微分方程的全局解，即适应任意初始条件的解。,1.斜率场图示微分方程的解,设微分方程为 。我们在平面上取一个(待求函数 y(x) 可能经过的) 点(a,b). 把这一点(a,b) 代入微分方程,得 y(a)=g(a,b) 。这就是说，如果y(x)经过这一点，则我们就可以推知 y(x) 在(a,b)点的切线斜率。 从实际问

21、题出发，估计 y(x) 的取值区域，比如,对其中每一个点(a,b),我们都可以预先求出 y(x) 在该点的切线斜率 y(a)=g(a,b) 根据斜率值，我们作出y(x)在该点处的短切线。这种短直线布满整个预定范围所形成的图形就称为斜率场或方向场。它对于推测函数y(x)的走势是十分重要的。参见P60图形3.4,对矩形区域中的每一个点 (a,b) , 都给定了方程 y(x) = g(x,y) 的一个初始条件 y(a)=b 因此,我们解出的是“一族函数”.是从全局的角度描绘 y(x) 。 此“一族”不同于不定积分的彼“一族”。为什么？,因为 y(x)=g(x,y(x) 所以有,但，此处的 y(x)

22、是未知函数，处在一个含积分运算的方程中实验表明，y(x) 们的确不是平移关系。就此也可写一篇毕业论文！,下面,进一步说明 y(x) 与斜率场的对应关系:,假设，得到了一条 y(x) 曲线。进行如下处理：,1) 分：将曲线分割成若干段首尾相连的弧线段。 2) 替：将每一段弧线用左端点处短切线代替。 3) 短: 直线段较短,使之与相邻直线段不相连（断）.,这样，由微分方程的全部解曲线就可以得到某一区域的斜率场。通常由计算机编程来完成这个任务。,由全局解给出斜率场,设针对微分方程 y(x)=g(x,y)给出了斜率场。从一个初值 y(a)=b 出发，沿该点切线所指明的方向前进,按切线方向一段接一段地走

23、下去,就会弯弯曲曲地走出一条”路”.这条路就是关于微分方程 y(x)=g(x,y) ， y(a)=b的解曲线。,反过来由斜率场给出全局解,例 3.2用斜率场求解下列微分方程:,解 用MATLAB编程实现,syms s,t % 定义两个符号变量,用它们构成代数表达式 f=snn(s)*snn(t) ; % 由s,t 构成的代数表达式,需要时再赋值,a=16.0; b=16.0; %斜率场的横向长度，纵向长度,x0= - 8; y0=-8; % 矩形区域的左下角坐标,m=40;n=40; % 绘图区域横纵向划分区间,每个小矩形,clear s,t,x0,y0,a,b % 清除以前可能用过的变量,下

24、面要使用,% 区域,的左下角做一条小切线,%每个小矩形区域的横向长度,%每个小矩形区域的纵向长度,% 在同一个图形窗口绘许多小斜线,% 外循环,% s从左向右取左下角横坐标数值,% 内循环,处理一列中的小区域,% t从下向上取小矩形左下角纵坐标,%在(s,t)点计算斜率y(s),预备画切线,%下面计算小切线的终点坐标,% 取小切线的横向长为2/3*h,利用斜率d计算高度,2/3*h,y1-t,2/3*v,(s,t),% 如果y1-t2/3*v,则重新确定终点坐标.,x1=s+1/d*(2/3*v);,%取小切线的纵向长为2/3*v,由d算横坐标,y1-t,(s,t),2/3*v,?,end,e

25、nd,end,% 取小切线的横向长为2/3*h,利用斜率d计算y1,然后绘图,plot( s , x1 , t , t+2/3*v );,MATLAB关于求导数的命令 DiFF Differentnate. 求导 DiFF(S) differentnates a symbolnc expressnon S with respect to its free variable as determnned by FiNDSYM. 对符号表达式s求导数,其自变量由fnndsym确定 DnFF(S,v) or DnFF(S,sym(v) dnfferentnates S wnth respect to

26、v. 对s求导,关于变量v DnFF(S,n), for a posntnve nnteger n, dnfferentnates S n tnmes. S的n次导数 DnFF(S,v,n) and DnFF(S,n,v) are also acceptable. 两种表达形式均可以 Examples; x = sym(x); t = sym(t); diff(sin(x2) is 2*cos(x2)*x diff(t6,6) is 720. See also nnt, jacobnan, fnndsym. Reference page in Help browser doc sym/diff

27、,上机实验内容,用向前欧拉公式解微分方程 y=sinx+y 请与f=dsolve(Dy=snn(x)+y,y(0)=0,x) 的结果做数值比较 2. 编程序用龙格-库塔3阶公式解微分方程 y=y*log(x) 请与 f=dsolve(Dy=y*log(x),y(0)=1,x) 的结果做数值比较 本周作业 用斜率场描绘微分方程的解函数曲线 (要求编出通用性较强的程序) y=sin(x*log(y) 问题:如何提高斜率场的可读性,课堂练习:用龙格库塔公式和斜率场解微分方程,取y(1)= 1,1) 图3.4是否可以看出该微分方程的通解或者积分曲线簇的形状? 2) 斜率场的另一个应用是对那些不能用初等

28、函数表示的积分,可作出它们的积分曲线.,提示,想,意义:,虽然函数,可积，但其原函数不是初等函数。,因此，采用数值解或斜率场描绘,很有意义。,相平面轨迹 上个世纪初，意大利科学家在研究地中海地区的渔业资源时，发现鲨鱼的捕获数量与鱼类的捕获数量有一定的关系。其后较长的时间，对这一问题的持续研究，就形成了我们今天的经典范例。 设 t 时刻鲨鱼的数量为y(t),鱼饵的数量为x(t)。客观上，鲨鱼的数量与鱼饵的数量是存在关系的。这一关系近似地描述为,鱼饵的增长率与鱼饵基数x(t)成正比,同时受制于当前数量 x(t)(自然资源和环境条件有制约作用)和鲨鱼数量y(t)的乘积; 鲨鱼的增长率受制于当前数量y

29、(t) (食物链中最高一环),同时获益于鱼饵x(t)数量与鲨鱼基数y(t)的乘积.,鲨鱼数量和鱼饵数量构成问题研究的两个相。 对每一个t 的取值, 求出数组(x( t ) ,y( t )，把它作为一系列平面点描绘在XOY平面上。以两个相为点得到的平面曲线叫做相平面轨迹图,x(t)= rx(t) - ax(t)y(t) y(t)= - sy(t)+bx(t)y(t),t x(t) y(t),t1 x1 y1,t0 x0 y0,t2 x(t2) y(t2),为认识x=x(t) , y=y(t) , 考虑以下三种方法作出其相平面轨迹图.,x(t) X相,y(t) Y相,1)直接用微分方程组的数值解法

30、 调用系统提供的ode45命令.这个命 令的原理是使用4阶5阶的龙格库塔 公式。设初始值x0=x(t0),y0=y(t0) 取参数点列t0,t1,t2,.tn,tn+1, 步长 为h,于是可相应解出数值解 (x(t) ,y(t) 其横、纵坐标分别定义为列向量X,Y； 利用X，Y的分量作为坐标，就可以 在XOY平面上画出相平面轨迹图。,X Y,tn x(tn) y(tn),tn+1 x(tn+1) y(tn+1),下面通过编程绘出当r=2,s=1,a=1,b=1时，初始值为 (x0,y0)=1 0.3,1 0.5,1 1.5 时绘出的相轨线,为了得到数值解，必须给出适当的初值条件。同时，为了得到

31、多条相平面轨迹，便于和其他方式得到的相平面轨迹图进行比较，这里我们给出了系列初值下的解及其相平面轨迹图，如图3.5.,x(t)= rx(t) - ax(t)y(t) y(t)= - sy(t)+bx(t)y(t),1)首先建立函数m文件(fnsh.m)，把具体的f,g反映到函数文件中 functnon xdot=fnsh(t,x) r=2;s=1;a=1,b=1; xdot=r*x(1)-a*x(1)*x(2); -s*x(2)+b*x(1)*x(2); 2)在工作空间执行以下程序 hold on for k=1:7 %对不同的初始值循环 ts=0:0.01:5; %时间序列 x0=1 0.1

32、+0.2*k; %对应x(t0),y(t0) t,x=ode45(fnsh,ts,x0); plot(x(:,1),x(:,2) end axis(0 4 0 4);xlabel(X);ylabel(Y),fnsh(t,x),ts x0,xdot,ode45的迭代示意图,输出,2)利用斜率场作图。 取定参数t的值,在t 点给出微分值dt.于是t和dt均是常数. dx(t) =f(x(t),y(t)dt=(r*x(t)-a*x(t)*y(t)dt dy(t) =g(x(t),y(t)dt=( -s*y(t)+b*x(t)*y(t)dt dy/dx=(-s*y+b*x*y)/(r*x-a*x*y)

33、 此式给出了关于y=f(x) 的微分方程 y(x)=g(x,y),它是一个一阶微分方程。由前一小节，在平面直角坐标系中可以作出斜率场图形，如图当参数是 r=2,s=1,a=1,b=1 时的情形。 y(x)=(-y+x*y)/(2*x-x*y) 修改绘制斜率场的通用程序并调用,得,图3.6通过斜率场得到的相平面轨迹图,3）利用二元函数的等值线作图,将以上微分方程通过分离变量 dy/dx=(-s*y+b*x*y)/(r*x-a*x*y) dy=(-s*+b*x*)*y/x*(r*-a*y)*dx (r-a*y) /ydy=(-s+b*x)/x*dx r/y*dy-a*dy=-s/x*dx+bdx

34、等式两边对x求不定积分,得 rlny-ay+slnx-bx=k (*),合理取定k的值,由此解出的(x,y)构成一条曲线,绘图后,从曲线中可以看出x与y的关系.,具体做法：,令 f(x,y)=rlny-ay+slnx-bx 适当地取x,y的值,生成网格X,Y,然后绘制空间曲面.,取定一个k值,让f(x,y)=k 便可得到一条平面曲线(x,y)|f(x,y)=k,该曲线(x,y,k)|f(x,y)=k反映了x与y的关系,X Y=meshgrnd(0.01:0.1:4 , 0.01:0.1:4); Z=2*log(Y)-Y+log(X)-X; mesh(X,Y,Z),contour(X,Y,Z,2

35、0); xlabel(x); ylabel(y);,图3.7通过等值线得到的相平面轨迹图,1）随着时间的推移，相平面轨迹图中轨迹的走向按逆时针还是顺 时针方向？即,随着t的变化,(x,y)如何变化? 2）三种方法作出的相平面轨迹图形在形式上似乎有一些差别，原 因何在？本质上是否一致？相平面轨迹图的中心在何处？,想,3.5 MATLAB软件求解,前面已经介绍了微分方程（组）的求解方法数值解法和图解法。这些方法计算工作量大，需要借助于计算机,3.5.1 解析解,求微分方程（组）的解析解的MATLAB命令,实现。下面简单介绍MATLAB软件在此领域的应用。,其中 表示第n个方程， 表示微分方程（组）

36、中的自变量，默认时自变量为t.,例3.3 求解一阶微分方程,求通解 输入：dsolve(Dy=1+y2) 输出：tan(t-c1) 求特解 输入：dsolve(Dy=1+y2,y(0)=1,x) 输出：tan(x+1/4*pi),例3.4 求解二阶微分方程 x2y+xy+(x2-n2)y=0, n=1/2 y(/2)=2,y(/2)=-2/, 方程两边同除以x2,然后输入： 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)*sin(x)/x(1/2) 化简输出结果，输入

37、： 计算结果为：,例3.5 求微分方程组 的解析解。求通解,输入：,输出：,求特解，输入：,输出：,3.5.2 数值解,设微分方程的形式为 其中t为自变量，y为因变量。 在MATLAB5.0版本中，用2阶3阶龙格-库塔公式和4阶5阶龙格-库塔公式的程序分别为,Options用于设定误差限（缺省时设定相对误差是 绝对误差 ）。程序为,3）对于等步长时用 则输出由区间 的等分点给出。,其中F是由微分方程(组)写成的M文件名。“时间序列” ts的取法有几种，1）当 分别表示自变量的初值和终值；2）若 则时间序列为指定时刻,为函数的初值。,ode23是微分方程组数值解的低阶方法，ode45为较 高阶方

38、法，与ode23类似。另外还有一些其他方法，如求 解非刚性微分方程组的可变阶方法ode123.,为输出矩阵，分别表示自变量t和因变量y的取值。,例3.6 求例3.4中方程 的数值解。,解 记 则转化为一阶微分方程组,两边同除以x2,得,y1(x)=f(x,y1(x),y2(x) y2(x)=g(x,y1(x),y2(x),一般形式,y=-1/x*y+(n/x)2-1)*y y(x)=?,求解如下：首先建立M函数记录微分方程组 functnon sty=jie3(x,y) sty=y(2) ; -y(2)/x+(0.5/x)2-1)*y(1);,存盘退出。 其中x是自变量序列,y用于接收初始条件

39、向量2 -2/pi.,输入：x,y=ode23(jie3,pi/2,pi, 2 -2/pi),自变量序列,初始条件向量,输出：,Y=2.000 -0.6366 1.9758 -0.6869 1.8518 -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,y2,作图程序：,输出结果见图3.8,无论是解析解还是数值解，都不如图形解直观明了， 在MATLAB下用

40、什么命令可以用图形显示解的解或 数值解？,如果微分方程（组）的解析解为 则可以 用MATLAB函数fplot作出函数曲线 的图 形，具体用法,想,提示,其中fun给出函数表达式； 限定坐标轴的大小。例,如果得到微分方程（组）的数值解 用MATLAB 函数 直接作出图形，其中x和y为向量（或矩 阵）。,fplot(xx, 0.01 1 0 2),3.6.1 问题及假设 我们将研究的计算机通信网络系统设定为无冗余的防火墙协议系统。,3.6 范例：状态转移方程组模型,计算机网络可靠性分析问题随着计算机通信网络系统特别是lnternet网络日益广泛的应用，提高系统的可靠性意义重大。研究和分析具有实用性

41、的高可靠计算机通信网络系统，是国际上非常活跃的一个研究方向。,无冗余是指无备份线路和计算机,防火墙协议系统在此是指一个相对独立的系统,计算机通信网络系统随时可能发生三个事件无故障、间歇故障和永久故障。,无故障工作,带故障工作,图3.9 无冗余的防火墙协议系统状态转移过程,不工作,因此，一般处于三种工作状态：无故障工作、带故障工作和不工作。这三种状态之间的转移过程如图3.9所示。要求建立该系统的状态转移模型，并进行可靠性分析。,3.6.2 分析与模型,该问题属于状态转移问题，用 分别表示系统处于无故障工作、带故障工作和不工作三种状态的概率，则p1(t)+p2(t)+p3(t)=1,利用马尔科夫状

42、态转移原理，有以下状态转移方程组：,结构与参数,初始条件 参数取值 这是一个带参数的微分方程组模型。,3.6.3 模型求解,求解这个带参数的微分方程组模型有两种方法， 一是用特征根法求解析解，另一个是用数值解法求 数值解。分别求解如下：,假定模型中的参数取下限值，即,1.解析解法 利用常系数线性微分方程组的特征根法。 MATLAB程序：,输出结果： 特征值（0，-0.0102, -0.0001) 与之对应的特征向量（0 0 1）；,根据特征值和特征向量写出其解析解（通解）， 并利用初始条件 得到特解,2.数值解法 MATLAB程序如下 首先编写微分方程的M文件 eqs0.m functnon

43、xdot=eqs0(t,p,flag,lp,lt,gm) xdot=-(lp+lt) gm 0 ; lt/2 -(gm+lp+lt ) 0; lp+lt/2 lp+lt 0*p(1) p(2) p(3);,在工作空间执行以下程序： ts=0 1000; p0=1;0;0; lp=10(-5); lt=10(-4); gm=0.01; t p=ode23(eqs0,ts,p0, ,lp,lt,gm); plot(t,1-p(:,3) %可靠度记为R(t)=p1(t)+p2(t),传送较多参数时,需要这种标志,输出结果如图3.10.,1,0.99,0.98,0.97,0.96,0.95,0.94,

44、可 靠 度 R(t),0,200,400,600,800,1000,参数取值,时间t(小时),图3.10 系统正常运行的可靠边度曲线,计算结果表明：在时间t=1000小时情况下，,显示系统工作的概率为94.18%.,从图3.10中还可以看出系统可靠性变化更加细致的情况，何时可靠达到99%，何时达到98%等等。,3.敏感性分析,当参数取值在以下范围,内变化时，系统的可靠性情况如何呢？留给读者自己分析。,3.7 实 验,3.7.1 实验一 求解下列微分方程(组),1.简单微分方程,实际对象的数学 模型包括两部分 模型结构 模型参数 当模型结构和模型 参数发生变化时， 模型求解结果也 应发生相应的变

45、 化。这种变化的 定量分析就叫做 敏感性分析,2.特殊微分方程,1）Duffings方程: u+u-u3=0,y-y=x-2 , y(0)=2 , y(1)=1,2) (1+x)y=2y-4 ,y(0)=0 , y(1)-2y(1)=0,记 u 为 x , 于是 方程可以变形为，,给出其相平面图，参数 的取值可以考虑-1/4,0.1,1/4 第一步 ,建立M文件 functnon f=u(t,x,flag,e) f=x(2);-x(1)+e*x(1)3; % x(1)表 x(t) , x(2)表 y(t),x = y y = - x + x3,第二步编一个程序调用ode45命令 hold on

46、 e=-1/4 0.1 1/4; color=r ,g, b; for n=1:3 ts=0:0.01:10; x0=1 1; t f=ode45(u,ts,x0, ,e(n); plot(f(:,1),f(:,2),color(n); end,2)单摆运动方程 u+8sinu=0 用图形表示其解u(t)函数。 令x=u(t) y=x得 x=y y= -8*sin(x) 建立M文件 function f=t3(t,x) f=x(2);-8*sin(x(1); 调用ode45命令并绘图t f=ode45(t3,0 10,0 1); plot(t,f(:,1),3 、微分方程组 1）线性微分方程组

47、x=Ax其中矩阵A的定义如下： 2 1;-3 6 0 1 0;4 3 -4;1 2 1 注意：x是由若干个函数做成的向量编写M文件,functnon f=fx(t,x) f=2 1;-3 6*x; 调用ode45函数，并绘图 ts=0,1; x=1;2; t x=ode45(fx,ts,x) plot(t,x(:,1),r, t,x(:,2),g),functnon f=fx2(t,x) f=0 1 0;4 3 -4;1 2 1*x; 调用ode45函数，并绘图 ts=0,1; x=1;2;3; t x=ode45(fx2,ts,x) plot(t,x(:,1),r, t,x(:,2),g,

48、t,x(:,3),b),2）非线性微分方程组 x=-y-x2 y=x 建立M文件 functnon y=l1(t,x) y=-x(2)-x(1)2 ; x(1); 调用ode45命令并绘图 t y=ode45(l1,0 1,1 2); plot(t,y(:,1),r,t,y(:,2),g),建立M文件 functnon y=l2(t,x) y=x(1)-4*sqrt(abs(x(1)*x(2) ;-x(2)+4*sqrt(abs(x(1)*x(2); 调用ode45命令并绘图 t y=ode45(l2,0 1,1 2); plot(t,y(:,1),r,t,y(:,2),g),b) x= x-

49、4x*sqrt(abs(xy) y=-y+4x*sqrt(abs(xy),3.7.2实验二 盐水的混合问题,一个圆柱形的容器，内装350升均匀混合的盐水溶液，如果纯水以每秒14升的速度从容器顶部流入，同时，容器内的混合的盐水以每秒10.5升的速度从容器底部流出。开始时，容器内盐的含量为7千克。求经过时间t后容器内盐的含量。 问题分析:记t时刻,容器内盐的含量为q(t)千克,溶液总量为p(t)升,我们考虑在微小时间段 dt 上含盐量和溶液总量的变化。 我们把dt取得足够小，使得不断注入纯水，流出盐溶液的动态过程可视为流出盐溶液，再注入纯水。,模型假设：盐溶液的变化是均匀的；在dt足够小时，把动态

50、过程静态化 建立模型： 随时间t变化的溶液总量为p(t)=350+(14-10.5)t (升) t时刻容器内的含盐比率为q(t)/p(t) (千克/升) dt时间段上减少的盐溶液为10.5*dt (升) dt时间段上减少的盐量为q(t)/p(t)*10.5*dt (千克) 因此，有 dq(t)=- q(t)/p(t)*10.5*dt,即q(t)适合微分方程 q(t)=- q(t)/p(t)*10.5=-q(t)/ (350+(14-10.5)t )*10.5 q(0)=7 模型求解： dsolve(Dq=-3*q/(100+t),q(0)=7,t) ans = 7000000/(100+t)3

51、 图形表示: fplot(7000000/(100+t)3,0 200 0 10),孟德尔（Mendel)第一定律：配子的基因是从其父辈的两个基因型中随机地选择的。 实际应用中，将比例(种群角度)作为概率(个体角度) Pk(A)=ProbAA或Aa;包含A Pk(a)=Probaa,仅含aa 并记Xk=Pk(a). 得到如下遗传模型： 1）致死基因遗传模型： Xk+1 = Xk / (1+Xk ) 讨论Xk随k变化的变化趋势 2）自然选择基因遗传模型： 其中： r1和r2分别表示在总人口数量中，新生儿基因 为(AA或Aa)和(aa)所占的比例。对不同的取值，讨论Xk的变化趋势. 选取初值X0=

52、0.9 3)突变基因遗传模型： 其中： 为A突变为a的概率（比例一般为10-5-10-6）。对不同的 讨论Xk的变化情 况，考虑初值X0=0.1,模型分析与假设: 设第k代从第k-1代继承含A的两个基因的概率(可能性)为 Pk(A)=ProbAA或Aa 第k代从第k-1代继承不含A的两个基因的概率(可能性)为 Pk(A)=Probaa aa基因组合构成致死基因 1)致死基因遗传模型 Xk+1 = Xk / (1+Xk ) 即第k+1代继承致死基因的概率为第k代的代数组合 三个认识以及模型建立: Xk为一个离散的数列.若有一个连续,可微的函数y(t),使得 y(t)= Xt t=0,1,2,n

53、那末,我们可以先研究y(t) t取连续值,然后利用y(t)的解析性质来认识Xt .实际贯彻这一认识时,不必提及y(t),把Xt改为X(t)并视为连续可微的即可.因此,致死基因遗传模型写为 x(t+1)=x(t)/(1+x(t) 将迭代公式改为增量公式 x(t+1)-x(t)=-x(t)2/(1+x(t) 这是为了建立关于x(t)的微分方程,dt=1 在人类繁衍这一问题上,给一微小的合理的”代”数增量dt,只能是dt=1.因此 x(t+1)-x(t)= -x(t)2/(1+x(t) = -x(t)2/(1+x(t) *dt 相应于右边的dt,左边有 x(t+1)-x(t)=dx(t) 于是得到关

54、于x(t)的微分方程 x(t)= -x(t)2/(1+x(t) x(0)=0.33 初始条件,假设第一子代从父代继承致死基因的概率为0.33 模型求解: dsolve(Dx= -x2/(1+x),x(0)=0.33,t) ans = exp(lambertw(exp(t+log(100/33)+100/33)-t-log(100/33)-100/33) 解的分析认识:,解函数中,引用了兰姆伯特函数。关于兰姆伯特函数，如果 w=lambertw(x) 则w*exp(w) = x. 这给出了兰姆伯特函数的自变量和函数值的对应法则 采用类似的方法可以做2),3).留给同学们做实验,3.7.4 实验四

55、 老鼠觅食 设有连续的很多个小老鼠笼子（正方形），它们首尾相连。在其前后两边的中央都开有一个洞，可供老鼠自由进出。并在右边放置鼠粮，左边未放鼠粮。老鼠在笼子里面只能够沿着笼子边沿（正方形的四条边）靠左边或右边向前通过。沿左边则吃不到鼠粮，只有沿右边才能够吃到鼠粮。在每个鼠笼子里，老鼠随机地选择左右之一向前行进。 1）奖励型(能感知奖励的老鼠)：如果老鼠沿右边吃到鼠粮后，则下次将毫不犹豫地沿右边，如果沿左边未吃到鼠粮，则下次将以1-的概率向左。,就这两种情况，分别建立并求解老鼠在第n次进入鼠笼子时向右能够吃到鼠粮的概率，并考察其无穷趋势。 当然有兴趣的读者可以考虑学习型(设置多次通过)，设想老鼠

56、具有一定的学习记忆能力。并将其在模型中反映出来。 1)模型分析与假设:,2）奖惩兼顾型：如果向右吃到鼠粮后，则下次向右的概率为1- ;如果向左未吃到鼠粮，则下次向左的概率为1-.,实验环境图示:,从进入第1个笼子到离开第n个笼子,老鼠在每一个笼子都有靠左或靠右两种选择, 因此,一共有2n种通过方式. 由于, 在笼子的右侧放置了鼠粮, 而老鼠能感知奖励,所以上述2n种通过方式不是 等可能的.我们下面考虑在第n个房间是靠右通过的那些方式. 把老鼠采用每一种通过方式的可能性即概率计算出来,然后相加,就得到”在第n个 笼子向右吃到鼠粮的概率”,根据实验设置,用字符串来表达老鼠可能选择的路径 1 2 3 n-2 n-1 n . L L L L L R . L L L L R R . L L L R R R . L L R R R R . L R R R R R . R R R R R