数1402 20143291 王文鹏

MATLAB程序代码：

clear; g=9.8;

p=4*10.^(-5); % B2/m

x1(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; %修正量参数 for i=1:1:2

theta(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轴 有密度修正的初始速度 for n=1:1:1000;

V1(n)=sqrt(Vx1(n)^2+Vy1(n)^2); %合速度 ax1(n)=-(p)*V1(n)*Vx1(n); %没有修正的Fx ay1(n)=-g-(p)*V1(n)*Vy1(n); %没有修正的Fy

Vx1(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; %没有修正的x y1(n+1)=y1(n)+Vy1(n)*dt+0.5*ay1(n)*dt^2; %没有修正的y if (y1(n+1)<0) break end end

r=-y1(n)/y1(n+1);

x1(n+1)=(x1(n)+r*x1(n+1))/(r+1); y1(n+1)=0; %保证数值大于0 %以下注释与上面类似，只是有修正向 for n=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)

break end end

r=-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'); hold on

xlabel 水平距离/m ylabel ?高度/m

title('有空气阻力的抛射体运动') end

作图如下

有空气阻力的抛射体运动8000700060005000with densitycorrection高度/m40003000200010000without densitycorrection00.511.5水平距离/m22.5x 104