数学实验第四次讲稿市公开课一等奖百校联赛特等奖课件.pptx

第1页 1)掌握微分方程求解三种解法：解析 法、数值解法以及图形表示解方法;2)学会使用MATLAB软件求解析解、数值解 和图形解；3）经过范例学习怎样建立微分方程模型和 分析问题思想；实验目第2页一，建立微分方程数学中有一大类问题，我们称之为微分方程，它是含有导数方程，比如 f(x,V(x),V(x)=0 是一阶微分方程，而 F(x,y,y,y(n)=0 是n阶微分方程，等等。第3页一，建立微分方程微分方程在当代科学每一个领域都有广泛应用，比如力学，运动学，电学，经济学，生物学，自动控制，化学等等，都能够看到大量利用微分方程表示事物改变规律，这也表达了微分方程主要性。在高等数学中我们对微分方程有一定认识，比如怎样求解一些简单微分方程，对于一些细心同学，可能想到对于复杂微分方程求解是一个主要问题，但 第4页一，建立微分方程 是，可能忽略一个一样主要问题，那就是怎样建立微分方程。一样是方程，其建立过程一定程度上和前一部分方程建立有相同之处，但不一样之处是微分方程建立过程中普通会包括到改变率概念（往往和时间相关）。比如运动学中速度是位移对时间改变率，加速度是速度改变率，电学中电流是电量改变率，人口学中人口增加率是人口数对时间改变率等等。第5页一，建立微分方程 下面用一些常见问题说明微分方程建模基本方法。例子1：放射性物质衰变。物理学家发觉，放射性物质衰变速度和物质量成正比。例子2：人口（生物数量）增加问题。经过观察发觉，人口增加和现有些人口数量成正比。第6页一，建立微分方程 例子3：热传导问题。牛顿定理表述热传导为：高温物体热量散失速度和温度差成百分比。例子4：物体运动问题。物体运动能够用牛顿三大定律描述。例子5：万有引力定理。开普勒发觉地球沿着以太阳为焦点椭圆 轨道运动，运动速度满足单位时间所扫过面积相等。第7页一，建立微分方程 例子6：传染病问题。医学经验表明，传染病有几个关键要素，易感人群，带病人群，而传染往往经过两种人群接触造成，不禁和易感人群和带病人群相关，还有他们接触相关。例子7：战争问题。交战双方伤亡即使比较复杂，不过大致遵照这么一些规律，比如武器杀伤力，防御力，还有双方人数对比等原因。第8页二，微分方程解法解析方法 对于一些比较简单微分方程，能够经过一些数学技巧解出，比如高等数学上接触一些方程：可分离变量方程，齐次方程，一阶线性微分方程，一些特殊二阶常系数微分方程等等。数值方法 能得到解析解微分方程毕竟只是少数，对于从实际问题中提取大量微分方程，无法得到方程解析解，数值方法能够说是实际问题中必不可少伎俩。第9页二，微分方程解法 在介绍详细方法之前，有必要了解普通方程形式和对解要求。常微分方程分成微分方程和微分方程组。比如：F(x,y,y,y(n)=0 是一个隐式高阶方程，而 是一阶微分方程组。第10页二，微分方程解法关于微分方程解，就是满足方程一个函数族（或者一条曲线族）。我们又称其为微分方程通解。用得愈加广泛是满足特定条件解，我们称其为特解。比如第11页二，微分方程解法之解析方法解析解dsolve(eqn1,eqn2,c1,var1,)微分方程组初值条件变量组注意：y Dy，y D2y 自变量名能够省略，默认变量名t。第12页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,.第13页例输入：y=dsolve(Dy=1+y2)y1=dsolve(Dy=1+y2,y(0)=1,x)输出：y=tan(t-C1)（通解，一簇曲线）y1=tan(x+1/4*pi)（特解，一条曲线）二，微分方程解法之解析方法第14页例 常系数二阶微分方程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)上述两例计算结果怎样？由此得出什么结论？例 非常系数二阶微分方程例无解析表示式！第15页x=dsolve(Dx)2+x2=1,x(0)=0)例 非线性微分方程x=sin(t)-sin(t)若欲求解某个数值解，怎样求解？t=pi/2;eval(x)二，微分方程解法之解析方法第16页输入：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软件实现软件实现返 回二，微分方程解法之解析方法第17页对于大量微分方程，只能得到其数值解，普通而言，得到解是方程一个特解近似。求微分方程数值解方法很多，比如：欧拉法，龙格库塔法等。其基本思想就是经过已知点得到函数值，并用该函数值代替一个小区间上函数导数，得到在该区间上一条直线，并用该直线作为方程特解近似。有兴趣同学能够参考微分方程数值解方面著作。二，微分方程解法之数值方法第18页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）第19页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)范例第20页 1、在解n个未知函数方程组时，x0和x均为n维向量，m-函数文件中待解方程组应以x分量形式写成.2、使用Matlab软件求数值解时，高阶微分方程必须等价地变换成一阶微分方程组.注意注意:第21页注意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;%列向量列向量第22页比如：比如：令令注意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)?返 回第23页例例2 Van der pol 方程方程：令令 y1=x(t),y2=x(t);该方程是否有解析解？范例第24页（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。第25页 蓝色曲线 y（1）；（原方程解）红色曲线 y（2）；计算结果范例第26页范例第27页例3 考虑Lorenz模型：其中参数=8/3，=10，=28解：1）编写M函数文件(lorenz.m);2)数值求解并画三维空间相平面轨线；(ltest.m)范例第28页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计算结果以下列图范例第29页图中，x1图形为实线(蓝），x2图形为“*”线（绿），x3图形为“+”线（红）.取t0，tf=0，10若自变量区间取0，20、0，40，计算结果以下范例第30页曲线呈震荡发散状三维图形混沌状ltest.m第31页观察结果 1、该曲线包含两个“圆盘”，每一个都是由螺线形轨道组成。一些轨道几乎是垂直地离开圆盘中一个而进入另一个。2、伴随t增加，x(t)先绕一个圆盘几圈，然后“跳”到另一个圆盘中。绕第二个圆盘几圈，又跳回原来圆盘。并以这么方式继续下去，在每个圆盘上绕圈数是随机。思索：该空间曲线与初始点x0选择相关吗？第32页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。第33页三，微分方程建模实例 1追击路线问题：一艘缉私舰雷达发觉距一艘缉私舰雷达发觉距c km处有一艘走私船正以匀速处有一艘走私船正以匀速 a 沿直线行驶。沿直线行驶。缉私舰马上以最大速度缉私舰马上以最大速度 b 追赶，若用雷达追赶，若用雷达进行跟踪，保持船瞬时速度方向一直指向进行跟踪，保持船瞬时速度方向一直指向走私船，试求缉私舰追逐路线和追上时间。走私船，试求缉私舰追逐路线和追上时间。第34页0(c,0)xyDR假设1、建立坐标系；2、缉私船在(c,0)处发觉走私船在(0,0)处；3、走私船逃跑方向为y轴方向；4、t时刻，走私船抵达R(0,at)缉私舰抵达D(x,y)；1 追击路线问题第35页0(c,0)xyDRyx(0,at)化简：对 t求导化简：（1）1 追击路线问题第36页(2)结合(1)、(2)得到以下微分方程：是否存在解析解？1 追击路线问题第37页1）当a b时，k=a/b b时，显然缉私船不可能追上走私船。讨论当a=b时呢？思索考虑数值求解法1 追击路线问题第38页敌艇逃跑速度：a=60公里/小时；缉私舰追击速度：b=80公里/小时缉私舰发觉敌艇时相距：c=500公里；计算得：k=(a/b)=0.75;编程计算，并画出追击路线和追击所花费时间。假设1 追击路线问题第39页转化成一阶微分方程组第40页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 追击路线问题第41页262.7时间：其中a=60(公里/小时)(小时)1 追击路线问题缉私舰追击路线敌艇逃跑路线第42页思索1、在本问题求解过程中，假定了走私艇逃跑方向是正北方向，而初始缉私艇位置在x轴正向。假如放宽这个假定，也就是当这个夹角是任意角度时，怎样建立方程进行求解。2，假如有多个走私艇在一个位置上进行交易，而缉私艇向该方向追赶。这些走私艇向不一样方向四散逃走，问怎样安排追赶路线？1 追击路线问题返 回第43页范例：范例：2 地中海鲨鱼问题地中海鲨鱼问题 意大利生物学家Ancona曾致力于鱼类种群相互制约关系研究，他从第一次世界大战期间,地中海各港口捕捉几个鱼类捕捉量百分比资料中，发觉鲨鱼百分比有显著增加（见下表）。第44页 但为何鲨鱼百分比大幅增加呢？生物学家Ancona无法解释这个现象，于是求援于著名意大利数学家V.Volterra，希望建立一个食饵捕食系统数学模型，定量地回答这个问题.因为战争使打鱼量下降，食用鱼增加，显然鲨鱼也随之增加。战争为何使鲨鱼数量增加？是什么原因？想2 地中海鲨鱼问题地中海鲨鱼问题第45页x1(t),x2(t)分别是食饵、捕食者在t时刻数量；r1,r2是食饵、捕食者固有增加率；1是捕食者掠取食饵能力，2是食饵对捕食者供养能力；1、符号说明：食饵捕食系统数学模型 捕食者存在使食饵增加率降低，假设降低程度与捕食者数量成正比，即2、基本假设：第46页食饵对捕食者数量x2起到增加作用，其程度与食饵数量x1成正比，即：综合以上和，得到以下模型：模型(一)：不考虑人工捕捉情况食饵捕食系统数学模型第47页 该 模型反应了在没有些人工捕捉自然环境中食饵与捕食者之间制约关系，没有考虑食饵和捕食者本身阻滞作用，是Volterra提出最简单模型.给定一组详细数据，用matlab软件求解。食饵：r1=1,1=0.1,x10=25;捕食(鲨鱼)：r2=0.5,2=0.02,x20=2;食饵捕食系统数学模型第48页 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)编制程序以下：食饵捕食系统数学模型第49页求解结果：求解结果：图2 反应了x1(t)与x2(t)关系。图2图1x1(t)x2(t)食饵捕食系统数学模型第50页食饵-捕食者模型数值解曲线初值：2x2(t)最大值最大值28.4x1(t)最大值最大值 99.325周期点10.7食饵捕食系统数学模型第51页食饵-捕食者模型相图（x1,x2）p 重心相轨线是封闭曲线等价于x1(t),x2(t)是周期函数。重心坐标怎样计算？(25,10)食饵捕食系统数学模型第52页模型（二）考虑人工捕捉情况 假设人工捕捉能力系数为e，相当于食饵自然增加率由r1 降为r1-e，捕食者死亡率由r2 增为 r2+e，所以模型(一)修改为：设战前捕捉能力系数e=0.3,战争中降为e=0.1,其它参数与模型(一)参数相同。观察结果会怎样改变？第53页即计算1）当e=0.3时：2）当e=0.1时：分别求出两种情况下鲨鱼在鱼类中所占百分比。画曲线：画曲线：plot(t,p1(t),t,p2(t),*)食饵捕食系统数学模型第54页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编程实现第55页结论：战争中鲨鱼百分比比战前高！e=0.1e=0.3食饵捕食系统数学模型第56页试验内容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)。第57页3、Rossler微分方程组：当固定参数b=2,c=4时，试讨论随参数a由小到大改变（如a(0,0.65)而方程解改变情况，而且画出空间曲线图形，观察空间曲线是否形成混沌状？试验内容返 回第58页 一横截面积为常数A，高为H水池内盛满水，由池底一横截面积为B小孔放水。设水从小 孔流出速度为v=(2gh)0.5，求在任意时刻 水面高度和将水放空所需时间。时间t 高度h。(要求，建立模型 选择适参数值并给出数值解。)4、水流出时间试验内容第59页 两种相同群体之间为了争夺有限同一个食物起源和生活空间而进行生存竞争时，往往是竞争力较弱种群灭亡，而竞争力较强种群到达环境允许最大数量。试验内容5、考虑相互竞争模型 假设有甲、乙两个生物种群，当它们各自生存于一个自然环境中，均服从 Logistic 规律。第60页1、x1(t),x2(t)是两个种群数量；2、r1,r2是它们固有增加率；3、n1,n2是它们最大容量；4、m2(m1)为种群乙(甲)占据甲(乙)位置 数量，而且 m2=x2;m1=x1。试验内容符号说明第61页 计算x(t),y(t),画出图形及相图。解释其改变过程2）改变r1,r2,n1,n2,x0,y0,而1，2不变，计算并分析结果；若1=1.5，2=0.7，再分析结果。由此能得到什么结论。试验内容返 回结结结结 束束束束第62页