用MATLAB进行抛体运动中探讨模拟.docx

用MATLAB进行抛体运动中探讨模拟.docx

《用MATLAB进行抛体运动中探讨模拟.docx》由会员分享，可在线阅读，更多相关《用MATLAB进行抛体运动中探讨模拟.docx（14页珍藏版）》请在冰豆网上搜索。

用MATLAB进行抛体运动中探讨模拟

MATLAB在抛体运动中的探讨

[摘要]计算机在大学物理中的应用已有二十多年的历史，MATLAB语言是一种集数值计算、符号运算、可视化建模、仿真和图形处理等多种功能的高级语言。

使用MATLAB模拟物理现象为我们解决问题提供了一种新的方法，利用其方便的数值计算和作图功能，可以方便的模拟一些物理过程。

对于处理非线性问题，既能进行数值求解，又能绘制有关曲线，方便实用，基于其功能强大，界面友善，语言自然，交互性强等优点，已成为教学和科研中最基础的软件之一，利用其解决复杂的数值计算问题，可以减少工作量，节约时间，图形绘制问题，真实直观，可以加深理解，提高工作效率。

[关键词]MATLAB语言抛体运动空气阻力力学图形绘制

一、问题的提出

MATLAB自推向市场以来，得到了广泛的应用和发展，在各高等院校中已经成为线性

代数、自动控制理论、数字信号处理、时间序列分析、动态系统仿真、图像处理等诸多课程的基本教学工具，成为大学生、硕士生、博士生必须掌握的基本技能，尤其在自动控制理论，是最具影响力、最有活力的软件。

MATLAB提供了强大的科学运算、灵活的程序设计流程、高质量的图形可视化与界面设计、便捷的与其它程序和语言接口的功能。

对于抛体运动问题，通常需要联立方程组，以及模拟它的路径，运动过程中不同的时间对应不同的位置，利用数学去计算很繁琐，手工绘图误差大，利用MATLAB可以很好地解决数值计算，模拟抛体运动的路径。

二、抛体运动的介绍

抛体运动：

将物体以一定的初速度向空中抛出，仅在重力作用下物体所作的运动，它的初速度不为零，可分为平抛运动和斜抛运动。

物理上提出的“抛体运动”是一种理想化的模型，即把物体看成质点，抛出后只考虑重力作用，忽略空气阻力。

抛体运动加速度恒为重力加速度，相等的时间速度变化量相等，并且速度变化的方向始终是竖直向下的。

一般的处理方法是将其分解为水平方向和竖直方向，平抛运动水平方向是匀速直线运动，竖直方向是自由落体运动，斜抛运动水平方向是匀速直线运动，竖直方向是竖直上抛运动，在任意方向上分解有正交分解和非正交分解两种情加速度及位移等进行相应分析。

无论怎样分解，都必须把运动的独立性和独立作用原理结合进行系统分解，即将初速度、受力情、加速度及位移等进行相应分析。

斜抛运动：

水平方向速度

（1）

竖直方向速度

（2）

水平方向位移（3）

竖直方向位移（4）

平抛运动：

水平方向速度（5）

竖直方向速度（6）

水平方向位移（7）

竖直方向位移（8）

合速度（9）

合速度方向与水平夹角：

（10）

合位移（11）

位移方向与水平夹角：

（12）

三、抛体运动的分析

1、斜抛运动的理论分析

忽略空气阻力情况下的抛射体运动是普通物理学中的一个常见问题，在高中物理教材中已有涉及，解决该问题的方法较多，分析的角度也有不同，运用的数学方法也是从初等数学到高等数学而不断深入，对该问题的分析往往是通过运用各种力学原理，推导出该运动的射程，飞行高度，飞行时间以及飞行路径曲线形状等公式，但其数值求解过程比较复杂，因此，对不具备较好高等数学基础的同学来说是较困难的。

设某一抛射体的初速度为，抛射角为，将其运动在X,Y轴上进行正交分解，水平方向速度（13）

竖直方向（14）

