常微分方程求解与定性分析.docVIP

学生实验报告（4） 学生姓名 陈宋 学号 008888 同组人：无 实验项目 Matlab程序设计与作图 实验仪器台号 D02 指导教师 谢建 实验日期及节次 2011-04-18第八九十节 一、实验综述 . 归纳和学习求解常微分方程(组)的基本原理和方法； 2. 掌握解析、数值解法，并学会用图形观察解的形态和进行解的定性分析； 3. 熟悉MATLAB软件关于微分方程求解的各种命令； 4. 通过范例学习建立微分方程方面的数学模型以及求解全过程； 通过该实验的学习，使学生掌握微分方程(组)求解方法（解析法、欧拉法、梯度法、改进欧拉法等），对常微分方程的数值解法有一个初步了解，同时学会使用MATLAB软件求解微分方程的基本命令，学会建立微分方程方面的数学模型。这对于学生深入理解微分、积分的数学概念，掌握数学的分析思维方法，熟悉处理大量的工程计算问题的方法是十分必要的。 二、实验过程（实验步骤、记录、数据、分析） 1．开启MATLAB软件平台，开启MATLAB编辑窗口； 2．根据问题，建立的线性规划模型，并编写求解规划模型的M文件； 3．保存文件并运行； 4．观察运行结果(数值或图形)，并不断地改变参数设置观察运行结果； 5．根据观察到的结果和体会，写出实验报告。 1．求微分方程的解析解，并画出它们的图形。 y = y + 2 x, y (0) = 1, 0 x 1； 程序如下：由y=dsolve(Dy=y+2*x,y(0)=1,x)得出解析解y =-2*x-2+3*exp(x) 建立函数m文件：function y=myfun4(x) y=-2*x-2+3*exp(x) 画图函数为fplot(myfun4,[0,1]) 图形如下： 2．求微分方程的数值解，要求编写求解程序。 首先建立函数M文件：function dy=myfun5(t,y) dy=zeros(2,1); dy(1)=y(2); dy(2)=0.1*(y(1).^3)-y(1); 输入命令：[T,Y]=ode15s(myfun5,[0,10],[0,1]); plot(T,Y(:,1),-,T,Y(:,2),*); 图形如下： 3．Rossler微分方程组： 当固定参数b=2，c=4时，试讨论随参数a由小到大变化（如 a∈(0，0.65))而方程解的变化情况，并且画出空间曲线图形，观察空间曲线是否形成混沌状？ 首先 建立函数M文件rossler.m，在其中用x(1)表示x，用x(2)表示y，用x(3)表示z. function r=rossler(t,x) global a; global b; global c; r=[-x(2)-x(3);x(1)+a*x(2);b+x(3)*(x(1)-c)]; 主程序如下： global a; global b; global c; b=2; c=4; t0=[0,200]; for a=0:0.02:0.65 [t,x]=ode45(rossler,t0,[0,0,0]); a subplot(1,2,1); plot(t,x(:,1),r,t,x(:,2),g,t,x(:,3),b); title(x(红色),y(绿色),z(篮色)随t变化情况);xlabel(t); subplot(1,2,2); plot3(x(:,1),x(:,2),x(:,3)) title(相图);xlabel(x);ylabel(y);zlabel(z); pause end 当a=0 时，(x,y,z)收敛于（0,0.5,0.5） 当a=0.05 时，(x,y,z)仍然收敛，但收敛速度较小。 当a=0.12 时，(x,y,z)仍然收敛（准确的说，只能说具有收敛的趋势，其收敛速度已 大大降低）。 a=0.16 时，(x,y,z)已经发散。但是，(x,y,z)并不是发散于无穷大，而是周期性变化。 4．炮弹发射 炮弹发射视为斜抛运动，已知初始速度为200m/s若要击中水平距离360m、垂直距离160m的目标，当忽略空气阻力时，发射角应为多大？此时炮弹的运行轨迹如何？ 要求： 1) 建立忽略空气阻力描述炮弹发射轨迹的数学模型； 2) 用Matlab 软件求解方程和微分方程； 3) 结合实际对解的合理性进行分析。如果要考虑水平方向的阻力，且设阻力与（水平方向）速度成正比，系数为0.1（1/s ），结果又如何？此时炮弹的运行轨迹如何？function y=myfun6(x) y=360*tan(x)-4.9*(360/200/cos(x)).^2-1