MATLAB数学实验.ppt

MATLAB数学实验 第六章 常微分方程 第六章 常微分方程 6.1 预备知识：常微分方程 6.2 解常微分方程的MATLAB指令 6.3 计算实验：Euler法和刚性方程组 6.4 建模实验: 导弹系统的改进 6.1 预备知识：常微分方程 数值解 自变量节点t0, t1, ….,tn, … 准确解y(tn)未知, 数值解yn作为y(tn)的近似 6.2解常微分方程MATLAB指令 将方程整理为标准形式y’ = y?2t/y。 odefun= @(t,y)y-2*t/y ; [t,y]=ode45(odefun,[0,4],1); [t,y] plot(t,y,o) ode45(odefun,0:1:4,1); [t,y]=ode45(odefun,0:1:4,1);[t,y] %M函数eg6_3fun.m function f=eg6_3fun(t,x) f(1)=-x(1)^3-x(2); f(2)=x(1)-x(2)^3; f=f(:); %注意要保证f为列向量 [t,x]=ode45(@eg6_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)); eg6_5.m(稍有修改) %根据z初始值预估：z=-1，z’=0 clear;close; sinit=bvpinit(0:4,[-1;0]) %[-1;0]是常数猜测值z=-1, z’=0 odefun=@(t,y)[y(2);-y(2)*sin(y(1))]; bcfun=@(ya,yb)[ya(1);yb(1)+2]; sol=bvp5c(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(解曲线,解点,粗略解) 6.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*odefun(t(i),y(i)); end t=t’;y=y’; 6.4 建模实验: 导弹系统的改进 例7 程序 close; fplot(‘exp(0.9*x)’, [0,10]); %模型1解析解 hold on; [t,x]=ode45(@(t,x)0.9*x*(1-x/1000),[0 10],1); %模型2数值解 plot(t,x); axis([0 10 0 1500]); hold off; 模型比较 短期预报二个模型相近, 但作为长期预报,后者较前者合理。 eg6_8fun function dy= eg6_8fun(t,y,a,b,d,theta) %编程技巧：为了防止分母为0，加了一个1e-8 dydx=(a*t*sin(theta)-y(2)+1e-8)/... (abs(d+a*t*cos(theta)-y(1))+1e-8); dy(1)=b/(1+dydx^2)^0.5; dy(2)=b/(1+dydx^(-2))^0.5; dy=dy(:); 在线算法eg6_8a 第六章 习题 ex1(1)(2)(3)(6) ex2 ex7 ex8 得下列微分方程组 初始条件为x(0)=0, y(0)=0, 对于给定的a,b,d, ?进行计算。当x(t)满足 x(t) ? d + a t cos? 出现交点，则认为已击中目标。 如果t T, 则敌舰在打击范围内，可以发射。 ? d 例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函数eg6_8fun，再用ode45求解 * * 1.微分方程的概念 常微分方程: f (t,y,y’,y’’,…,y(n))=0 微分方程组: 联系一些未知函数的一组微 分方程 2.微分方程解析求解：某些特殊的方程 (1)初等积