第五讲Matlab求解微分方程-内江师范学院数值仿真与数学实验教学

第五讲 Matlab求解微分方程教学目的：学会用MATLAB求简单微分方程的解析解、数值解，加深对微分方程概念和应用的理解；针对一些具体的问题，如追击问题，掌握利用软件求解微分方程的过程；了解微分方程模型解决问题思维方法及技巧；体会微分方程建摸的艺术性．教学重点：利用机理分析建模微分方程模型，掌握追击问题的建模方法，掌握利用MATLAB求解数值解．教学难点：利用机理分析建模微分方程模型，通过举例，结合图形以及与恰当的假设突破教学难点．1微分方程相关函数（命令）及简介因为没有一种算法可以有效地解决所有的ODE 问题，为此，Matlab 提供了多种求解器Solver，对于不同的ODE 问题，采用不同的Solver．阶常微分方程(组)的初值问题的解的 Matlab 的常用程序，其中：ode23 采用龙格-库塔2 阶算法，用3 阶公式作误差估计来调节步长，具有低等的精度．ode45 则采用龙格-库塔4 阶算法，用5 阶公式作误差估计来调节步长，具有中等的精度．2 求解微分方程的一些例子2.1 几个可以直接用 Matlab 求微分方程精确解的例子：例1：求解微分方程22x xe xy dxdy-=+，并加以验证． 求解本问题的Matlab 程序为：syms x y %line1 y=dsolve('Dy+2*x*y=x*exp(-x^2)','x') %line2 diff(y,x)+2 *x*y-x*exp(-x^2) %line3 simplify(diff(y,x)+2*x*y-x*exp(-x^2)) %line4说明：(1) 行line1是用命令定义x,y 为符号变量．这里可以不写，但为确保正确性，建议写上；(2) 行line2是用命令求出的微分方程的解：1/2*exp(-x^2)*x^2+exp(-x^2)*C1(3) 行line3使用所求得的解．这里是将解代入原微分方程，结果应该为0，但这里给出：-x^3*exp(-x^2)-2*x*exp(-x^2)*C1+2*x*(1/2*exp(-x^2)*x^2+exp(-x^2)*C1)(4) 行line4 用 simplify() 函数对上式进行化简，结果为 0， 表明)(x y y =的确是微分方程的解．例2：求微分方程0'=-+x e y xy 在初始条件e y 2)1(=下的特解，并画出解函数的图形．求解本问题的 Matlab 程序为： syms x yy=dsolve('x*Dy+y-exp(x)=0','y(1)=2*exp(1)','x') ezplot(y)微分方程的特解为：y=1/x*exp(x)+1/x* exp (1) (Matlab 格式)，即xe e y x+=，此函数的图形如图 1：图1 y 关于x 的函数图象2.2 用ode23、ode45等求解非刚性的标准形式的一阶常微分方程(组)的初值问题的数值解(近似解)．例３：求解微分方程初值问题⎪⎩⎪⎨⎧=++-=1)0(2222y xx y dx dy 的数值解，求解范围为区间[0, 0.5]．fun=inline('-2*y+2*x^2+2*x','x','y'); [x,y]=ode23(fun,[0,0.5],1); x'; y'; plot(x,y,'o-') >> x' ans =0.0000 0.0400 0.0900 0.1400 0.1900 0.2400 0.2900 0.3400 0.3900 0.4400 0.4900 0.5000>> y' ans =1.0000 0.9247 0.8434 0.7754 0.7199 0.6764 0.6440 0.6222 0.6105 0.6084 0.6154 0.6179 图形结果为图 2．图2 y 关于x 的函数图像3 常微分在实际中的应用 3.1 导弹追踪问题1、符号说明，w ，乙舰的速率恒为v 0；设时刻t 乙舰的坐标为((),())X t Y t ，导弹的坐标为((),())x t y t ；当零时刻，((0),(0))(1,0)X Y =，((0),(0))(0,0)x y =．建立微分方程模型．202,01(0)0,'(0)0v d yk x dx w y y ⎧⎪⎪= =<<⎨⎪⎪==⎩由微分方程模型解得11(1)(1)11(),12(1)2(1)2(1)2(1)k k x x y x k k k k k +-+--=-+-≠+-+-++代入题设的数据1/5k =，得到导弹的运行轨迹为4655555(1)(1)81224y x x =--+-+当1=x 时245=y ，即当乙舰航行到点)245,1(处时被导弹击中. 被击中时间为:00245v v y t ==. 若v 0=1, 则在t =0.21处被击中. 利用MALAB 作图如图３．clear, x=0:0.01:1; y=-5*(1-x).^(4/5)/8+5*(1-x).^(6/5)/12+5/24; plot(x,y,'*')图３导弹运行轨迹（解析法） 图４两种方法对比的导弹运行轨迹2、数值方法求解．设导弹速率恒为w ，则得到参数方程为⎪⎪⎩⎪⎪⎨⎧--+-=--+-=)()()()()()(2222y Y y Y x X wdt dy x X y Y x X w dt dx因乙舰以速度0v 沿直线1x =运动，设01v =，5,1,w X Y t ===，因此导弹运动轨迹的参数方程为：⎪⎪⎪⎩⎪⎪⎪⎨⎧==--+-=--+-=0)0(,0)0()()()1(5)1()()1(52222y x y t y t x dtdyx y t x dt dxMATLAB 求解数值解程序如下，结果见图４ t0=0,tf=0.21;[t,y]=ode45('eq2',[t0 tf],[0 0]); X=1;Y=0:0.001:0.21;plot(X,Y,'-') plot(y(:,1),y(:,2),'*'),hold onx=0:0.01:1; y=-5*(1-x).^(4/5)/8+5*(1-x).^(6/5)/12+5/24; plot(x,y,'r')很明显数值计算的方法比较简单而且适用． 3.2 蚂蚁追击问题（1）建立平面直角坐标系．以正多边形的中心为原点, 设正多边形的一个顶点为起始点, 连接此点和原点作x 轴. 根据x 轴作出相应的y 轴; 选取足够小的t ∆进行采样．（2）在每一时刻t ，计算每只蚂蚁在下一时刻t t +∆时的坐标．不妨设甲追逐对象是乙，在时间t 时，甲的坐标为A 11(,)x y ，乙的坐标为B 22(,)x y ．甲在t t +∆时在'A 点(如图1), 其坐标为11(cos ,sin )x v t y v t θθ+∆+∆，其中2121cos ,sin ,x x y yd d dθθ--===. 同理，依次计算下一只蚂蚁在t t +∆时的坐标．通过间隔t ∆进行采样，得到新一轮各个蚂蚁在一个新的正多边形位置坐标．（4）重复2）步，直到d 充分小为止．（5）连接每只蚂蚁在各时刻的位置，就得到所求的轨迹．用MALAB 求解并作图，函数zhuJi(x,y)在附录一定义，如图６ t=[1:8]; s=7*exp(t.*2*pi/length(t)*i); x=real(s); y=imag(s);zhuJi(x,y)图6当蚂蚁为7只时的图形习题1. 求微分方程0sin 2')1(2=-+-x xy y x 的通解．2. 求微分方程x e y y y x sin 5'2''=+-的通解．3. 求微分方程组11,(,)A x y 22,(,)B x y 图1 在采样时间内，相连蚂蚁追击⎪⎪⎩⎪⎪⎨⎧=-+=++00y x dtdy y x dtdx在初始条件0|,1|00====t t y x 下的特解，并画出解函数()y f x =的图形．4. 分别用 ode23、ode45 求上述第 3 题中的微分方程初值问题的数值解(近似解)，求解区间为[0,2]t ∈．利用画图来比较两种求解器之间的差异．4 参考文献:[1] Mastering Matlab 6, D. Hanselman, B. Littlefield, 清华大学出版,2002[2] 赵静等编．数学建模与数学实验（第3版）．北京：高等教育出版社．2008. [3] 姜启源编. 数学模型（第二版）．北京：高等教育出版社.1993.[4] 石勇国. 蚂蚁追击问题与等角螺线. 宜宾学院学报. 2008，(6): 23-25. [5] 张伟年，杜正东，徐冰．常微分方程．北京：高等教育出版社．2006．5 附录附录一:zhuji(x,y)的M 文件 function zhuji(x,y) clf v=1; dt=0.05;x(length(x)+1)=x(1);y(length(y)+1)=y(1); plot(x,y,'*') holdfor it=1:1000for i=1:length(x)-1d=sqrt((x(i)-x(i+1))^2+(y(i)-y(i+1))^2); x(i)=x(i)+v*dt*(x(i+1)-x(i))/d; y(i)=y(i)+v*dt*(y(i+1)-y(i))/d; endplot(x,y,'.') hold onx(length(x))=x(1); y(length(y))=y(1); end。