数学实验第四次讲稿.ppt

数学实验第四次讲稿 Still waters run deep.流静水深流静水深,人静心深人静心深 Where there is life,there is hope。有生命必有希望。有生命必有希望 1)掌握微分方程求解的三种解法：解析 法、数值解法以及图形表示解的方法;2)学会使用MATLAB软件求解析解、数值解 和图形解；3）通过范例学习怎样建立微分方程模型和 分析问题的思想；实验目的一，建立微分方程数学中有一大类问题，我们称之为微分方程，它是含有导数的方程，比如 f(x,V(x),V(x)=0 是一阶微分方程，而 F(x,y,y,y(n)=0 是n阶微分方程，等等。一，建立微分方程微分方程在现代科学的每一个领域都有广泛的应用，比如力学，运动学，电学，经济学，生物学，自动控制，化学等等，都可以看到大量利用微分方程表示的事物变化规律，这也体现了微分方程的重要性。在高等数学中我们对微分方程有一定的认识，比如怎样求解一些简单的微分方程，对于一些细心的同学，可能想到对于复杂的微分方程的求解是一个重要的问题，但 一，建立微分方程 是，可能忽略一个同样重要的问题，那就是如何建立微分方程。同样是方程，其建立的过程一定程度上和前一部分方程的建立有相似之处，但不同之处是微分方程的建立过程中一般会涉及到变化率的概念（往往和时间有关）。比如运动学中速度是位移对时间的变化率，加速度是速度的变化率，电学中电流是电量的变化率，人口学中人口增长率是人口数对时间的变化率等等。一，建立微分方程 下面用一些常见的问题说明微分方程建模的基本方法。例子1：放射性物质的衰变。物理学家发现，放射性物质衰变的速度和物质的量成正比。例子2：人口（生物数量）增长问题。通过观察发现，人口增长和现有人口数量成正比。一，建立微分方程 例子3：热传导问题。牛顿定理表述热传导为：高温物体热量的散失速度和温度差成比例。例子4：物体的运动问题。物体运动可以用牛顿三大定律描述。例子5：万有引力定理。开普勒发现地球沿着以太阳为焦点的椭圆 轨道运动，运动速度满足单位时间所扫过的面积相等。一，建立微分方程 例子6：传染病问题。医学经验表明，传染病有几个关键要素，易感人群，带病人群，而传染往往通过两种人群的接触导致，不禁和易感人群和带病人群有关，还有他们的接触有关。例子7：战争问题。交战双方的伤亡虽然比较复杂，但是大致遵循这样一些规律，比如武器的杀伤力，防御力，还有双方的人数对比等因素。二，微分方程的解法解析方法 对于一些比较简单的微分方程，可以通过一些数学技巧解出，比如高等数学上接触的一些方程：可分离变量的方程，齐次方程，一阶线性微分方程，一些特殊的二阶常系数微分方程等等。数值方法 能得到解析解的微分方程毕竟只是少数，对于从实际问题中提取的大量微分方程，无法得到方程的解析解，数值方法可以说是实际问题中必不可少的手段。二，微分方程的解法 在介绍具体的方法之前，有必要了解一般方程的形式和对解的要求。常微分方程分成微分方程和微分方程组。比如：F(x,y,y,y(n)=0 是一个隐式的高阶方程，而 是一阶微分方程组。二，微分方程的解法关于微分方程的解，就是满足方程的一个函数族（或者一条曲线族）。我们又称其为微分方程的通解。用得更加广泛的是满足特定条件的解，我们称其为特解。比如二，微分方程的解法之解析方法解析解dsolve(eqn1,eqn2,c1,var1,)微分方程组初值条件变量组注意：y Dy，y D2y 自变量名可以省略，默认变量名t。dsolve(eq1,eq2,.,cond1,cond2,.,v)symbolically solves the ordinary differential equation(s)specified by eq1,eq2,.using v as the independent variable and the boundary and/or initial condition(s)specified by cond1,cond2,.The default independent variable is t.The letter D denotes differentiation with respect to the independent variable;with the primary default,this is d/dx.A D followed by a digit denotes repeated differentiation.For example,D2 is d2/dx2.Any character immediately following a differentiation operator is a dependent variable.For example,D3y denotes the third derivative of y(x)or y(t).Initial/boundary conditions are specified with equations like y(a)=b or Dy(a)=b,where y is a dependent variable and a and b are constants.If the number of initial conditions specified is less than the number of dependent variables,the resulting solutions will contain the arbitrary constants C1,C2,.例输入：y=dsolve(Dy=1+y2)y1=dsolve(Dy=1+y2,y(0)=1,x)输出：y=tan(t-C1)（通解，一簇曲线）y1=tan(x+1/4*pi)（特解，一条曲线）二，微分方程的解法之解析方法例 常系数的二阶微分方程y=dsolve(D2y-2*Dy-3*y=0,x)y=dsolve(D2y-2*Dy-3*y=0,y(0)=1,Dy(0)=0,x)输入：x=dsolve(D2x-(1-x2)*Dx+x=0,x(0)=3,Dx(0)=0)上述两例的计算结果怎样？由此得出什么结论？例 非常系数的二阶微分方程例无解析表达式！x=dsolve(Dx)2+x2=1,x(0)=0)例 非线性微分方程x=sin(t)-sin(t)若欲求解的某个数值解，如何求解？t=pi/2;eval(x)二，微分方程的解法之解析方法输入：x,y=dsolve(Dx=3*x+4*y,Dy=-4*x+3*y)x,y=dsolve(Dx=3*x+4*y,Dy=-4*x+3*y,x(0)=0,y(0)=1)例MATLABMATLAB软件实现软件实现返 回二，微分方程的解法之解析方法对于大量的微分方程，只能得到其数值解，一般而言，得到的解是方程的一个特解的近似。求微分方程数值解的方法很多，比如：欧拉法，龙格库塔法等。其基本思想就是通过已知点得到函数值，并用该函数值代替一个小区间上函数的导数，得到在该区间上的一条直线，并用该直线作为方程特解的近似。有兴趣的同学可以参考微分方程数值解方面的著作。二，微分方程的解法之数值方法ode23 ode45 ode113ode15sode23s由待解方程写成的m-函数文件ts=t0，tf，t0、tf为自变量的初值和终值函数的初值ode23：组合的2/3阶龙格-库塔算法ode45：运用组合的4/5阶龙格-库塔算法自变量值函数值用于设定误差限(缺省时设定相对误差10-3,绝对误差10-6),命令为：options=odeset（reltol,rt,abstol,at）,rt，at：分别为设定的相对误差和绝对误差.Matlab软件计算数值解t，x=solver（f,ts,x0,options）1）首先建立M-文件 （weif.m）function f=weif(x,y)f=-y+x+1;2）求解：x，y=ode23(weif,0,1,1)3)作图形:plot(x,y,r);4)与精确解进行比较 hold on ezplot(x+exp(-x),0,1)例1 y=-y+x+1，y(0)=1标准形式：y=f(x,y)范例 1、在解n个未知函数的方程组时，x0和x均为n维向量，m-函数文件中的待解方程组应以x的分量形式写成.2、使用Matlab软件求数值解时，高阶微分方程必须等价地变换成一阶微分方程组.注意注意:注意11、建立、建立M文件函数文件函数 function xdot=fun(t,x)xdot=f1(t,x(1),x(2)；f2(t,x(1),x(2)；2、数值计算、数值计算（执行以下命令）（执行以下命令）t,x=ode23(t,x=ode23(funfun,t,t0 0,t,tf f,x,x0 0,y,y0 0)注意：执行命令不能写在注意：执行命令不能写在MM函数文件中。函数文件中。xd(1)=f1(t,x(1),x(2);xd(2)=f2(t,x(1),x(2);xdot=xd;%列向量列向量例如：例如：令令注意2：function xdot=fun1(t,x)（fun1.m)xdot=f(t,x(1),x(2)；x(1)；t,x=ode23(t,x=ode23(fun1fun1,t,t0 0,t,tf f,x,x0 0,y,y0 0)M-M-文件函数如文件函数如何写呢？何写呢？注意：注意：y(t)y(t)是是原方程的解。原方程的解。x(t)x(t)只是中间只是中间变量。变量。如果方程形式是：如果方程形式是：z=f(t,z,z)?返 回例例2 Van der pol 方程方程：令令 y1=x(t),y2=x(t);该方程是否有解析解？范例（1）编写M文件(文件名为 vdpol.m):function yp=vdpol(t,y);global a;yp=y(2);a*(1-y(1)2)*y(2)-y(1);（2）编写程序如下:（vdj.m）global a;a=1;t,y=ode23(vdpol,0,20,3,0);y1=y(:,1);%原方程的解 y2=y(:,2);plot(t,y1,t,y2,-)%y1(t),y2(t)曲线图 pause,plot(y1,y2),grid,%相轨迹图，即y2(y1)曲线这里使用了全局变这里使用了全局变量进行参数传递，量进行参数传递，对于变化的参数对于变化的参数a，可以在主程序中加可以在主程序中加以修改。另外，当以修改。另外，当a的取值比较大时，的取值比较大时，曲线非常陡峭，不曲线非常陡峭，不能使用能使用ode45求解，求解，这是使用这是使用ode15s。蓝色曲线 y（1）；（原方程解）红色曲线 y（2）；计算结果范例范例例3 考虑Lorenz模型：其中参数=8/3，=10，=28解：1）编写M函数文件(lorenz.m);2)数值求解并画三维空间的相平面轨线；(ltest.m)范例1、lorenz.mfunction xdot=lorenz(t,x)xdot=8/3*x(1)+x(2)*x(3);-10*(x(2)-x(3);-x(2)*x(3)+28*x(2)-x(3);2、ltest.mx0=0 0 0.1;t,x=ode45(lorenz,0,10,x0);plot(t,x(:,1),-,t,x(:,2),*,t,x(:,3),+)pauseplot3(x(:,1),x(:,2),x(:,3),grid on计算结果如下图范例图中，x1的图形为实线(蓝），x2的图形为“*”线（绿），x3的图形为“+”线（红）.取t0，tf=0，10若自变量区间取0，20、0，40，计算结果如下范例曲线呈震荡发散状三维图形的混沌状ltest.m观察结果 1、该曲线包含两个“圆盘”，每一个都是由螺线形轨道构成。某些轨道几乎是垂直地离开圆盘中一个而进入另一个。2、随着t的增加，x(t)先绕一个圆盘几圈，然后“跳”到另一个圆盘中。绕第二个圆盘几圈，又跳回原来的圆盘。并以这样的方式继续下去，在每个圆盘上绕的圈数是随机的。思考：该空间曲线与初始点x0的选择有关吗？1）x0=0 0.1 0.1;t0,tf=0,30;解向量y2）x00=0.01 0.11 0.11;t0,tf=0,30;解向量x y x=(y1-x1,y2-x2,y3-x3).注：这是这两个向量注：这是这两个向量必须是同维数，并且所取的时间节点应该相同。解决的办法是必须是同维数，并且所取的时间节点应该相同。解决的办法是将将t0,tf直接定义为如下的形式，比如直接定义为如下的形式，比如0：0.1：30。三，微分方程建模实例 1追击路线问题：一艘缉私舰雷达发现距一艘缉私舰雷达发现距c km处有一艘走私船正以匀速处有一艘走私船正以匀速 a 沿直线行驶。沿直线行驶。缉私舰立即以最大的速度缉私舰立即以最大的速度 b 追赶，若用雷追赶，若用雷达进行跟踪，保持船的瞬时速度方向始终达进行跟踪，保持船的瞬时速度方向始终指向走私船，试求缉私舰追逐路线和追上指向走私船，试求缉私舰追逐路线和追上的时间。的时间。0(c,0)xyDR假设1、建立坐标系；2、缉私船在(c,0)处发现走私船在(0,0)处；3、走私船逃跑方向为y轴方向；4、t时刻，走私船到达R(0,at)缉私舰到达D(x,y)；1 追击路线问题0(c,0)xyDRyx(0,at)化简：对 t求导化简：（1）1 追击路线问题(2)结合(1)、(2)得到如下微分方程：是否存在解析解？1 追击路线问题1）当a b时，k=a/b b时，显然缉私船不可能追上走私船。讨论当a=b时呢？思考考虑数值求解法1 追击路线问题敌艇逃跑速度：a=60公里/小时；缉私舰追击速度：b=80公里/小时缉私舰发现敌艇时的相距：c=500公里；计算得：k=(a/b)=0.75;编程计算，并画出追击路线和追击所花费的时间。假设1 追击路线问题转化成一阶微分方程组2.Zhui.mx,y=ode23(zhuiji,500,1,0,0);y1=y(:,1);plot(x,y1)1.Zhuiji.mfunction f=zhuiji(x,y)f=y(2);0.75*sqrt(1-y(2)2)/x;MATLAB程序1 追击路线问题262.7时间：其中a=60(公里/小时)(小时)1 追击路线问题缉私舰追击路线敌艇逃跑路线思考1、在本问题的求解过程中，假定了走私艇的逃跑方向是正北方向，而初始缉私艇的位置在x轴正向。如果放宽这个假定，也就是当这个夹角是任意角度时，如何建立方程进行求解。2，如果有多个走私艇在一个位置上进行交易，而缉私艇向该方向追赶。这些走私艇向不同方向四散逃走，问如何安排追赶路线？1 追击路线问题返 回范例：范例：2 地中海鲨鱼问题地中海鲨鱼问题 意大利生物学家Ancona曾致力于鱼类种群相互制约关系的研究，他从第一次世界大战期间,地中海各港口捕获的几种鱼类捕获量百分比的资料中，发现鲨鱼的比例有明显增加（见下表）。但为何鲨鱼的比例大幅增加呢？生物学家Ancona无法解释这个现象，于是求助于著名的意大利数学家V.Volterra，希望建立一个食饵捕食系统的数学模型，定量地回答这个问题.因为战争使捕鱼量下降，食用鱼增加，显然鲨鱼也随之增加。战争为什么使鲨鱼数量增加？是什么原因？想2 地中海鲨鱼问题地中海鲨鱼问题x1(t),x2(t)分别是食饵、捕食者在t时刻的数量；r1,r2是食饵、捕食者的固有增长率；1是捕食者掠取食饵的能力，2是食饵对捕食者的供养能力；1、符号说明：食饵捕食系统的数学模型 捕食者的存在使食饵的增长率降低，假设降低的程度与捕食者数量成正比，即2、基本假设：食饵对捕食者的数量x2起到增长的作用，其程度与食饵数量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.02*x(1);2、建立主程序、建立主程序shark.m如下：如下：t,x=ode45(shier,0 15,25 2);plot(t,x(:,1),-,t,x(:,2),*)plot(x(:,1),x(:,2)编制程序如下：食饵捕食系统的数学模型求解结果：求解结果：图2 反映了x1(t)与x2(t)的关系。图2图1x1(t)x2(t)食饵捕食系统的数学模型食饵-捕食者模型的数值解曲线初值：2x2(t)最大值最大值28.4x1(t)的最大值的最大值 99.325周期点10.7食饵捕食系统的数学模型食饵-捕食者模型的相图（x1,x2）p 重心相轨线是封闭曲线等价于x1(t),x2(t)是周期函数。重心坐标如何计算？(25,10)食饵捕食系统的数学模型模型（二）考虑人工捕获的情况 假设人工捕获能力系数为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),*)食饵捕食系统的数学模型function dx=shier1(t,x)dx=zeros(2,1);dx(1)=x(1)*(0.7-0.1*x(2);dx(2)=x(2)*(-0.8+0.02*x(1);t1,x=ode45(shier1,0 15,25 2);t2,y=ode45(shier2,0 15,25 2);x1=x(:,1);x2=x(:,2);x3=x2./(x1+x2);y1=y(:,1);y2=y(:,2);y3=y2./(y1+y2);(shark1.m)plot(t1,x3,-,t2,y3,*)function dy=shier2(t,y)dy=zeros(2,1);dy(1)=y(1)*(0.9-0.1*y(2);dy(2)=y(2)*(-0.6+0.02*y(1);Matlab编程实现结论：战争中鲨鱼的比例比战前高！e=0.1e=0.3食饵捕食系统的数学模型实验内容1、求微分方程的解析解,并画出它们的图形，y=y+2x,y(0)=1,0 x1;y+y=cos(x),y(0)=1,y(0)=0;2、求方程 y=y-2x/y,y(0)=1的数值解(0 x1)。3、Rossler微分方程组：当固定参数b=2,c=4时，试讨论随参数a由小到大变化（如a(0,0.65)而方程解的变化情况，并且画出空间曲线图形，观察空间曲线是否形成混沌状？实验内容返 回 一横截面积为常数A，高为H的水池内盛满水，由池底一横截面积为B的小孔放水。设水从小 孔流出的速度为v=(2gh)0.5，求在任意时刻的 水面高度和将水放空所需的时间。时间t 高度h。(要求，建立模型 选择适的参数值并给出数值解。)4、水的流出时间实验内容 两种相似的群体之间为了争夺有限的同一种食物来源和生活空间而进行生存竞争时，往往是竞争力较弱的种群灭亡，而竞争力较强的种群达到环境容许的最大数量。实验内容5、考虑相互竞争模型 假设有甲、乙两个生物种群，当它们各自生存于一个自然环境中，均服从 Logistic 规律。1、x1(t),x2(t)是两个种群的数量；2、r1,r2是它们的固有增长率；3、n1,n2是它们的最大容量；4、m2(m1)为种群乙(甲)占据甲(乙)的位置的 数量，并且 m2=x2;m1=x1。实验内容符号说明 计算x(t),y(t),画出图形及相图。解释其变化过程2）改变r1,r2,n1,n2,x0,y0,而1，2不变，计算并分析结果；若1=1.5，2=0.7，再分析结果。由此能得到什么结论。实验内容返 回结结结结 束束束束