MATLAB数学实验-第六章-常微分方程

MATLAB数学实验,第六章常微分方程,第六章常微分方程,6.1预备知识：常微分方程6.2解常微分方程的MATLAB指令6.3计算实验：Euler法和刚性方程组6.4建模实验:导弹系统的改进,1.微分方程的概念,常微分方程:f(t,y,y,y,y(n)=0微分方程组:联系一些未知函数的一组微分方程线性常微分方程:y(n)+a1(t)y(n-1)+an-1(t)y+an(t)y=b(t)若ai(t)(i=1,n)与t无关,称为常系数的若b(t)=0,称为齐次的,6.1预备知识：常微分方程,2.初等积分法分离变量法等3.常系数线性微分方程线性常微分方程的解为一个特解和相应的齐次微分方程通解的叠加。齐次微分方程的解可用特征根法求得,例1求x+0.2x+3.92x=0的通解,解特征方程为2+0.2+3.92=0roots(10.23.92求得共轭复根+i=-0.11.9774i,通解为x(t)=Aetcos(t)+Betsin(t),4.初值问题数值解,数值解法:寻求解y(t)在一系列离散节点t0t1odefun=inline(y-2*t/y,t,y);t,y=ode45(odefun,0,4,1);t,yplot(t,y,o-)t,y=ode45(odefun,0:1:4,1);t,y,例2解微分方程y=y2t/y,y(0)=1,0t4,例3解微分方程组0t,x=ode45(eg6_3fun,030,1;0.5),编程器窗口,指令窗口,例4(1)解微分方程组0t,x=ode45(odefun,030,1;0,2),编程器窗口,指令窗口,例4（1）中的函数亦可用匿名函数表示dfun=(t,x,mu)x(2);-mu*(x(1)2-1)*x(2)-x(1)t,x=ode45(dfun,020,1;0,2);plot(t,x)或者用内联函数dfun1=inline(x(2);mu*(1-x(1)2)*x(2)-x(1),t,x,mu)t,x=ode45(t,x)dfun1(t,x,2),020,1;0);,作图,subplot(1,2,1);plot(t,x(:,1),t,x(:,2),:);subplot(1,2,2);plot(x(:,1),x(:,2);,例4(2)求解微分方程组已知当=0时,f=0,T=1,解首先引入辅助变量，t=,y1=f,y4=T,化为一阶方程组,functionf=eg6_4fun(t,y)f(1)=y(2);f(2)=y(3);f(3)=-3*y(1)*y(3)+2*y(2)2-y(4);f(4)=y(5);f(5)=-2.1*y(1)*y(5);f=f(:);,y0=0,0,0.68,1,-0.5;t,y=ode45(eg6_4fun,05,y0);plot(t,y(:,1),t,y(:,4),:);,2.边值问题解法常微分方程边值问题,sinit=bvpinit(tinit,yinit)由在粗略节点tinit的预估解yinit生成粗略解网络sinit,这里y是2维向量，表示解函数x(t)及其导函数x(t)。,求解三部曲：粗网络、解结构、数值化,其Matlab标准形式为,sol=bvp4c(odefun,bcfun,sinit)odefun是微分方程组函数bcfun为边值条件函数sinit是由bvpinit得到的粗略解网络sol是一个结构,sol.x为求解节点.sol.y是y(t)的数值解,sx=deval(sol,ti)计算由bvp4c得到的解在ti的值。,例5求解边值问题解首先改写为标准形式。令y(1)=z,y(2)=z,则方程为y(1)=y(2),y(2)=|y(1)|边界条件为ya(1)=0,yb(1)+2=0,eg6_5.m,clear;close;sinit=bvpinit(0:4,1;0)%注意sinit的域名odefun=inline(y(2);-abs(y(1),t,y);bcfun=inline(ya(1);yb(1)+2,ya,yb);sol=bvp4c(odefun,bcfun,sinit)%注意sol的域名t=linspace(0,4,101);y=deval(sol,t);plot(t,y(1,:),sol.x,sol.y(1,:),o,sinit.x,sinit.y(1,:),s)legend(解曲线,解点,粗略解),1.Euler法,Euler法：在节点处用差商近似代替导数,Euler格式,k=0,1,2,6.3计算：Euler法和刚性方程组,M函数euler.m,functiont,y=euler(odefun,tspan,y0,h)t=tspan(1):h:tspan(2);y(1)=y0;fori=1:length(t)-1,y(i+1)=y(i)+h*feval(odefun,t(i),y(i);endt=t;y=y;,odefun:f(t,y)的函数句柄或Inline函数，tspan=t0,tf表示自变量初值t0和终值tfy0:表示初值向量y0h:步长输出列向量t表示节点(t0,t1,tn)输出矩阵y表示数值解，每一列对应y的一个分量。,M函数euler.m给出Euler法计算程序使用格式为t,y=euler(odefun,tspan,y0,h),odefun=inline(y-2*t/y,t,y);t,y=euler(odefun,0,4,1,0.01),1.产品销售量的增长例7经调查发现,电饭锅销售速度与当时的销量成正比.建立一个数学模型以预测销量.,其中k为常数。解得x(t)=x0ek(t-t0)。,6.4建模实验:导弹系统的改进,解：设x(t)表示t时刻的销量,x0为初始时刻t0的销量,模型1（指数增长模型）,当k0,t时，x(t)，这对于销售初期可认为是合适的，长期显然不合适。,设x为全部需要量，那么销售速度与当时的潜在需要量(1-x/x)成正比模型2（阻滞增长模型）,设t0=0(年),x0=1(万台),x=100(万台)k=0.9(年-1万台-1)，,例7程序,close;fplot(exp(0.9*x),0,10);%模型1解析解holdon;t,x=ode45(inline(0.9*x*(1-x/1000),t,x),010,1);%模型2数值解plot(t,x);axis(01001500);holdoff;,模型比较,短期预报二个模型相近,但作为长期预报,后者较前者合理。,目前的电子系统能迅速测出敌舰的种类、位置以及敌舰行驶速度和方向，导弹自动制导系统能保证在发射后任一时刻都能对准目标。根据情报，这种敌舰能在我军舰发射导弹后T分钟作出反应并摧毁导弹。要求：改进电子导弹系统使能自动计算出敌舰是否在有效打击范围之内。,2.导弹系统的改进,设我舰发射导弹时位置在坐标原点，敌舰在x轴正向d（km）处,其行驶速度为a(km/h),方向与x轴夹角为,导弹飞行线速度b(km/h)。,易知t时刻敌舰位置为(d+atcos,atsin)。,设t时刻导弹位置为(x(t),y(t),则,为了保持对准目标，导弹轨迹切线方向应为,由上面两个方程得下列微分方程,初始条件为x(0)=0,y(0)=0,对于给定的a,b,d,进行计算。当x(t)满足x(t)d+atcos出现交点，则认为已击中目标。,如果tT,则敌舰在打击范围内，可以发射。,例8在导弹系统中设a=90km/h,b=450km/h,T=0.1h.求d,的有效范围?,解有两个极端情形容易算出。若=0,即敌