常微分方程的求解与定性分析

学生实验报告（4）
一、实验综述
. 归纳和学习求解常微分方程(组)的基本原理和方法；
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．求微分方程⎪⎩
⎪⎨⎧====-+]100[0)0(;0)0(01.03t u
u u u u 的数值解，要求编写求解程序。

首先建立函数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 微分方程组：
⎪⎩
⎪⎨⎧-+=+=--=)('''c x z b z ay
x y z y x 当固定参数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 ），结果又如何？此时炮弹的运行轨迹如何？
（1）
建立模型
x=200cosθ*t
y=200sinθ*t-1/2*9.8t2
函数myfun6.m文件如下：
function y=myfun6(x)
y=360*tan(x)-4.9*(360/200/cos(x)).^2-160;
建立Newton2.m函数如下：
function y=newton2(x0,x1,n,tol)
x(1)=x0;
x(2)=x1;
b=1;
i=2;
while(abs(b)>eps*x(i))
x(i+1)=x(i)-myfun6(x(i))*(x(i)-x(i-1))/(myfun6(x(i))-myfun6(x(i-1)));
b=x(i+1)-x(i);
i=i+1;
if(i>n) error('n is full');
end
end
disp(i-2);
y=x(i);
运行结果如下：
（2）
建立模型如下：
代入初始条件可以得出x=-10*200cosθ*exp(-0.1t)+10*200cosθ
建立myfun7函数如下：function y=myfun7(x)
y=200*sin(x)*(-10*log(1-360/2000/cos(x)))-4.9*((-10*log(1-360/2000/cos(x))).^2)-160 建立newton2函数如下：
function y=newton2(x0,x1,n,tol)
x(1)=x0;
x(2)=x1;
b=1;
i=2;
while(abs(b)>eps*x(i))
x(i+1)=x(i)-myfun7(x(i))*(x(i)-x(i-1))/(myfun7(x(i))-myfun7(x(i-1)));
b=x(i+1)-x(i);
i=i+1;
if(i>n) error('n is full');
end
end
disp(i-2);
y=x(i);
运行结果如下：
错后进行调制、再进行运行得到结果，结束。

心得体会：通过自己的一次次锻炼，总算对matlab的线性规划有所了解了，虽然操作中有很多错误，尤其是模型没有建好所带来的错误，在一次次的修改调试中不断提升自己的能力，希望自己进一步加强学习Matlab。

四、指导教师评语及成绩：
成绩：指导教师签名：
批阅日期：。