有空气阻力的斜抛运动MATLAB代码

数140220143291王文鹏MATLAB程序代码：clear;g=9.8;p=4*10.^(-5);%B2/mx1(1)=0;y1(1)=0;V1(1)=700;x2(1)=0;y2(1)=0;V2(1)=700;%各个值的初始量dt=0.1;theta=[35,45];%角度a=6.5*10^(-3);alfa=2.5;t=273;%修正量参数fori=1:1:2theta(i)=theta(i)*pi/180;Vx1(1)=V1(1)*cos(theta(i));%x轴没有密度修正的初始速度Vy1(1)=V1(1)*sin(theta(i));%Y轴没有密度修正的初始速度Vx2(1)=V2(1)*cos(theta(i));%x轴有密度修正的初始速度Vy2(1)=V2(1)*sin(theta(i));%Y轴有密度修正的初始速度forn=1:1:1000;V1(n)=sqrt(Vx1(n)^2+Vy1(n)^2);%合速度ax1(n)=-(p)*V1(n)*Vx1(n);%没有修正的Fxay1(n)=-g-(p)*V1(n)*Vy1(n);%没有修正的FyVx1(n+1)=Vx1(n)+ax1(n)*dt;%x轴没有密度修正的瞬时速度Vy1(n+1)=Vy1(n)+ay1(n)*dt;%y轴没有密度修正的瞬时速度x1(n+1)=x1(n)+Vx1(n)*dt+0.5*ax1(n)*dt^2;%没有修正的xy1(n+1)=y1(n)+Vy1(n)*dt+0.5*ay1(n)*dt^2;%没有修正的yif(y1(n+1)0)breakendendr=-y1(n)/y1(n+1);x1(n+1)=(x1(n)+r*x1(n+1))/(r+1);y1(n+1)=0;%保证数值大于0%以下注释与上面类似，只是有修正向forn=1:1:1000;V2(n)=sqrt(Vx2(n)^2+Vy2(n)^2);ax2(n)=-(p)*V2(n)*Vx2(n)*(1-a*y2(n)/t).^alfa;ay2(n)=-g-(p)*V2(n)*Vy2(n)*(1-a*y2(n)/t).^alfa;Vx2(n+1)=Vx2(n)+ax2(n)*dt;Vy2(n+1)=Vy2(n)+ay2(n)*dt;x2(n+1)=x2(n)+Vx2(n)*dt+0.5*ax2(n)*dt^2;y2(n+1)=y2(n)+Vy2(n)*dt+0.5*ay2(n)*dt^2;if(y2(n+1)0)breakendendr=-y2(n)/y2(n+1);x2(n+1)=(x2(n)+r*x2(n+1))/(r+1);y2(n+1)=0;plot(x1,y1,'--',x2,y2,'b');holdonxlabel水平距离/mylabel¸高度/mtitle('有空气阻力的抛射体运动')end作图如下00.511.522.5x104010002000300040005000600070008000水平距离/m高度/m有空气阻力的抛射体运动withoutdensitycorrectionwithdensitycorrection