《MATLAB计算平抛阻尼运动.docx》由会员分享，可在线阅读，更多相关《MATLAB计算平抛阻尼运动.docx（8页珍藏版）》请在三一办公上搜索。
1、MATLAB计算平抛阻尼运动小球在空气中平抛运动的轨迹 问题一小球在空气中作平抛运动，初速度为v0，所受的阻力与速率成正比：f = -kv，k称为阻力系数。画出小球运动的轨迹。 数学模型 如B7.1图所示，小球受到重力mg，方向向下；空气阻力f，方向与速度方向相反。根据牛顿第二定律可列出直角坐标方程 x O v0 f d2xdxd2ydy，m2=mg-k。 (7.1.1) m2=-kdtdtdtdt由于vx = dx/dt，vy = dy/dt，上式可化为 mg y B7.1图 dvdvxkk=-vx，y=g-vy。 (7.1.2) dtmdtm分离变量得 dvydvxk=-dt，=dt， v
2、xmg-kvy/m积分得 lnvx+lnCx=-kvkmt，-ln(g-y)+lnCy=t， mkmkmgkt)，vy=1-exp(-t)。 (7.1.3) mkm当t = 0时，vx = v0，vy = 0，可得Cx = -v0，Cy = -g，因此 vx=v0exp(-当t = 0时，x = 0，y = 0，积分上式可得 mv0kmgm2gkx=1-exp(-t)，y=t+2exp(-t)-1。 (7.1.4) kmkkm这是小球的运动方程，也是以时间t为参数的轨道方程。 算法方法一：用解析式。取t0 = m/k为时间单位，取s0 = m2g/k2为坐标单位，坐标方程可表示为 x = s0
3、v0*1 exp(-t*)，y = s0t* + exp(-t*) - 1。 (7.1.4*) 其中，t* = t/t0，为约化时间；v* = kv0/mg，是小球的约化初速度。 程序zqy6_3ode.m的第一部分如下。 %阻力与速度成正比的平抛运动的轨迹(用解析解) clear %清除变量 v0=input(请输入水平初速度kv0/mg:); %键盘输入水平初速度 t=0:0.1:4; %时间向量 x=v0*(1-exp(-t); %有空气阻力的横坐标 y=t+exp(-t)-1; %有空气阻力的纵坐标 fs=16; %字体大小 figure %开创图形窗口 1 plot(x,y) %画轨
4、迹 grid on %加网格 xlabel(itx/srm_0,FontSize,fs) %标记横坐标 ylabel(ity/srm_0,FontSize,fs) %标记纵坐标 title(阻力与速度成正比的平抛运动的轨迹,FontSize,fs)%标题 text(0,max(y),itkvrm_0/itmgrm=,num2str(v0),FontSize,fs)%约化初速度 axis ij equal %原点设在左上角并使坐标间隔相等 说明程序执行时，从键盘输入约化初速度，例如1等值。 方法二：用两个微分方程的数值解。微分方程组(7.1.1)式可化为 d(vy/V0)gt0kt0vyd(vx
5、/V0)ktv。 =-0x，=-d(t/t0)mV0d(t/t0)V0mV0其中，t0 = m/k，V0 = gt0 = mg/k。取t* = t/t0，vx* = vx/V0，vy* = vy/V0，可得 *dvydv*x，=-v=1-vxy。 (7.1.1*) *dtdt取v(1) = vx*，v(2) = vy*，可得两个一阶方程组 dv(1)dv(2)=-v(1)=1-v(2)。 ，dt*dt*在初始时刻，小球约化初速度为 v(1)=vxkv0*， =v0V0mg而v(2) = 0。在任意时刻，v(1)和v(2)表示约化速度。根据速度值，利用指令cumtrapz指令可求坐标。 程序zq
6、y6_3ode.m的第二部分如下。 %阻力与速度成正比的平抛运动的轨迹(求两个一阶微分方程的数值解) t0,V=ode45(zqy6_3fun1,t,v0,0);%求微分方程的数值解 x=cumtrapz(V(:,1)*t(2); %横坐标 y=cumtrapz(V(:,2)*t(2); %纵坐标 hold on %保持图像 plot(x,y,r.) %画轨迹 程序在执行时将调用一个函数zqy6_3fun1.m。 %阻力与速度成正比的平抛运动的加速度函数 function f=fun(t,v) f=-v(1); %水平加速度 1-v(2); %竖直加速度 说明用两个常微分方程的数值解可求出速度
7、，还需要求位置坐标才能画轨迹。 方法三：用四个微分方程的数值解。微分方程组(7.1.1)式可化为 2d2(x/s0)kt0d(x/s0)d2(y/s0)gt0kt0d(y/s0)，。 =-=-d(t/t0)2md(t/t0)d(t/t0)2s0md(t/t0)取t* = t/t0，x* = x/s0，y* = y/s0，由于t0 = m/k，s0 = m2g/k2，可得 2 d2x*dx*d2y*dy*=-*，*2=1-*。 (7.1.1*) dt*2dtdtdt取r(1) = x*，r(2) = y*，r(3) = dx*/dt*，r(4) = dy*/dt*，可得四个一阶方程组 dr(1)
8、dr(3)=r(3)=-r(3)； ， dt*dt*dr(2)dr(4)=r(4)=1-r(4)。 ；dt*dt*在初始时刻，小球约化位移为r(1) = r(2) = 0，初始约化速度为 dx*t0dx1dxv0*， r(3)=*=v0dts0dtmg/kdtV0r(4) = 0。 在任意时刻，r(1)和r(2)表示约化坐标，r(3)和r(4)表示约化速度。 程序zqy6_3ode.m的第三部分如下。 %阻力与速度成正比的平抛运动的轨迹(求四个一阶微分方程的数值解) t0,R=ode45(zqy6_3fun2,t,0,0,v0,0);%求微分方程的数值解 plot(R(:,1),R(:,2),
9、ko) %画轨迹 程序在执行时将调用一个函数zqy6_3fun2.m。 %阻力与速度成正比的平抛运动的函数 function f=fun(t,r) f= r(3); %水平速度 r(4); %竖直速度 -r(3); %水平加速度 1-r(4); %竖直加速度 说明物体在平面上运动时，运动的微分方程组一般可化四个一阶常微分方程，要注意速度和加速度的排列顺序。 注意 如果取r(1) = x*，r(2) = dx*/dt*，r(3) = y*，r(4) = dy*/dt*，四个一阶方程组需要修改如下 dr(1)dr(2)=r(2)=-r(2)； ， *dtdtdr(3)dr(4)=r(4)=1-r(
10、4)。 ；dt*dt*在初始时刻，小球约化位移为r(1) = r(3) = 0，初始约化速度为 dx*t0dx1dxv0*， r(2)=*=v0dts0dtmg/kdtV0r(4) = 0。 在任意时刻，r(1)和r(3)表示约化坐标，r(2)和r(4)表示约化速度。 函数文件需要修改如下(zqy6_3fun3.m) %阻力与速度成正比的平抛运动的函数 function f=fun(t,r) f= r(2); %水平速度 -r(2); %水平加速度 3 r(4); %竖直速度 1-r(4); %竖直加速度 主程序需要修改如下 t0,R=ode45(zqy6_3fun3,t,0,v0,0,0);
11、%求微分方程的数值解 plot(R(:,1),R(:,3),ko) %画轨迹 方法四：用微分方程的符号解。公式(7.1.1*)可化为二阶微分方程 d2x*dx*d2y*dy*+=0，*2-1+*=0。 (7.1.1*) dt*2dt*dtdt据此可求解微分方程的符号解。 程序zqy6_3ode.m的第四部分如下。 %阻力与速度成正比的平抛运动的轨迹(用微分方程的符号解) sx,sy=dsolve(D2x+Dx,D2y-1+Dy,Dx(0)=v0,Dy(0)=0,x(0)=0,y(0)=0)%求微分方程的符号解 x=subs(sx,v0,v0); %替换初速度 x=subs(x,t,t); %替换时间形成横坐标 y=subs(sy,v0,v0); %替换初速度 y=subs(y,t,t); %替换时间形成纵坐标 plot(x,y,ms) %画轨迹 说明符号解有两个结果 sx = v0-v0*exp(-t) sy = exp(-t)+t-1 这与手工推导的公式相同。将数值替换符号时，数值和向量要分别替换。 图示如P7.1图所示，约化速度kv0/mg取为2，与无阻力的情况相比，在有空气阻力的情况下，小球的横坐标和纵坐标都小一些。如果约化速度取其他值，也可画出类似的坐标曲线和轨迹(图略)。 P7.1图 湖南大学物电院 周群益 4
