《数学实验》第五讲.ppt

1、MATLAB数学实验,第五章 常微分方程,第五章 常微分方程,5.1 预备知识：常微分方程 5.2 解常微分方程的MATLAB指令 5.3 计算实验：Euler法和刚性方程组 5.4 建模实验: 导弹系统的改进,本章主要MATLAB指令,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,称为齐次的,5.1 预备知识：常微分方

2、程,2.初等积分法 分离变量法等 3.常系数线性微分方程 线性常微分方程的解为一个特解和相应 的齐次微分方程通解的叠加。 齐次微分方程的解可用特征根法求得,例1 求x+ 0.2 x+3.92x = 0的通解,解 特征方程为2 + 0.2 +3.92=0 roots(1 0.2 3.92 求得共轭复根 +i=-0.11.9774i, 通解为 x(t) = Aet cos(t) +Bet sin(t),4.初值问题数值解,数值解法:寻求解y(t)在一系列离散节点 t0 t1 tn tf 上的近似值yk (k=0,1,n)。 hk = tk+1 tk 为步长，通常取为常量h 。,其中,数值解,准确解

3、y(xn)未知, 数值解yn作为y(xn)的近似,高阶常微分方程初值问题可以化为一阶常微分方程组，已给一个n阶方程 y(n) = f (t, y, y,y(n-1) 设y1 = y, y2 = y, , y n= y(n-1), 化为一阶方程组,y0： 表示初值向量y0； t： 表示节点列向量 (t0 , t1 , , tn)T； y： 数值解矩阵,每一列对应y的一个分量 若无输出参数，则作出图形,初值问题求解 常用格式t,y = ode45(odefun, tspan, y0),odefun：表示f (t, y)的函数句柄或Inline函数 t是标量，y是标量或向量；,tspan： 若为t0

4、, tf, 表示自变量初值t0和终值tf 若为t0,t1, tn, 表示输出节点列向量,5.2解常微分方程MATLAB指令,完整格式 t,y = ode45(odefun, tspan, y0,options, p1,p2,) options为计算参数(如精度要求)设 置，默认可用空矩阵表示； p1,p2, 为附加传递参数，这时 odefun的表示为f (t, y, p1,p2, ), odefun= inline(y-2*t/y,t,y); t,y=ode45(odefun,0,4,1); t,y plot(t,y,o-) t,y=ode45(odefun,0:1:4,1);t,y,例2 解

5、微分方程y = y2t/y, y(0)=1, 0t4,例3 解微分方程组 0t30,%M函数eg5_3fun.m function f=eg5_3fun(t,x) f(1)=-x(1)3-x(2); f(2)=x(1)-x(2)3; f=f(:);, t,x=ode45(eg5_3fun,0 30,1;0.5),编程器窗口,指令窗口,作图, subplot(1,2,1);plot(t,x(:,1),t,x(:,2),:); subplot(1,2,2);plot(x(:,1),x(:,2);,例4 求解微分方程组 已知当=0时, f=0, T=1,解 首先引入辅助变量， t=, y1=f, ,

6、 , y4=T, 化为一阶方程组,function f=eg5_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(eg5_4fun,0 5,y0); plot(t,y(:,1),t,y(:,4),:);,2. 边值问题解法 常微分方程边值问题,sinit=bvpinit(tinit,yinit) 由在粗略节点tinit的预估解yinit生成粗略解网络sinit,这里y是2

7、维向量，表示解函数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

8、)+2=0,eg5_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

9、法：在节点处用差商近似代替导数,Euler格式,k=0,1,2,5.3 计算：Euler法和刚性方程组,M函数euler.m,function t, y = euler(odefun, tspan, y0, h) t=tspan(1):h:tspan(2); y(1)=y0; for i=1:length(t)-1, y(i+1)=y(i)+h*feval(odefun,t(i),y(i); end t=t;y=y;,odefun: f(t, y)的函数句柄或Inline函数，tspan=t0, tf 表示自变量初值t0和终值tf y0: 表示初值向量y0 h : 步长 输出列向量t 表示节点

10、 (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) = x0 e k(t-t0)。,5.4 建模实验: 导弹系统的改进,解：设x(t)表示t时刻的销量, x0为初始时刻t0 的销量, 模型

11、1（指数增长模型）,当k 0, 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解析解 hold on; t,x=ode45(inline(0.9*x*(1-x/1000),t,x),0 10,1); %模型2数值解 plot(t,x); axis(0 10 0 1500); hold o

12、ff;,模型比较,短期预报二个模型相近, 但作为长期预报,后者较前者合理。,目前的电子系统能迅速测出敌舰的种类、位置以及敌舰行驶速度和方向， 导弹自动制导系统能保证在发射后任一时刻都能对准目标。 根据情报，这种敌舰能在我军舰发射导弹后T分钟作出反应并摧毁导弹。 要求：改进电子导弹系统使能自动计算出敌舰是否在有效打击范围之内。,2. 导弹系统的改进,设我舰发射导弹时位置在坐标原点，敌舰在x轴正向d（km）处, 其行驶速度为 a (km/h), 方向与x轴夹角为 , 导弹飞行线速度b(km/h) 。,易知 t 时刻敌舰位置为 (d+atcos,atsin)。,设t 时刻导弹位置为(x(t),y(t

13、) , 则,为了保持对准目标，导弹轨迹切线方向应为,由上面两个方程得下列微分方程,初始条件为x(0)=0, y(0)=0, 对于给定的a,b,d, 进行计算。当x(t)满足 x(t) d + a t cos 出现交点，则认为已击中目标。,如果t T, 则敌舰在打击范围内，可以发射。,例8 在导弹系统中设 a=90km/h, b=450km/h, T=0.1h. 求d, 的有效范围?,解 有两个极端情形容易算出。 若 =0, 即敌舰正好背向行驶, 即x轴正向。 那么导弹直线飞行, 击中时间 t=d/(b-a)T 得d=T(b-a)=36km。 若 =, 即迎面驶来, 类似地有d=T(a+b)=54km 一般地, 有36d54。,在线算法：对于测定的d 和，可用,计算出t。如d=50, =/2，写出M函数eg5_8fun，再用od