阻尼斜抛运动ode解法

练习七阻尼斜抛运动1.实验题目研究受空气阻力作用的抛体运动，设抛体质量为m，初速度为，所受空气阻力的大小R与速率v的n次方成正比，即，其中b是阻尼系数。2．实验目的和要求⑴画出三种情况下即当没有空气阻力，空气阻力分别与速度一次方，二次方成正比时，抛体的运动轨迹图及水平速度随时间变化的图形。⑵求出三种情况的抛体轨迹的最高点，到达最高点所需的时间及所计算的抛体轨迹的终点的水平速度⑶学习使用解常微分方程的指令ode45，尤其注意如何编写一个单独的函数文件并向函数传递参数b,n的值⑷学习分区作图的方法，掌握在图形上加上各种标注文字的方法3．解题分析将抛体视为质点，根据牛顿运动定律，抛体的运动微分方程可统一写为(2.1.1)以抛出点为原点建立直角坐标系Oxy,Ox延水平方向，Oy垂直向上。于是由方程(2.1.1)可得两个投影方程(2.1.2)当b=0即空气阻力为零时，抛体做匀变速运动，方程组(2.1.2)可以用代数方法求解当时，如果n=1，则空气阻力与速率的一次方成正比（它适用于低速情况即v约为），方程组(2.1.2)可以求出解析解。做法如下：0vnRbv212nndrmmgbvmgbvdtv122222ndxbxxydtmvvv122222ndybgyxydtmvvv0b10.1ms将方程组(2.1.2)改写为(2.1.3)取初始条件为t=0时，，，将上式分离变量积分，得(2.1.4)取初始条件为t=0时,x=y=0,再积分一次，得到抛体的运动学方程为(2.1.5)消去(2.1.5)式中的时间t，得到抛体的轨道方程(2.1.6)由(2.1.4)，(2.1.5)式可知，当时，，，它的物理图像为：在水平方向抛体受空气阻力的水平分力作用，水平速度不断减小而趋于零；在竖直方向抛体受重力和空气阻力的竖直分力作用，上升阶段竖直速度逐渐减小直至为零，下降阶段的初期重力大于阻力分力使竖直速率增加，速率增加则阻力22xxdxdvbvdtdtm22yydydvbgvdtdtm0xxvv0yyvv0btmxxdxvvedt0btmyydymgmgvvedtbb01btxmmvxeb2021btymmvmgmgyetbbb202000ln1yxxxvmgmgbxyxvbvbmvt0xvymgvb增大，直到重力与阻力大小相等且方向相反，速率达到一极限值，并以匀速下降，此时轨道趋近于渐近线由式(2.1.4)和(2.1.5)知，y=0时的x值为水平射程；时的x,y值为轨迹最高点位置；轨迹最高点位置亦可由(2.1.6)式，按求出当n=2时，求方程组(2.1.2)的解析解非常困难。下面对三种情况分别计算数值解，初始条件都取为t=0时，x=0,,y=0,设，，，，将方程组(2.1.2)化为4个一阶微分方程(2.1.7)再将方程组(2.1.7)写成一个单独的函数文件，而b和p=n-1作为函数文件中的参数，然后编写主程序文件znxpyd.m。编程的基本思路是：根据三种情况设置参数b和n的三个值，用for循环对三种情况重复解三遍常微分方程，解微分方程的指令是ode45，每次解微分方程都用题目给定的相同的初始条件，但要将不同的b和n的值传递给函数文件。解微分方程的时间范围取为0到10，步长为0.01。在图形窗口用分区作图指令画了两幅图，一幅是抛体的运动轨迹；另一幅是抛体的水平速度分量随时间变化，作图用彗星轨迹的指令comet。为了比较，1mgvb1v0xmvxb0ydyvdt0dydx3dxdt5dydt1yx2dxydt3yy4dyydt12dyydt12222224ndybyyydtm34dyydt14222424ndybgyyydtm设置了坐标轴的范围，并用holdon将三种运动轨迹画在一幅图内。为了便于区分三条轨迹，在图中加注了文字说明。图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，0.4060，6.1888。从这个结果和图形都可以看出，在给定的时间范围内，第一种情况的水平速度保持不变；第二种情况的水平速度虽然在不断地下降，但尚未达到极限值零；而第三种情况地水平速度已经基本达到极限值零。如果适当增加解微分方程的时间范围，也能使第二种情况水平速度达到极限值，这与前面解析分析的结果一致。注意由于计算的时间范围较大，抛体上升的高度较小，所以上升段的曲线不明显，如果缩短计算的时间，就可以更明显地表现上升段的曲线。4．思考题⑴将解微分方程的时间范围缩短为0～1s，适当改变坐标轴的设定范围，重新画图，结果有什么不同？⑵画出抛体的垂直速度随时间变化的图形，以观察它达到极限值的情况⑶将作图指令comet改为plot,结果有什么不同？5．参考程序主程序的文件名是znxpyd.mm=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';figurefori=1:3[t,y]=ode45('znxpydfun',[0:0.01:10],[0,3,0,5],[],b(i),p(i),m);610H{i}=max(y(:,3))T{i}=t(find(y(:,3)==H{i}))vx0{i}=y(end,2)subplot(2,1,1)axis([06-702]);holdonxlabel('x')ylabel('y')plot(y(:,1),y(:,3));subplot(2,1,2)axis([01004])holdonxlabel('t')ylabel('dx/dt')text(px(i),py(i),strdd{i});plot(t,y(:,2))end函数文件是一个独立的文件，文件名为znxpydfun.mfunctionydot=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)];