抛体运动的计算机模拟及数值求解
山东省东阿县第一中学 王海军
引 论
探究有空气阻力时的抛体运动轨迹以及影响射程的因素。研究理想状态（忽略空气阻力）的抛体运动，可以通过初等数学，求解运动学的一元二次方程。当有空气阻力时，速度与阻力相互影响，物体受变力的作用，上述方法无法应用。能否找到一种方法，使我们仅用初等数学，并借助计算机的辅助，解决此问题。
本次探究活动用初等数学结合物理知识推导出Feynman-Newton法。运用中学所学知识建立更接近“物理真实”的物理模型。尝试运用初等数学与计算机程序语言相结合解决某一物理模型。
关键词：抛体运动 空气阻力 轨迹 射程 计算机模拟 数值求解 程序设计
正 文
空气动力学知识：
当考虑到空气阻力存在的情况下，如果空气密度为 ，作为垂直于运动方向的截面都是圆形而且最大横截面积为A的光滑落体，同样可取拖曳系数C 的值为0.46，当其速度为v 时所受到的空气阻力的量值为
Feynman-Newton 法的初等数学与物理概念相结合的理解与应用：
设A 、B为物体运动过程中相邻的极小时间间隔为h的两点，我们可以用A点的速度v求出B的速度v+h*a,A与B间的位移可以用v*h来求出。对于整个运动过程：编制程序时，在循环过程内进行迭代计算，可用赋值语句将上面的各式写为
vx + = h*ax；vy + = h*ay；
x + = h*vx；y + = h*vy；
由此可以看出，在每一次迭代计算中，实际上是将上一次迭代已经求得的加速度和速度作为 时间内的平均加速度和平均速度来计算下一时刻的速度和位置．那么，若初始条件为:t = 0 时x = 0、y = 0、 =0、 =0, 当要计算第一步末的位置x,y时, 按上面的顺序应先算出第一步末的速度再计算位移, 这就意味着要以刚计算出的第一步末的速度作为此区间的平均速度代入后两式中算位移；若将上面的位置和速度赋值语句交换顺序, 计算位置时就是以t=0 时的速度为此区间的平均速度进行运算, 这样就会导致所算出的第一步末的位移为零．显然这两种算法都不太合理, 若要想作更为精确一些的处理，应该采用新的计算方法．先计算出A、B 中间时刻的速度v=v+ h/2*a，用它算出两点间的位移h*v. 对于整个运动过程：编制程序时，vx + = h/2*ax；vy + = h/2*ay；计算第一半步长速度
在循环过程内进行迭代计算，可用赋值语句将上面的各式写为
x += h*vx；y += h*vy；
vx += h*ax；vy += h*ay；这样修改过的方法即称为Feynman-Newton 法
1．无阻力时的抛体运动
step1．根据所选取的初速度和抛射角确定显示器水平方向和竖直方向坐标取值范围．例如令X方向为[0,70]m，Y方向为[0,20] m，在640×480显示方式下，大约变量x 的值应乘以倍率8，而变量y的值应乘以倍率22．
step2．给出初始条件．例如选取 t=0, x= y= 0, 时间递进步长可选0.01s，初速度v和抛射角θ可由scanf()语句输入
step3．因X 方向无外力作用, 为恒定值, 只需计算Y 方向的速度改变． Y 方向速度的第一个半步长值
step4．建立循环．因为抛体回到地面时的Y 坐标为0, 则可检测 y值，一旦y<0，循环结束．
step5．为得到抛体的空间轨道曲线, 用连线语句lineto(20+x*8,460-y*22)；作图．
step6．更新x, y, , 在达循环终止条件前返回step4．
源程序（TC2.0编译通过）
#include
#include //定义图形函数原型的头文件
#include //定义调用dos控制台I/O例程的头文件
#include
float f1();//加速度ax的函数原型
float f2();//加速度ay的函数原型
main()
{
int driver=DETECT,mode;//初始化图形系统
float v,h=0.02,theta,x,y,vx,vy;
x=0;y=0;//赋位置初值
printf("Please input v0=");scanf("%f",&v);//按提示输入初速度
printf("Please input theta=");scanf("%f",&theta);//按提示输入抛射角
vx=v*cos(theta*M_PI/180);vy=v*sin(theta*M_PI/180);
initgraph(&driver,&mode," ");setbkcolor(WHITE);setcolor(BLUE);//初始化图形系统
rectangle(20,20,620,460);moveto(20,460);//沿边界画线框，将图像点移到原点
vy+=h*f2()/2;//计算第一半步长Y方向速度
do{
lineto(20+x*8,460-y*22);//画无空气阻力时的抛体运动曲线
x+=h*vx;y+=h*vy;vx+=h*f1();vy+=h*f2();//计算t时刻的位移、速度
}while(y>=0);//当y<0终止循环
getch();
return (0);
}
float f1()//加速度ax的函数原型
{
float ax;ax=0;
return (ax);
}
float f2()//加速度ay的函数原型
{
float ay;ay=-9.81;
return (ay);
}
2有空气阻力时的抛体运动
当考虑到空气阻力存在的情况下，如果空气密度为 ，作为垂直于运动方向的截面都是圆形而且最大横截面积为A的光滑落体，同样可取拖曳系数C 的值为0.46，当其速度为v 时所受到的空气阻力的量值为
设在运动过程中任一位置时, 抛体速度v 与X 轴向间交角为φ, 则沿X 和Y 方向的运动方程分别为
取参数
则沿两个坐标轴方向的加速度分量可表示为
加速度分量的函数形式应分别改为
ax = -k*sqr(vx*vx+vy*vy)*vx；
ay = -9.81-k*sqr(vx*vx+vy*vy)*vy；
还要给出运动物体的质量、密度、截面半径和拖曳系数，然后求得参量k 的值
k=0.5*c*pi*r*r*rho/m；
有空气阻力和无阻力时的抛体运动轨道的比较
源程序（TC2.0编译通过）
#include
#include
#include
int main(void)
{
int gdriver=DETECT,gmode;//定义初始化图形系统用到的整型变量
float c=0.46,r=0.0732/2,m=0.145,rho=1.2,h=0.01,k,t,v,theta,x,y,vx,vy;
float f1(float k,float vx,float vy);//定义加速度ax的函数原型，有参数传递
float f2(float k,float vx,float vy);//定义加速度ay的函数原型，有参数传递
initgraph(&gdriver,&gmode," ");
setbkcolor(WHITE);setcolor(BLUE);//初始化图形系统
t=0;x=0;y=0;v=30;theta=40*M_PI/180;k=0.5*c*M_PI*rho*r*r/m;
vx=v*cos(theta);vy=v*sin(theta);
rectangle(20,20,620,460);moveto(20,460);//沿边界画线框，将图像点移到原点
do
{
lineto(20+vx*t*6,460-(vy*t-0.5*9.81*t*t)*22);//无阻力运动轨迹图
t+=h;
}
while(vy*t-0.5*9.81*t*t>=0);//当y<0终止循环
vx+=h*f1(k,vx,vy)/2;vy+=h*f2(k,vx,vy)/2;//计算第一半步长速度
moveto(20,460);
do{
lineto(20+x*6,460-y*22);//画有空气阻力时运动曲线
x+=h*vx;y+=h*vy;//计算t时刻的位移
vx+=h*f1(k,vx,vy);vy+=h*f2(k,vx,vy);//计算t时刻的速度
}while(y>=0);//当y<0终止循环
getch();
return 0;
}
float f1(float k,float vx,float vy)//加速度ax的函数原型
{
float ax;
ax=-k*sqrt(vx*vx+vy*vy)*vx;
return (ax);
}
float f2(float k,float vx,float vy)//加速度ay的函数原型
{
float ay;
ay=-9.81-k*sqrt(vx*vx+vy*vy)*vy;
return (ay);
}
不同初速度抛体的运动轨道曲线簇
源程序（TC2.0编译通过）
#include
#include
#include
int main(void)
{
int gdriver=DETECT,gmode,i;//定义初始化图形系统用到的整型变量
float c=0.46,r=0.0732/2,m=0.145,rho=1.2,h=0.01,theta=40*M_PI/180,k,v,x,y,vx,vy;
float f1(float k,float vx,float vy);//定义加速度ax的函数原型，有参数传递
float f2(float k,float vx,float vy);//定义加速度ay的函数原型，有参数传递
initgraph(&gdriver,&gmode," ");
setbkcolor(WHITE);setcolor(BLUE);//初始化图形系统
rectangle(20,20,620,460);
v=20;k=0.5*c*M_PI*rho*r*r/m;
for(i=0;i<=6;i++)//建立循环，随参数i逐次改变初速度的值
{
x=0;y=0;vx=(v+4*i)*cos(theta);vy=(v+4*i)*sin(theta);
vx+=h*f1(k,vx,vy)/2;vy+=h*f2(k,vx,vy)/2;//计算第一个半步长的速度
moveto(20,460);
do
{
lineto(20+x*6,460-y*12);
x+=h*vx;y+=h*vy;//计算t时刻的位移
vx+=h*f1(k,vx,vy);vy+=h*f2(k,vx,vy);//计算t时刻的速度
}
while(y>=0);//当y<0终止循环
}
getch();
return 0;
}
float f1(float k,float vx,float vy)
{
float ax;
ax=-k*sqrt(vx*vx+vy*vy)*vx;
return (ax);
}
float f2(float k,float vx,float vy)
{
float ay;
ay=-9.81-k*sqrt(vx*vx+vy*vy)*vy;
return (ay);
}
1．无阻力时的抛体运动，由运动学方程代入求根公式求解，显示固定数值的抛射高度、抛射速度条件下能达到的最大水平位移及对应的抛射角。
源程序（TC3.0编译通过）
#include
#include
#include
#define PI 3.1415926
int main()
{
float v,g=9.8,h;//定义速度、高度、重力加速度并赋值
float m,n;//定义存储最大水平位移的变量和对应此位移的抛射角
float x,a;//定义位移及抛射角
for (h=1.4;h<=2;h+=0.05)//让抛射高度递增
{
for (v=1;v<=12;v++)//让抛射速度递增
{ m=0;n=0;
for (a=0;a<=PI/2;)//让抛射角度从0递增到PI/2
{
x=v*v*cos(a)*sin(a)/g+v*cos(a)*sqrt(v*v*sin(a)*sin(a)+2*g*h)/g;//由运动方程求解
if (x>=m)
{
m=x;
n=a*180/PI;//判断并记录最大水平位移及对应的抛射角
}
a+=0.01;//让抛射角递增
}
printf("h=%3.2f\n",h);
printf("v=%2.0f\n",v);
printf("x=%5.2f,a=%5.2f\n\n",m,n);//显示固定数值的抛射高度、抛射速度条件下能达到的最大水平位移及对应的抛射角
}
}
return 0;
}
2．无空气阻力时的抛体运动，利用Feynman-Newton法，显示固定数值的抛射高度、抛射速度条件下能达到的最大水平位移及对应的抛射角。
源程序（TC3.0编译通过）
#include
#include
#include
#define PI 3.1415926
int main()
{
float h,v;//定义速度、高度
float t=0.01;//时间递增的步长
float vx,vy;//定义两个方向的速度
float x,y;//两个方向的位移
float m,n;//定义存储最大水平位移的变量和对应此位移的抛射角
float f1();
float f2();
for (h=1.4;h<=2;h+=0.05)//让抛射高度递增
{
for (v=1;v<=12;v++)//让抛射速度递增
{
m=0;n=0;
for (a=0;a<=PI/2;)//让抛射角度从0递增到PI/2
{
x=0;y=h;
vx=v*cos(a);vy=v*sin(a);vy+=t*f2()/2;
do
{
x+=t*vx;y+=t*vy;vy+=t*f2();//利用时间递增逐步求解
}while(y>=0);
if (x>=m)
{
m=x;
n=a*180/PI;//判断并记录最大水平位移及对应的抛射角
}
a+=0.01;//让抛射角递增
}
printf("h=%3.2f\n",h);
printf("v=%2.0f\n",v);
printf("x=%5.2f,a=%5.2f\n\n",m,n);//显示固定数值的抛射高度、抛射速度条件下能达到的最大水平位移及对应的抛射角
}
}
return 0;
}
float f1()//定义加速度ax的函数原型，无参数传递
{
float ax;ax=0;
return (ax);
}
float f2()
{
float ay;ay=-9.81;//定义加速度ay的函数原型，无参数传递
return (ay);
}
3. 有空气阻力时的抛体运动，利用Feynman-Newton法，显示固定数值的抛射高度、抛射速度条件下能达到的最大水平位移及对应的抛射角。
源程序（TC3.0编译通过）
#include
#include
#include
#define PI 3.1415926
int main()
{
float h,v,a;//定义速度、高度、加速度
float t=0.01;//时间递增的步长
float vx,vy;//定义两个方向的速度
float x,y;//两个方向的位移
float m,n;//定义存储最大水平位移的变量和对应此位移的抛射角
float c=0.46,r=0.05,z=4,rho=1.2;
float k=0.5*c*PI*rho*r*r/z;//与空气阻力有关的几个量
float f1(float k,float vx,float vy);
float f2(float k,float vx,float vy);
for (h=1.4;h<=2;h+=0.05)
{
for (v=1;v<=12;v++)
{
m=0;n=0;
for (a=0;a<=PI/2;)
{