质点的坐标是（15）

（16）

从上两式消去，便得质点的轨迹运动方程（17）

抛射体能达到的最大高度为（18）

其到达最大高度所需时间为（19）

空中飞行时间为（20）

抛射体的最大射程为（21）

它跟初速度和抛射角有关，在抛射角不变的情况下，射程与成正比，所以射程随初速度的增大而增大。

在初速度不变的情况下，随着抛射角的增大，射程也增大，当度时，，射程达到最大值，以后随着抛射角的增大，射程减小。

利用MATLAB的绘图功能，可以更直观的体现上述结论。

x=linspace（0,pi/2,100）;%产生行向量发射角

g=10;%重力加速度

v1=10;%初速度取10

v3=20;%初速度取20

v4=25;%初速度取25

y1=v1^2*sin（2*x）/g;%初速度为10下的射程

y2=v2^2*sin（2*x）/g;%初速度为15下的射程

y3=v3^2*sin（2*x）/g;%初速度为20下的射程

y4=v4^2*sin（2*x）/g;%初速度为25下的射程

subplot（2,2,1）;%选择2*2个区的一号区

plot（x,y1）;%输出初速度为10下的射程曲线

title（'v0=10'）;%加图形标题

text（pi/4,10,'射程为10'）;%在最大射程处加图形说明

subplot（2,2,2）;%选择2*2个区的二号区

plot（x,y2）;%输出初速度为15下的射程曲线

title（'v0=15'）;%加图形标题

text（pi/4,22.5,'射程为22.5'）;%在最大射程处加图形说明

subplot（2,2,3）;%选择2*2个区的三号区

plot（x,y3）;%输出初速度为20下的射程曲线

title（'v0=20'）;%加图形标题

text（pi/4,40,'射程为40'）;%在最大射程处加图形说明

subplot（2,2,4）;%选择2*2个区的四号区

plot（x,y4）;%输出初速度为25下的射程曲线

title（'v0=25'）;%加图形标题

text（pi/4,62.5,'射程为62.5'）;%在最大射程处加图形说明

程序运行结果如图1所示。

2、斜抛运动解决实际问题

求解最大飞行路径所对应的抛射角问题（空气阻力忽略不计），如图2所示，X,Y坐标轴分别代表抛射体的射程与射高，在处，设在某一微小时段抛射体的路径变量

图1射程与抛射角、初速度的关系

为，其对应的水平及竖直方向的变量为与，

则（22）

设射程为R，则飞行路径长度（23）

根据前面的推论，（24）

其中为抛射的初始速度，为抛射角，

根据运动学原理，有（25）

（26）

从（24）、（25）中消除，我们可得到该运动的抛物线方程：

（27）

从（24）中可知，为求解L,先得求出，因此在（4）式两边同时对求导，得：

（28）

将（27）代入式（24），等式两边同时积分，便得到了飞行路径长度与抛射角之间的关系：

（29）

根据式（28），为求得L的最大值，将（28）两边同时对求导（30）

令，可得到最大飞行路径所对应的抛射角的大小，但解此方程是比较困难的。

为此，我们采用MATLAB的函数运算功能来解决这一问题。

程序如下，设其中的抛射初速度，。

x=（0:

pi/100:

pi/2）;%产生行向量x

y1=（sin（x）+（cos（x）.*cos（x））.*log（1+sin（x））./cos（x））*100/9.8;

%飞行路径长度与抛射角之间的函数关系

y2=cos（x）.*（1-sin（x）.*log（（1+sin（x））./cos（x）））*200/9.8;

%飞行路径对抛射角的一阶导数的函数关系

m=（sin（pi/6）+（cos（pi/6）*cos（pi/6））*log（1+sin（pi/6））/cos（pi/6））*100/9.8;

%抛射角取某一特定值时飞行路径值

n=cos（pi/3）*（1-sin（pi/3）*log（（1+sin（pi/3））/cos（pi/3）））*200/9.8;

