常微分方程的求解与定性分析实验报告

一、实验综述

1、实验目的及要求

 归纳和学习求解常微分方程(组)的基本原理和方法；

 掌握解析、数值解法，并学会用图形观察解的形态和进行解的定性分析；  熟悉MATLAB软件关于微分方程求解的各种命令；

 通过范例学习建立微分方程方面的数学模型以及求解全过程；

 通过该实验的学习，使学生掌握微分方程(组)求解方法（解析法、欧拉法、

梯度法、改进欧拉法等），对常微分方程的数值解法有一个初步了解，同时学会使用MATLAB软件求解微分方程的基本命令，学会建立微分方程方面的数学模型。这对于学生深入理解微分、积分的数学概念，掌握数学的分析思维方法，熟悉处理大量的工程计算问题的方法是十分必要的。

2、实验仪器、设备或软件

电脑 、matlab7.0

二、实验过程（实验步骤、记录、数据、分析） 实验内容：

根据实验内容和步骤，完成以下实验，要求写出实验报告（实验目的→问题→数学模型→算法与编程→计算结果→分析、检验和结论） 1．求微分方程的解析解，并画出它们的图形。 y '= y + 2 x, y (0) = 1, 0< x <1； m=dsolve('Dy=y+2*x','y(0)=1','x') ezplot(m,[0 1]) m =

3*exp(x) - 2*x – 2

3 exp(x) - 2 x - 243.532.521.5100.10.20.30.40.5x0.60.70.80.91

uu0.1u30(0)0 的数值解，要求编写求解程序。 1．求微分方程u(0)0;ut[010]function dy=vdp1000(t,y) dy=zeros(2,1); dy(1)=y(2);

dy(2)=-y(1)+0.1*y(1)^3;

[T,Y]=ode15s('vdp1000',[0 10],[0 0]); plot(T,Y(:,1),'-')

10.80.60.40.20-0.2-0.4-0.6-0.8-1012345678910

3．Rossler微分方程组：

x'yz' yxayz'bz(xc)当固定参数b=2，c=4时，试讨论随参数a由小到大变化（如 a∈(0，0.65))而方程解的变化情况，并且画出空间曲线图形，观察空间曲线是否形成混沌状？ 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.1:0.6

[t,x]=ode45('rossler',t0,[0,0,0]); 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(篮色)随t变化情况0.60.40.80.20.60-0.2-0.4-0.6-0.8-1z相图0.40.2000.5-0.5050100t150200y-1-0.50xa=0.1:

x(红色),y(绿色),z(篮色)随t变化情况0.60.40.80.20.60-0.2-0.4-0.6-0.8-1z

相图0.40.2001-0.50050100t150200y-1-0.5x0.5

a=0.2:

x(红色),y(绿色),z(篮色)随t变化情况43221.510-1-2-3-4z相图10.5020-2050100t150200y-4-50x5

a=0.3:

x(红色),y(绿色),z(篮色)随t变化情况543210-1-2-3-4-5050100t150200yz相图32.521.510.50550-5-50x

a=0.4:

x(红色),y(绿色),z(篮色)随t变化情况6相图46524z320-2105-405-50-10-5x10-6050100t150200y

a=0.5:

x(红色),y(绿色),z(篮色)随t变化情况121086420-2-4-6-8050100t150200z相图121086420505-5y-10-50x10

结果分析：从图像可以看出，当a=0时，微分方程的解（x,y,z）收敛与（0，0.5，0.5）；当a=0.1时，（x,y,z）仍收敛与（0，0.5，0.5），只是收敛速度减慢；当a=0.2时，（x,y,z）已发散，周期性变化；随着a的增大，(x,y,z)接近其极限环的速度加快，空间曲线成混沌状。

4．炮弹发射角的确定 炮弹发射视为斜抛运动，已知初始速度为200m/s，若要击中水平距离360m、在垂直距离160m的目标，当忽略空气阻力时，发射角应为多大？此时炮弹的运行轨迹如何？

要求：

(1) 建立在忽略空气阻力情况下的描述炮弹发射轨迹的数学模型； (2) 用Matlab 软件求解方程和微分方程； (3) 结合实际对解的合理性进行分析。 进一步思考：

如果要考虑水平方向的阻力，且设阻力与（水平方向）速度成正比，系数为0.1（1/s ），结果又如何？此时炮弹的运行轨迹如何？

解：

（1）忽略空气阻力时，设发射角为a，炮弹的飞行时间为t 水平方向：Vx=V0*cos a=200*cos（a） 竖直方向：Vy=V0*sin a=200*sin（a） 得到：x=Vx*t=200*cos（a）*t=360

y=Vy*t-1/2*g*t^2=200*sin（a）*t-1/2*9.8*t^2=160 得到炮弹的路程为Y=360*tan(a)-4.9*(360/200/cos(a)).^2-160

编程为：

function Y=fun2(a)

Y=360*tan(a)-4.9*(360/200/cos(a)).^2-160;

function Y=fun3(a0,a1,n,tol) a(1)=a0; a(2)=a1; b=1; i=2;

while(abs(b)>eps*a(i))

a(i+1)=a(i)-fun2(a(i))*(a(i)-a(i-1))/(fun2(a(i))-fun2(a(i-1))); b=a(i+1)-a(i); i=i+1;

if(i>n) error('n is full'); end end

disp(i-2); Y=a(i);

fun3(0.5,1,100,1e-6)

结果为：

ans =

0.4633

（3）结果合理 ，符合实际

进一步思考： 建立模型如下：

错误!未找到引用源。= -0.1dx/dt

代入初始条件可以得出x=-10*200cosθ*exp(-0.1t)+10*200cosθ 建立myfun6函数如下： function Y=fun6(a)

Y=200*sin(a)*(-10*log(1-360/2000/cos(a)))-4.9*((-10*log(1-360/2000/cos(a))).^2)-160

建立fun7函数如下：

function Y=fun7(a0,a1,n,tol) a(1)=a0; a(2)=a1; b=1; i=2;

while(abs(b)>eps*a(i))

a(i+1)=a(i)-fun6(a(i))*(a(i)-a(i-1))/(fun6(a(i))-fun6(a(i-1))); b=a(i+1)-a(i); i=i+1;

if(i>n) error('n is full'); end end

disp(i-2); Y=a(i);

输入：>> fun7(0.5,1,100,1e-6)

结果： ans =

0.4297

三、结论

1、实验结果

编程及实验结果分析如上

2、分析讨论

① 通过此次实验， 学习了微分方程求解的方法，并学会了建立微分方程的数学模型，使用求解微分方程的基本指令；

② 实验题目中的第三四道题，综合性较强，需要考虑多方面的程序，感觉较难，通过查找相关例题以及和同学讨论得以解决；

③ 微分方程的求解问题较常见，因为许多的微分方程人为解起来较难，而使用matlab可以很轻松得解答。因此掌握使用matlab求解微分方程很重要，在今后的学习中要多练习，熟练使用matlab软件。