微分方程组Matlab的解析解数值解.ppt

单击此处编辑母版标题样式,单击此处编辑母版文本样式,第二级,第三级,第四级,第五级,*,数学实验,Experiments in Mathematics,重庆邮电学院基础数学教学部,微 分 方 程,实验目的,实验内容,MATLAB,2,、学会用,Matlab,求微分方程的数值解,.,实验软件,1,、学会用,Matlab,求简单微分方程的解析解,.,1,、,求简单微分方程的解析解,.,4,、实验作业,.,2,、求微分方程的数值解,.,3,、数学建模实例,求微分方程的数值解,（,一）常微分方程数值解的定义,（二）建立数值解法的一些途径,（,三）用,Matlab,软件求常微分方程的数值解,返 回,1,目标跟踪问题一：导弹追踪问题,2,、目标跟踪问题二：慢跑者与狗,3,、地中海鲨鱼问题,返 回,数学建模实例,微分方程的解析解,求微分方程（组）的解析解命令,:,dsolve,(,方程,1,方程,2,方程,n,初始条件,自变量,),To,Matlab,（,ff1,）,结 果：,u=,tg(t-c,),解,输入命令,:y=dsolve(D2y+4*Dy+29*y=0,y(0)=0,Dy(0)=15,x),结 果 为,:y=3e,-2x,sin,（,5x,）,To,Matlab,（,ff2,）,解,输入命令：,x,y,z=,dsolve(Dx,=2*x-3*y+3*,z,Dy,=4*x-5*y+3z,Dz,=4*x-4*y+2*z,t),；,x=simple(x)%,将,x,化简,y=simple(y),z=simple(z),结 果 为：,x=(c,1,-c,2,+c,3,+c,2,e,-3t,-c,3,e,-3t,)e,2t,y=-c,1,e,-4t,+c,2,e,-4t,+c,2,e,-3t,-c,3,e,-3t,+(c,1,-c,2,+c,3,)e,2t,z=(-c,1,e,-4t,+c,2,e,-4t,+c,1,-c,2,+c,3,)e,2t,To,Matlab,（,ff3,）,返 回,微分方程的数值解,（,一）常微分方程数值解的定义,在生产和科研中所处理的微分方程往往很复杂且大多得不出一般解。而在实际上对初值问题，一般是要求得到解在若干个点上满足规定精确度的近似值，或者得到一个满足精确度要求的便于计算的表达式。,因此，研究常微分方程的数值解法是十分必要的,。,返 回,（,二）建立数值解法的一些途径,1,、用差商代替导数,若步长,h,较小，则有,故有公式：,此即,欧拉法,。,2,、使用数值积分,对方程,y=f(x,y),两边由,x,i,到,x,i+1,积分，并利用梯形公式，有：,实际应用时，与欧拉公式结合使用：,此即,改进的欧拉法,。,故有公式：,3,、使用泰勒公式,以此方法为基础，有,龙格,-,库塔法,、,线性多步法,等方法。,4,、数值公式的精度,当一个数值公式的截断误差可表示为,O,（,h,k+1,）,时（,k,为正整数，,h,为步长），称它是一个,k,阶公式,。,k,越大，则数值公式的精度越高。,欧拉法是一阶公式，改进的欧拉法是二阶公式。,龙格,-,库塔法有二阶公式和四阶公式。,线性多步法有四阶阿达姆斯外插公式和内插公式。,返 回,（,三）用,Matlab,软件求常微分方程的数值解,t，x=solver（f,ts,x,0,options）,ode45 ode23 ode113ode15sode23s,由待解方程写成的,m-,文件名,ts,=t,0,，,t,f,，,t,0,、,t,f,为自变量的初值和终值,函数的初值,ode23,：,组合的,2/3,阶龙格,-,库塔,-,芬尔格算法,ode45,：,运用组合的,4/5,阶龙格,-,库塔,-,芬尔格算法,自变量值,函数值,用于设定误差限,(,缺省时设定相对误差,10,-3,绝对误差,10,-6,),命令为：,options=,odeset,（,reltol,rt,abstol,at,）,rt,，,at,：,分别为设定的相对误差和绝对误差,.,1,、在解,n,个未知函数的方程组时，,x,0,和,x,均为,n,维向量，,m-,文件中的待解方程组应以,x,的分量形式写成,.,2,、使用,Matlab,软件求数值解时，高阶微分方程必须等价地变换成一阶微分方程组,.,注意,:,解,:,令,y,1,=x,，,y,2,=y,1,1,、建立,m-,文件,vdp1000.m,如下：,function,dy,=vdp1000(t,y),dy,=zeros(2,1);,dy(1)=y(2);,dy(2)=1000*(1-y(1)2)*y(2)-y(1);,2,、取,t,0,=0,，,t,f,=3000,，,输入命令：,T,Y=ode15s(vdp1000,0 3000,2 0);,plot(T,Y(:,1),-),3,、结果如图,To,Matlab,（,ff4,）,解,1,、建立,m-,文件,rigid.m,如下：,function,dy,=rigid(t,y),dy,=zeros(3,1);dy(1)=y(2)*y(3);,dy(2)=-y(1)*y(3);,dy(3)=-0.51*y(1)*y(2);,2,、,取,t,0,=0,，,t,f,=12,，,输入命令：,T,Y=ode45(rigid,0 12,0 1 1);,plot(T,Y(:,1),-,T,Y(:,2),*,T,Y(:,3),+),3,、结果如图,To,Matlab,（,ff5,）,图中，,y,1,的图形为实线，,y,2,的图形为“*”线，,y,3,的图形为“,+”,线,.,返 回,导弹追踪问题,设位于坐标原点的甲舰向位于,x,轴上点,A(1,0),处的乙舰发射导弹，导弹头始终对准乙舰,.,如果乙舰以最大的速度v,0,(,是常数,),沿平行于,y,轴的直线行驶，导弹的速度是,5v,0,，,求导弹运行的曲线方程,.,又乙舰行驶多远时，导弹将它击中？,解法一,（解析法）,由,(1),(2),消去,t,整理得模型,:,To Matlab(chase1),轨迹图见程序,chase1,解法二,(,数值解,),1.,建立,m-,文件,eq1.m,function,dy,=eq1(x,y),dy,=zeros(2,1);,dy(1)=y(2);,dy(2)=1/5*sqrt(1+y(1)2)/(1-x);,2.,取,x,0,=0,，,x,f,=0.9999,，,建立主程序,ff6.m,如下,：,x,0,=0,，x,f,=0.9999,x,y=ode15s(,eq1,x0 xf,0 0);,plot(x,y(:,1),b.,),hold on,y=0:0.01:2;,plot(1,y,b*,),结论,:,导弹大致在（,1,，,0.2,）处击中乙舰,To Matlab(ff6),令,y,1,=y,y,2,=y,1,，,将方程（,3,）化为一阶微分方程组。,解法三,(,建立参数方程求数值解,),设时刻,t,乙舰的坐标为,(X(t),，,Y(t),，,导弹的坐标为,(x(t),，,y(t).,3,因乙舰以速度,v,0,沿直线,x=1,运动，设,v,0,=1,，则,w=5,，,X=1,，,Y=t,4.,解导弹运动轨迹的参数方程,建立,m-,文件,eq2.m,如下：,function,dy,=eq2(t,y),dy,=zeros(2,1);,dy(1)=5*(1-y(1)/sqrt(1-y(1)2+(t-y(2)2);,dy(2)=5*(t-y(2)/sqrt(1-y(1)2+(t-y(2)2);,取,t,0,=0,，,t,f,=2,，,建立主程序,chase2.m,如下：,t,y=ode45(,eq2,0 2,0 0);,Y=0:0.01:2;,plot(1,Y,-,),，,hold on,plot(y(:,1),y(:,2),*,),To Matlab(chase2),5.,结果见图,1,导弹大致在（,1,，,0.2,）处击中乙舰，与前面的结论一致,.,图1,图2,返 回,在,chase2.m,中，按二分法逐步修改,t,f,，,即分别取,t,f,=1,，,0.5,0.25,直到,t,f,=0.21,时，得图,2.,结论：时刻,t=0.21,时，导弹在（,1,，,0.21,）处击中乙舰。,To Matlab(chase2),慢跑者与狗,一个慢跑者在平面上沿椭圆以恒定的速率,v=1,跑步,设椭圆方程为,:x=10+20cost,y=20+5sint.,突然有一只狗攻击他,.,这只狗从原点出发,以恒定速率,w,跑向慢跑者,狗的运动方向始终指向慢跑者,.,分别求出,w=20,w=5,时狗的运动轨迹,.,1.,模型建立,设时刻,t,慢跑者的坐标为,(X(t),，,Y(t),，,狗的坐标为,(x(t),，,y(t).,则,X=10+20cost,Y=20+15sint,狗从,(0,0),出发,与导弹追踪问题类似，建立狗的运动轨迹的参数方程,:,2.,模型求解,(1)w=20,时,建立,m-,文件,eq3.m,如下,:,function,dy,=eq3(t,y),dy,=zeros(2,1);,dy(1)=20*(10+20*cos(t)-y(1)/sqrt,(10+20*cos(t)-y(1)2+(20+15*sin(t)-y(2)2);,dy(2)=20*(20+15*sin(t)-y(2)/sqrt,(10+20*cos(t)-y(1)2+(20+15*sin(t)-y(2)2);,取,t,0,=0,，,t,f,=10,，,建立主程序,chase3.m,如下：,t0=0;tf=10;,t,y=ode45(,eq3,t0 tf,0 0);,T=0:0.1:2*pi;,X=10+20*,cos(T,);,Y=20+15*sin(T);,plot(X,Y,-,),hold on,plot(y(:,1),y(:,2),*,),在,chase3.m,，,不断修改,t,f,的值,分别取,t,f,=5,2.5,3.5,至,3.15,时,狗刚好追上慢跑者,.,To Matlab(chase3),建立,m-,文件,eq4.m,如下,:,function,dy,=eq4(t,y),dy,=zeros(2,1);,dy(1)=5*(10+20*cos(t)-y(1)/sqrt,(10+20*cos(t)-y(1)2+(20+15*sin(t)-y(2)2);,dy(2)=5*(20+15*sin(t)-y(2)/sqrt,(10+20*,cos(t,)-y(1)2+(20+15*sin(t)-y(2)2);,取,t,0,=0,，,t,f,=10,，,建立主程序,chase4.m,如下：,t0=0;tf=10;,t,y=ode45(,eq4,t0 tf,0 0);,T=0:0.1:2*pi;,X=10+20*,cos(T,);,Y=20+15*sin(T);,plot(X,Y,-,),hold on,plot(y(:,1),y(:,2),*,),在,chase3.m,，,不断修改,t,f,的值,分别取,t,f,=20,40,80,可以看出,狗永远追不上慢跑者,.,To Matlab(chase4),(2)w=5时,返 回,地中海鲨鱼问题,意大利生物学家,Ancona,曾致力于鱼类种群相互制约关系的研究，他从第一次世界大战期间,地中海各港口捕获的几种鱼类捕获量百分比的资料中，发现鲨鱼等的比例有明显增加（见下表），而供其捕食的食用鱼的百分比却明显下降,.,显然战争使捕鱼量下降，食用鱼增加，鲨鱼等也随之增加，但为何鲨鱼的比例大幅增加呢？,他无法解释这个现象，于是求助于著名的意大利数学家,V.Volterra,，,希望建立一个食饵,捕食系统的数学模型，定量地回答这个问题,.,该 模型反映了在没有人工捕获的自然环境中食饵与捕食者之间的制约关系，没有考虑食饵和捕食者自身的阻滞作用，是,Volterra,提出的最简单的模型,.,首先，建立,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);,其次，建立主程序,shark.m,如下：,t,x=ode45(shier,0 15,25 2);,plot(t,x(:,1),-,t,x(:,2),*),plot(x(:,1),x(:,2),To,Matlab(shark,),求解结果：,左图反映了,x,1,（,t,）与,x,2,（,t,）,的关系。,可以猜测：,x,1,（,t,）与,x,2,（,t,）,都是周期函数。,模型（二）考虑人工捕获,设表示捕获能力的系数为,e,，,相当于食饵的自然增长率由,r,1,降为,r,1,-e,，,捕食者的死亡率由,r,2,增为,r,2,+e,设战前捕获能力系数,e=0.3,战争中降为,e=0.1,则战前与战争中的模型分别为,:,模型求解,:,1,、分别用,m-,文件,shier1.m,和,shier2.m,定义上述两个方程,2,、建立主程序,shark1.m,求解两个方程，并画出两种情况下鲨鱼数在鱼类总数中所占比例,x,2,(t)/x,1,(t)+x,2,(t),To Matlab(shark1),实线为战前的鲨鱼比例，,“,*,”,线为战争中的鲨鱼比例,结论：战争中鲨鱼的比例比战前高！,返 回,实 验 作 业,1.,一个小孩借助长度为,a,的硬棒，拉或推某玩具,.,此小孩沿某曲线行走，计算并画出玩具的轨迹,.,2.,讨论资金积累、国民收入、与人口增长的关系,.,（,1,）若国民平均收入,x,与按人口平均资金积累,y,成正比，说明仅当总资金积累的相对增长率,k,大于人口的相对增长率,r,时，国民平均收入才是增长的,.,（,2,）作出,k(x),和,r(x),的示意图，分析人口激增会引起什么后果,.,