%抛射角取某一特定值时飞行路径一阶导的值

plot（x,y1,'b:

'）;%输出飞行路径长度与抛射角之间的函数表达式

holdon;%设置图形保持状态

plot（x,y2,'k'）;%输出飞行路径对抛射角的一阶导数的函数表达系

holdoff;%关闭图形保持

text（pi/6,m,'y1'）;%在指定位置添加图例说明

text（pi/3,n,'y2'）;%在指定位置添加图列说明

grid;%网格线控制

运行结果如图2所示。

图2给出了飞行路径随抛射角的变化曲线及飞行路径曲线的斜度，从图中可以得到，当（弧度）时，即度时，飞行路径最大，

此时（31）

我们知道，在不考虑空气阻力的情况下，当抛射角度时，其射程最远，但此时其飞行路径并不是最远，而是当抛射角度时，其飞行路径最远，且其长度约为，实际上，由于空气阻力的存在，抛射体在空中是沿导弹曲线（弹头飞行时其重心所经过的路线）飞行的，它与抛物线不同，它的升弧与降弧不对称，在重力与空气阻力的共同影响下，弹道形成不均等的圆弧，升弧较长而直伸，降弧较短而弯曲.斜抛射出的炮弹的射程和射高都没有按抛体计算得到的值那么大，路线也不是理想曲线。

图2抛射角与飞行路径及其一阶导数曲线图

物体在空气中受到的阻力，与物体运动速度大小有密切联系，速度越小，越接近理想情况，当物体速度低于200米每秒时，阻力与物体速度大小的平方成正比，速度介于400至600米每秒之间时，空气阻力与速度大小的三次方成正比，在速度很大的情况下，阻力与速度大小的高次方成正比。

3、抛射角为90度的特殊抛体运动

一弹性小球，初始高度h=10m,向上初速度v0=15米每秒，与地面碰撞的速度衰减系数k=0.8，试计算任意时刻球的位置和速度。

高度与时间的关系：

，（32）

速度与时间关系：

（33）

对等式两边积分，有，（34）

，（35）

由此可得数学方程：

第一次落地前：

（36）

（37）

（38）

第二次落地前：

（39）

（40）

（41）

（42）

第三次落地前：

（43）

（44）

（45）

（46）

......

第n次落地前：

（47）

（48）

（49）

（50）

如用手工进行计算，计算量极大，利用MATLAB编程如下：

v0=15;%初速度

h=10;%初始高度

g=-9.8;%重力加速度

k=0.8;%衰减系数

T=0;%落地时间

fort=0:

0.05:

20%产生时间的行向量

v=v0+g*（t-T）;%求速度

y=h+v0*（t-T）+g*（t-T）^2/2;%求高度

ify<=0%循环判断条件

v0=-k*v;%衰减的速度

T=t;%求球每次落地所用时间

h=0;%将高度变零

end%选择结构结束

subplot（1,2,1）;%选择1*2中的一号区

pause（0.1）;%延缓

plot（1,y,'or','MarkerSize',10,'Markerface',[1,0,0]）;

%输出求球的运动图像

title（'运动变化图'）;%图形名称

axis（[0,2,0,25]）;%坐标控制

subplot（2,2,2）;%选择2*2中的二号区

axis（[0,20,-25,30]）;%坐标控制

gridon;%不画网格线

plot（t,v,'*r','MarkerSize',2）;%画球的速度曲线

xlabel（'时间t'）;%坐标轴说明

ylabel（'速度v'）;%坐标轴说明

title（'速度变化趋势图'）;%图形名称

holdon;%设置图形保持状态

subplot（2,2,4）;%选择2*2中的四号区

axis（[0,20,0,25]）;%坐标控制

gridon;%不加网格线

plot（t,y,'*b','MarkerSize',2）;%画球的位置曲线

xlabel（'时间t'）;%坐标轴说明

ylabel（'高度y'）;%坐标轴说明

title（'位置变化图'）;%图形名称

gridon%不加网格线

holdon%设置图形保持状态