练习七阻尼斜抛运动
1.实验题目
研究受空气阻力作用的抛体运动，设抛体质量为m，初速度为v0，所受空气
nR?bv阻力的大小R与速率v的n次方成正比，即，其中b是阻尼系数。
2．实验目的和要求
⑴画出三种情况下即当没有空气阻力，空气阻力分别与速度一次方，二次方成正比时，抛体的运动轨迹图及水平速度随时间变化的图形。
⑵求出三种情况的抛体轨迹的最高点，到达最高点所需的时间及所计算的抛体轨迹的终点的水平速度
⑶学习使用解常微分方程的指令ode45，尤其注意如何编写一个单独的函数文件并向函数传递参数b,n的值
⑷学习分区作图的方法，掌握在图形上加上各种标注文字的方法 3．解题分析
将抛体视为质点，根据牛顿运动定律，抛体的运动微分方程可统一写为
d2r?m2?mg?bvn?mg?bvn?1?dtv(2.1.1)
以抛出点为原点建立直角坐标系Oxy,Ox延水平方向，Oy垂直向上。于是由方程(2.1.1)可得两个投影方程
dxb??dt2m2?vx2?vy2??n?1????2?vx
dyb??g?dt2m2?vx2?vy2??n?1????2?vy(2.1.2)
当b=0即空气阻力为零时，抛体做匀变速运动，方程组(2.1.2)可以用代数方法求解
当b?0时，如果n=1，则空气阻力与速率的一次方成正比（它适用于低速情况即v约为0.1m?s），方程组(2.1.2)可以求出解析解。做法如下：
?1将方程组(2.1.2)改写为
d2xdvxb???vxdt2dtm
d2ydvyb???g?vydt2dtm(2.1.3)
取初始条件为t=0时，vx?vx0，vy?vy0，将上式分离变量积分，得
b?tdx?vx?vx0emdt
btdymg??mmg??vy??vy0?e??dtb?b(2.1.4) ?
取初始条件为t=0时,x=y=0,再积分一次，得到抛体的运动学方程为
b?t?mvx0?mx??1?e?b??
b?t??mvy0m2g??mgy???2??1?em??tb???b?b(2.1.5)
消去(2.1.5)式中的时间t，得到抛体的轨道方程
m2g?bx??vy0mg?y???x?ln1????b2?vx0bvx0??mvx0?(2.1.6)
由(2.1.4)，(2.1.5)式可知，当t??时，vx?0，为：
在水平方向抛体受空气阻力的水平分力作用，水平速度不断减小而趋于零；在竖直方向抛体受重力和空气阻力的竖直分力作用，上升阶段竖直速度逐渐减小直至为零，下降阶段的初期重力大于阻力分力使竖直速率增加，速率增加则阻力
vy??mgb，它的物理图像
增大，直到重力与阻力大小相等且方向相反，速率达到一极限值
x?mvx0b
vy?v1?mgb，并以v1匀速下降，此时轨道趋近于渐近线
由式(2.1.4)和(2.1.5)知，y=0时的x值为水平射程；
dy?0dt时的x,y值为
dy?0轨迹最高点位置；轨迹最高点位置亦可由(2.1.6)式，按dx求出
当n=2时，求方程组(2.1.2)的解析解非常困难。下面对三种情况分别计算数值
dxdy?3?5dtdt解，初始条件都取为t=0时， x=0, ,y=0,
y2?dxdyy4?dt，y3?y，dt，将方程组(2.1.2)化为4个一阶微分方程
设y1?x，
dy1?y2dt
dy2b22??y2?y4dtm??n?12y2
dy3?y4dt
dy4b22??g?y2?y4dtm??n?12y4(2.1.7)
再将方程组(2.1.7)写成一个单独的函数文件，而b和p=n-1作为函数文件中的参数，然后编写主程序文件znxpyd.m。编程的基本思路是：根据三种情况设置参数b和n的三个值，用for循环对三种情况重复解三遍常微分方程，解微分方程的指令是ode45，每次解微分方程都用题目给定的相同的初始条件，但要将不同的b和n的值传递给函数文件。解微分方程的时间范围取为0到10，步长为0.01。在图形窗口用分区作图指令画了两幅图，一幅是抛体的运动轨迹；另一幅是抛体的水平速度分量随时间变化，作图用彗星轨迹的指令comet。为了比较，
设置了坐标轴的范围，并用hold on将三种运动轨迹画在一幅图内。为了便于区分三条轨迹，在图中加注了文字说明。
图2.1即为程序在一个图形窗口画出的两个分区图形
图2.1 抛体的运动轨迹图和抛体水平速度随时间变化的图形
轨迹的最高点也就是函数文件中y(3)的最大值，也就是微分方程的解y的第三个分量y(:,3),可用指令max求得，每次求出的最高点的值存放在基本矩阵H中，最高点对应的时间也就是上升到最高点所需的时间，求法是用指令find求出最高点在y(3)的位置即编号，再求对应这个编号的t值，每次求出的上升到最高点所需的时间的值存放在基元矩阵T。根据定义，y(2)是水平速度，所以要求的轨迹终点的水平速度就是y(2)的最后一个元素y(end,2).最后在屏幕上显示的结果是 H =
[1.2755] [1.1948] [0.9826] T =
[0.5100] [0.4900] [0.4300] vx0 =
[3] [0.4060] [6.1888e-006]
这表明，三种情况下的抛体的最大高度分别是1.2755，1.1948，0.9826；上升时间分别是0.5100，0.4900，0.4300；所求的轨迹终点的水平速度是3，
?6?100.4060，6.1888。从这个结果和图形都可以看出，在给定的时间范围内，
第一种情况的水平速度保持不变；第二种情况的水平速度虽然在不断地下降，但尚未达到极限值零；而第三种情况地水平速度已经基本达到极限值零。如果适当增加解微分方程的时间范围，也能使第二种情况水平速度达到极限值，这与前面解析分析的结果一致。注意由于计算的时间范围较大，抛体上升的高度较小，所以上升段的曲线不明显，如果缩短计算的时间，就可以更明显地表现上升段的曲线。
4．思考题
⑴将解微分方程的时间范围缩短为0～1s，适当改变坐标轴的设定范围，重新画图，结果有什么不同？
⑵画出抛体的垂直速度随时间变化的图形，以观察它达到极限值的情况 ⑶将作图指令comet改为plot,结果有什么不同？ 5．参考程序
主程序的文件名是znxpyd.m m=1;b=[0,0.2,0.2];p=[0,0,1]; px=[4.6;4.5;4.5]; py=[3.5;1.8;0.4]; strdd{1}='无阻尼'; strdd{2}='阻力正比于v'; strdd{3}='阻力正比于v^2'; figure for i=1:3
[t,y]=ode45('znxpydfun',[0:0.01:10],[0,3,0,5],[],b(i),p(i),m);
H{i}=max(y(:,3))
T{i}=t(find(y(:,3)==H{i})) vx0{i}=y(end,2) subplot(2,1,1) axis([0 6 -70 2]); hold on xlabel('x') ylabel('y')
plot(y(:,1),y(:,3)); subplot(2,1,2) axis([0 10 0 4]) hold on xlabel('t') ylabel('dx/dt')
text(px(i),py(i),strdd{i}); plot(t,y(:,2)) end
函数文件是一个独立的文件，文件名为znxpydfun.m functionydot=znxpydfun(t,y,flag,b,p,m) ydot=[y(2);
-b/m*y(2)*(y(2).^2+y(4).^2)^(p/2); y(4);
-9.8-b/m*y(4)*(y(2).^2+y(4).^2)^(p/2)];